A new model for the Xray continuum of the magnetized accreting pulsars
^{1}
INAF–Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio
5,
35122
Padova,
Italy
email:
ruben.farinelli@oapd.inaf.it
^{2}
ISDC, Chemin d’Écogia 16, 1290
Versoix,
Switzerland
^{3}
Department of Physics and Astronomy, George Mason
University, Fairfax,
VA
22030,
USA
Received: 26 August 2015
Accepted: 27 January 2016
Context. Accreting highly magnetized pulsars in binary systems are among the brightest Xray emitters in our Galaxy. Although a number of highquality broadband (0.1–100 keV) Xray observations are available, the spectral energy distribution of these sources is usually investigated by adopting pure phenomenological models rather than models linked to the physics of accretion.
Aims. In this paper, a detailed spectral study of the Xray emission recorded from the highmass Xray binary pulsars Cen X3, 4U 0115+63, and Her X1 is carried out by using BeppoSAX and joined Suzaku +NuStar data, together with an advanced version of the compmag model, which provides a physical description of the highenergy emission from accreting pulsars, including the thermal and bulk Comptonization of cyclotron and bremsstrahlung seed photons along the neutron star accretion column.
Methods. The compmag model is based on an iterative method for solving secondorder partial differential equations, whose convergence algorithm has been improved and consolidated during the preparation of this paper.
Results. Our analysis shows that the broadband Xray continuum of all considered sources can be selfconsistently described by the compmag model. The cyclotron absorption features (not included in the model) can be accounted for by using Gaussian components. From the fits of the compmag model to the data we inferred the physical properties of the accretion columns in all sources, finding values reasonably close to those theoretically expected according to our current understanding of accretion in highly magnetized neutron stars.
Conclusions. The updated version of the compmag model has been tailored to the physical processes that are known to occur in the columns of highly magnetized accreting neutron stars and it can thus provide a better understanding of the highenergy radiation from these sources. The availability of broadband highquality Xray data, such as those provided by BeppoSAX in the past and currently from NuStar and other facilities, is crucial to fully exploit the potentialities of the model. The advent of the AstroH mission, endowed with an unprecedented combination of high sensitivity and Xray broadband coverage, provides good perspectives to improve our understanding of accretion onto highly magnetized neutron stars through physical models like the one adopted here.
Key words: Xrays: binaries / magnetic fields / radiative transfer / pulsars: individual: Cen X3 / pulsars: individual: 4U 0115+63 / pulsars: individual: Her X1
© ESO, 2016
1. Introduction
Xray binary pulsars (XRBPs) were discovered more than forty years ago with the pioneering observations of the bright sources Her X1 (Giacconi et al. 1971) and Cen X3 (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^{11−12} 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 highmass Xray binaries (HMXB; see e.g. Walter et al. 2015, for a recent review). Some binary systems with intermediate properties between HMXBs and lowmass Xray binaries (LMXBs) also show evidence of extended accretion columns (IMXBs; see e.g. the case of Her X1; 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 Xrays as the plasma decelerates and settles onto the stellar surface. Xray pulsations are generated owing to the disalignment between the NS magnetic and rotational axes. A peculiar feature that is observed in the Xray spectra of many HMXBs and IMXBs is the cyclotron resonant scattering absorption line (CRSF), which provides a direct measurement of the NS magneticfield strength. This can be estimated by using the relation E_{cyc} = 11.6 B_{12} × (1 + z)^{1} keV, where E_{cyc} is the centroid energy of the fundamental CRSF, B_{12} is the NS magnetic field strength in units of 10^{12} G, and z is the gravitational redshift in the lineforming 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 broadband Xray emission of accreting Xray pulsars is usually described by using phenomenological models, including an absorbed power law extending up to ~100 keV with a rollover 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 highquality observations (Ferrigno et al. 2009, hereafter F09). The development of refined spectral models linking the emission from Xray 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 angulardependent 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 freefalling 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 Xray pulsars (e.g. the supergiant fast Xray 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/Xray 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 X1, Cen X3, and 4U 0115+63, which are known to possess prominent cyclotron lines in their Xray spectra. The Xray spectrum of 4U 0115+63 is also known to be characterized by the presence of a very strong and broad emissionlike feature at 10 keV, i.e. right below the energy of the fundamental cyclotron absorption line (F09). This “10keV feature” (Coburn et al. 2002) has been observed in a number of Xray 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 10keV 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 Xray pulsars Her X1, Cen X3, 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 secondorder bulkComptonization 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 FokkerPlanck approximation of the radiative transfer equation (RTE) in the accretion column is based on a general finitedifference 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 Xray transients as these sources host pulsars which usually display a relatively low Xray luminosity (≃10^{33−34} erg s^{1}; Farinelli et al. 2012b). As shown by BW07, in brighter accreting Xray 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 Xray 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 integrodifferential Boltzmann equation with a particle diffusion term and a Coulomb scattering kernel (Riffert et al. 1999). The protonelectron 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_{e}v^{2} 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 electronproton 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 nonmagnetized bremsstrahlung emissivity spectrum. However, a fully selfconsistent treatment is particularly difficult and barely feasible in a runtime 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 phaseaveraged 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 secondorder bulk Comptonization effects are taken into account, a spacedependent term containing a β^{2} factor appears in the energy operator. In this case, the variableseparation 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^{12} 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 electronscattering cross section for Ophotons is (1)where (2)The approximate cross section of the Ephotons 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 Ophotons and Ephotons 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 transitionrate term should thus be solved. As outlined by BW07, the angular dependence of the cross sections must be replaced by some angleaverage 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 energyaveraged cross sections. On the one hand, Ephotons 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, Ophotons 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 Ophotons, while the contribution of Ephotons is phenomenologically taken into account by using Gaussian absorption features in the Xray spectra.
The RTE for the zerothmoment photon occupation number n(x) in the cylindrically symmetric case can be written as
where x = E/kT_{e}, Θ = kT_{e}/m_{e}c^{2}, f_{b} = 1 + m_{e}v^{2}/ 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 secondorder 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 freefree 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^{2}/c^{2} and can thus be neglected (see Eq. (B4) in PL97). The freefree 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_{0}(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 Ophotons α_{ff} is a smooth function of energy which can be closely approximated by that of unmagnetized case. For Ephotons, 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 Ophotons, the expression for the absorption coefficient in Eq. (6) is a good approximation, while the freefree absorption resonances of Ephotons 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 Ephotons mainly interact with the electrons via the cyclotron resonance, but this is effectively an elastic scattering process that preserves the angleaveraged 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_{0}τ_{⊥}/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^{2}/c^{2} (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 zerothorder term for the emission coefficient.
Using the definition dτ = n_{e}σ_{∥}dz for the optical depth along the zaxis, performing a logarithmic sampling of the adimensional energy through the change of variable x = e^{q}, and defining n = J/x^{3}, after some lengthy calculations we obtain the adimensional equation for the zero moment intensity, (9)where , ξ = 15.5r_{0}/ṁ, and , with H = (σ/σ_{∥})Θ.
The parameter r_{0} 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_{12} = B/ 10^{12}. 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 freefree absorption coefficient and are defined as (13)and (14)where c_{1} = 2.1 × 10^{13} and c_{2} = 4c_{1}.
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_{0}). 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 Xray 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 β_{0} is the matter velocity in units of c at the NS surface, and z_{0} is the adimensional NS radius. Given the definition of r_{0}, we obtain km. Here and throughout the paper we assume m = 1.4 and R_{ns} = 10 km, implying z_{0} = 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 τ_{0} 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, τ_{0} 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 conelike shape of halfangle θ_{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 Xray emission. We now differentiate between the pencil beam case, when the Xray 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 topcolumn surface , and divided by 4πD^{2}, where D is the source distance. Using Eqs. (22) and (23), we obtain (24)where r_{0} and d_{10} are the adimensional column radius and source distance in units of 10 kpc, respectively.
In the fanbeam 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 zerothmoment specific intensity J(x,τ) = x^{3}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_{0} 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_{0}, to be inferred from the fit to the Xray data. This means that, in principle, N_{comp} should be kept fixed during the fit to a value of the order of ~, where d_{10} 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 spacetime. In a more realistic situation, it is likely that hydrodynamical and General Relativity (GR) effects will affect the spectral energy distribution of the Xray radiation. In particular, the gravitationallensing effect – due to the curved spacetime close to the NS surface – is expected to enhance the emerging Xray 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 lowquality data that were available to F12. The data that are used in the present paper to study some of the brightest Xray 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 twopoint linear interpolation largely uncertain, resulting in the nonconvergence 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 predefined 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 χ^{2} 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 quadcore Linux machine with a 2.4 GHz processor).
3. Sources
3.1. Cen X3
The highmass Xray binary Cen X3 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 68 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 Xray 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 Otype star, providing convincing indications for a strong perturbation of the wind by the NS Xray radiation.
The source Xray spectrum is generally well described by a powerlaw model with exponential cutoff and a Gaussianlike 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 hourlong dips in the source lightcurve (see e.g. Naik et al. 2011). The Xray variability of Cen X3 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 Xray radiation from 4U 0115+63 was first discovered in 1969 (Johns et al. 1978) during one of the gianttype II outbursts displayed by the source. On that occasion, the peak luminosity of the event reached a few 10^{37} 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_{X}sini = 140.1 lt−s. Xray 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 Xray spectrum of the source above ~10 keV is usually well described by powerlaw model with an exponential cutoff (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 broadband Xray spectrum of 4U 0115+63 with the BW07 model, aiming at constraining the properties of the Xray emitting region and the highenergy 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 X1
Her X1 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 Xray eclipses, lasting about 20% of the binary orbit (Scott et al. 2000). The superorbital modulation of the source Xray 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 NSdisk system (Staubert et al. 2013). During the superorbital period, Her X1 displays a “mainon” state lasting 10–11 days where the source achieves the highest Xray luminosity and a short slightly dimmer state which typically last 5 to 7 days. At other superorbital phases, the source is still visible in the Xray domain, although at a much lower luminosity level.
Her X1 is characterized by low Galactic absorption and is situated at a distance of 6.6 ± 0.4 kpc (Reynolds et al. 1997). The source Xray spectrum shows emission lines produced from the photoionized accretion disk and on the companion surface, as shown by Chandra and XMMNewton grating observations (see e.g. Zane et al. 2004; JimenezGarate 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 Xray 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 socalled 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 X1 where the cyclotron line is typically monitored, these authors showed that the source should be in the supercritical 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 X1 is strongly nondipolar, resulting in multiple accretion sites on the NS surface with noncylindrical sections (Shakura et al. 1991). It is worth noting that the combination of a pencil and a fanbeamed emission from a misaligned dipole would also be able to reproduce the pulse profiles of Her X1 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 Xray 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 X1 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 mainon 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 highquality 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 X3 and 4U 0115+63, together with a combined NuSTAR (Harrison et al. 2013) and Suzaku (Mitsuda et al. 2007) observation of Her X1 during one of the source mainon 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 X3 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.
Log of the observations of the sources analysed in this paper.
For the BeppoSAX observations, we analysed data from the highenergy 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 offsource 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 signaltonoise ratio.
We adopted standard reduction procedures for the Suzaku and NuSTAR data, devoting particular attention to evaluating the effect of pileup 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 phaseresolved spectral analysis.
All spectral analyses were performed with the Xspec version 12.8 (Arnaud 1996), using the compmag model described in Sect. 2. Intercalibration constants were included in the fits and found to be compatible with typically expected values^{1}. 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 X3 and Her X1, 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 Xray pulsars emerging in both a pencilbeam or a fanbeam 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 Xray flux. The parameter β_{0} defines in this geometry the terminal velocity of the infalling material at the NS surface. The corresponding equations for the fanbeam geometry are Eqs. (18) and (26).
A preliminary fit to the available spectra considered in this work revealed that reasonably good values of the χ^{2} 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 fanbeam emission geometry is expected to set in, we first fixed the normalization of the compmag model (N_{comp}) and concentrate on the fanbeam 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 pencilbeam geometry.
In all cases, we modelled the effect of the interstellar photoelectric absorption on the Xray emission from each of the considered sources with the phabs spectral component in Xspec, assuming solar abundances.
5.1. Cen X3
We first fit the low luminosity (LL) and high luminosity (HL) spectra of Cen X3 obtained from the two available BeppoSAX observations (see Table 1) with the compmag model in the fanbeam geometry case, plus a cyclotron Gaussian absorption line centred at ~30 keV. The results of the fits were not yet acceptable (χ^{2}/d.o.f. = 4049/318 and χ^{2}/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 (χ^{2}/d.o.f. = 509/315) and HL (χ^{2}/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 χ^{2}/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 phaseresolved 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^{22} to ~1.3 × 10^{23} 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 (χ^{2}/d.o.f. = 369/311) with N_{H} ~ 5 × 10^{22} 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 X3. Also show are the bestfit compmag model and the residuals from the fit in units of σ (bottom panel). All bestfit parameters determined from the fit are given in Table 2. 

Open with DEXTER 
Although the model discussed above provided a reasonably good fit to both the available observations of Cen X3, we obtained a similar accretion rate in both observations from the compmag bestfit parameters. This result does not seem to be consistent with the eight times higher Xray 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 χ^{2}/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 phaseaveraged 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 Xray 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 bestfit 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 pencilbeam geometry (assuming a freefall velocity profile for the material within the accretion column), we have to introduce the flow terminal velocity β_{0} 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 χ^{2}. 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_{0}, i.e. and (as expected according to the parameter definitions; see Sect. 2). We measured the slopes b_{1} = 1.12 ± 0.03 and b_{2} = 0.44 ± 0.01 from the fits to the LL observation, b_{1} = 0.94 ± 0.02 and b_{2} = 0.49 ± 0.04 from the HL data. By assuming for the pencilbeam geometry the same value of R_{0} measured for the fanbeam 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 pencilbeam 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 fanbeam 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}< 0.1 and z_{max} ≲ 4 km for either the LL and HL data.
5.2. 4U 0115+63
Bestfit models for the sample of sources analysed in this paper.
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 fanbeam geometry, plus two absorption Gaussian lines at energies corresponding to the fundamental and the second harmonic of the CRSF. This gave χ^{2}/d.o.f. = 710/422 and χ^{2}/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 χ^{2}/d.o.f. = 518/420 and χ^{2}/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 χ^{2}/d.o.f. = 444/418 for Obs. I and χ^{2}/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 χ^{2}/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 Xray 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. 

Open with DEXTER 
Concerning the applicability of the compmag model in the pencilbeam geometry to 4U 0115+63 (assuming a freefall velocity for the material inside the accretion column), the same considerations of Cen X3 apply (see Sect. 5.1). From the fit to the source spectrum in Obs. I, we determined b_{1} = 1.06 ± 0.02 and b_{2} = 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_{1} = 0.93 ± 0.02 and b_{2} = 0.48 ± 0.01. The corresponding values of the mass accretion rate for a column radius R_{0} equal to the value obtained in the fanbeam 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 X3, the flow terminal velocity and column height determined for the two 4U 0115+63 observations remain roughly constant, with β_{0} ~ 0.2 and z_{max} ~ 3 km.
5.3. Her X1
A fit with an absorbed compmag model, together with a cyclotron absorption line, to the combined Suzaku and NuStar observation of Her X1 provided an unacceptable result (χ^{2}/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 (χ^{2}/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 cutoff powerlaw model, and thus the second line is likely to have a modeldependent rather than a physical origin.
Fig. 3 Same as Fig. 1 but for the combined Suzaku and NuStar observation of Her X1. 

Open with DEXTER 
Fig. 4 From top to bottom: comptonized bremsstrahlung and cyclotron spectra of Cen X3, 4U 0115+63, and Her X1 obtained from the bestfit models given in Table 2. 

Open with DEXTER 
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 (χ^{2}/d.o.f. = 3769/3727). The estimated contributions of the cyclotron and bremsstrahlung Comptonized component to the total Xray emission of the source in the 3–100 keV energy range were 60% and 40%, respectively. We report all values of the bestfit model to the data of Her X1 in Table 2, and show the unfolded spectrum with the residuals from the fit in Fig. 3. The high statistical quality of the Her X1 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 X1, 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 highenergy cutoff (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 X1 it was not possible to find a formally acceptable fit to the data by using the compmag model in the pencilbeam geometry and considering a freefall 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 Xray luminosity emerging from a highly magnetized accreting pulsar in the fanbeam and pencilbeam 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 Xray 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 Xray 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 Xray 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 Xray 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 broadband Xray observations of three bright Xray pulsars: 4U 0115+63, Cen X3, and Her X1. We used publicly available BeppoSAX archival data for Cen X3 and 4U 0115+63, and a recent combined Suzaku +NuStar observation of Her X1. In all cases, the broadband 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 selfconsistent 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 Xray 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 Xray radiation arising from the pulsar accretion column. A similar conclusion applies to the case of Cen X3, 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 Xray 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 X1, 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 broadband Xray 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 Xray 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 FokkerPlanck 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^{2} 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 Xray 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_{0} and its height H, would probably need to be scaled in a nontrivial 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 Xray astrophysical community. A wellknown 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 Xray 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 X1, and Cen X3, provides a value of R_{0} ~ 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^{12} 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 electrondensity 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^{3} of the BW07 model to a broadband spectrum of an Xray 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 Xray 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_{0} had to be fixed in the fit (6 × 10^{16} 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 broadband 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 firstorder bulk Comptonization term in the energydiffusion operator, while we also introduced the secondorder term. This component, being spacedependent, cannot be included when performing the analytical treatment of Comptonization through the variableseparation 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 Xray 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 Xray continuum and thus the detection of particularly weak wiggles in the highenergy end of an exponentially decaying spectrum can easily turn out to be statistically nonsignificant 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 selfconsistent 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 Xray 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 electronproton 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 protonelectron bremsstrahlung, there is a second channel for the production of cyclotron emission photons due to the photonelectron 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 Xray data (see Fig. 4) would thus support the idea that the eγ RS dominates over the ep 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 ep processes dominate over the eγ interactions. A more realistic solution to the problem would require a selfconsistent treatment where both processes are taken into account simultaneously via a nonlinear 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 magnetohydrodynamical effects. First, we emphasize that the cyclotron emission term in the compmag model (Eq. (10)) is selfconsistently weighted over the vertical profile of the accretion column because we assume B_{z} = B_{0}(z_{0}/z)^{3}. 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^{2} = x^{2} + y^{2}. 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.
Magnetohydrodynamical simulations of hypercritical 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 protonelectron 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 pe 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_{e}v^{2}, E′ = 1 / 2m_{e}v^{′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^{13} G. The pe 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 pe 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 zeromagneticfield 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 pe 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 rollover of the spectra, the correlation we found is most likely related to the scaling law between the cutoff 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 Xray 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 X1, and Cen X3 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 Xray pulsars, the possible correlations we identified in Sect. 5 between the different geometrical parameters of the compmag model (i.e. R_{0} and z_{max}) should instead be taken with caution. These relations were derived by assuming a purely Newtonian treatment of the Xray 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 raytracing 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 Xray 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). 

Open with DEXTER 
8. Geometry of the emission region
In modelling the source spectra, we have considered both a pencilbeam geometry with an accelerating velocity profile in the accretion column and a fanbeam geometry where matter decelerates towards the NS surface. The first case is expected to occur in less luminous sources (L_{X}< 10^{37} 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 (fanbeam). At higher luminosities (L_{X}> 10^{37} 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 electronphoton 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 (fanbeam), 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 phaseaveraged 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 Xray luminosities L_{X}> 10^{37} erg s^{1}, and thus the usage of the fanbeam geometry in the compmag model was considered more appropriate to fit their Xray spectra. We have shown in principle that fits with the compmag model and the pencilbeam geometry are possible for Cen X3 and 4U 0115+63, while no solution could be found in the case of Her X1. In the pencilbeam case, the empirical dependence discussed in Sects. 5.1 and 5.2 is derived from the definition of the compmag normalization for the topcolumn 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 Xray flux, a decrease of N_{comp} has to be compensated by a corresponding linear increase in ṁ. We verified that the changes in ṁ and R_{0} derived as a consequence of a variation in the model normalization, N_{comp}, are sufficient to keep the value of χ^{2}unchanged. As we reported during the analysis of Cen X3 and 4U 0115+63, both β_{0} and z_{max} are insensitive to variations of N_{comp}. This can be explained considering that the emerging flux (and thus the source Xray 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 pencilbeam case than in the fanbeam 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 fanbeam or the pencilbeam geometry is better suited to describing the Xray spectrum of an Xray pulsar. On the other hand, the fact that the pencilbeam geometry could not be used to satisfactorily fit the Xray spectrum of Her X1 gives us confidence that the compmag model might be able to distinguish between the different emission states of these sources, at least when highquality and wide broadband data are available. Additional observations of these sources with the next generation of Xray instruments providing broadband Xray coverage (like those onboard AstroH; 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 Xray pulsars for which broadband highquality data were available. We found that the model is well suited to describing quantitatively the Xray spectral energy distribution of these systems and that it provides reasonable values of the physical and geometrical parameters that are used to compute the highenergy 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 photonelectron 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 phaseaveraged spectra. Despite its limitations, the compmag model represents an additional significant step forward in the study of highly magnetized accreting Xray pulsars.
See https://heasarc.gsfc.nasa.gov/docs/sax/abc/saxabc/saxabc.html
An application of the BW07 model to other sources was first discussed in a conference proceedings by Marcu et al. (2015).
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 X1 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 Xray astronomy satellite BeppoSAX, a project of the Italian Space Agency (ASI) with participation of the Netherlands Agency for Aerospace Programs (NIVR).
References
 Arnaud, K. A. 1996, Astronomical Data Analysis Software and Systems V, 101, 17 [NASA ADS] [Google Scholar]
 Arons, J., Klein, R. I., & Lea, S. M. 1987, ApJ, 312, 666 [NASA ADS] [CrossRef] [Google Scholar]
 Ash, T. D. C., Reynolds, A. P., Roche, P., et al. 1999, MNRAS, 307, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Basko, M. M., & Sunyaev, R. A. 1976, MNRAS, 175, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, P. A., & Wolff, M. T. 2007, ApJ, 654, 435 (BW07) [NASA ADS] [CrossRef] [Google Scholar]
 Becker, P. A., Klochkov, D., Schönherr, G., et al. 2012, A&A, 544, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bernal, C. G., Lee, W. H., & Page, D. 2010, Rev. Mex. Astron. Astrofis., 46, 309 [NASA ADS] [Google Scholar]
 Bildsten, L., Chakrabarty, D., Chiu, J., et al. 1997, ApJS, 113, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1981, MNRAS, 194, 1033 (BP81) [NASA ADS] [Google Scholar]
 Blum, S., & Kraus, U. 2000, ApJ, 529, 968 [NASA ADS] [CrossRef] [Google Scholar]
 Boella, G., Butler, R. C., Perola, G. C., et al. 1997a, A&AS, 122, 299 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boella, G., Chiappetti, L., Conti, G., et al. 1997b, A&AS, 122, 327 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Burderi, L., di Salvo, T., Robba, N. R., La Barbera, A., & Guainazzi, M. 2000, ApJ, 530, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Caballero, I., Kretschmar, P., Santangelo, A., et al. 2007, A&A, 465, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Canuto, V., Lodenquai, J., & Ruderman, M. 1971, Phys. Rev. D, 3, 2303 [NASA ADS] [CrossRef] [Google Scholar]
 Coburn, W., Heindl, W. A., Rothschild, R. E., et al. 2002, ApJ, 580, 394 [NASA ADS] [CrossRef] [Google Scholar]
 Cominsky, L., Clark, G. W., Li, F., Mayer, W., & Rappaport, S. 1978, Nature, 273, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Daugherty, J. K., & Harding, A. K. 1986, ApJ, 309, 362 [NASA ADS] [CrossRef] [Google Scholar]
 Davidson, K. 1973, Nat. Phys. Sci., 246, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Davidson, K., & Ostriker, J. P. 1973, ApJ, 179, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Enoto, T., Makishima, K., Terada, Y., et al. 2008, PASJ, 60, 57 [NASA ADS] [Google Scholar]
 Farinelli, R., Titarchuk, L., Paizis, A., & Frontera, F. 2008, ApJ, 680, 602 [NASA ADS] [CrossRef] [Google Scholar]
 Farinelli, R., Ceccobello, C., Romano, P., & Titarchuk, L. 2012a, A&A, 538, A67 (F12) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Farinelli, R., Romano, P., Mangano, V., et al. 2012b, MNRAS, 424, 2854 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrigno, C., Becker, P. A., Segreto, A., Mineo, T., & Santangelo, A. 2009, A&A, 498, 825 (F09) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics (Cambridge, UK: Cambridge University Press), 398 [Google Scholar]
 Frontera, F., Costa, E., dal Fiume, D., et al. 1997, A&AS, 122, 357 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fürst, F., Grefenstette, B. W., Staubert, R., et al. 2013, ApJ, 779, 69 (F13) [NASA ADS] [CrossRef] [Google Scholar]
 Giacconi, R., Gursky, H., Kellogg, E., Schreier, E., & Tananbaum, H. 1971, ApJ, 167, L67 [NASA ADS] [CrossRef] [Google Scholar]
 Harding, A. K., & Daugherty, J. K. 1991, ApJ, 374, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Heindl, W. A., Coburn, W., Gruber, D. E., et al. 1999, ApJ, 521, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Herold, H. 1979, Phys. Rev. D, 19, 2868 [NASA ADS] [CrossRef] [Google Scholar]
 Hutchings, J. B., Cowley, A. P., Crampton, D., van Paradijs, J., & White, N. E. 1979, ApJ, 229, 1079 [NASA ADS] [CrossRef] [Google Scholar]
 Iaria, R., di Salvo, T., Robba, N. R., et al. 2005, ApJ, 634, L161 [NASA ADS] [CrossRef] [Google Scholar]
 Ji, L., Schulz, N., Nowak, M., Marshall, H. L., & Kallman, T. 2009, ApJ, 700, 977 [NASA ADS] [CrossRef] [Google Scholar]
 JimenezGarate, M. A., Raymond, J. C., Liedahl, D. A., & Hailey, C. J. 2005, ApJ, 625, 931 [NASA ADS] [CrossRef] [Google Scholar]
 Johns, M., Koski, A., Canizares, C., et al. 1978, IAU Circ., 3171 [Google Scholar]
 Klochkov, D., Staubert, R., Santangelo, A., Rothschild, R. E., & Ferrigno, C. 2011, A&A, 532, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Klochkov, D., Suleimanov, V., Pühlhofer, G., et al. 2015, A&A, 573, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krzeminski, W. 1974, ApJ, 192, L135 [NASA ADS] [CrossRef] [Google Scholar]
 Leahy, D. A. 2004, MNRAS, 348, 932 [NASA ADS] [CrossRef] [Google Scholar]
 Lyubarskii, Y. E., & Syunyaev, R. A. 1982, Sov. Astron. Lett., 8, 330 [NASA ADS] [Google Scholar]
 Makishima, K., Mihara, T., Nagase, F., & Tanaka, Y. 1999, ApJ, 525, 978 [NASA ADS] [CrossRef] [Google Scholar]
 Manzo, G., Giarrusso, S., Santangelo, A., et al. 1997, A&AS, 122, 341 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marcu, D. M., Pottschmidt, K., Gottlieb, A. M., et al. 2015, ArXiv eprints [arXiv:1502.03437] [Google Scholar]
 Mihara, T., Makishima, K., & Nagase, F. 2004, ApJ, 610, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Mitsuda, K., Inoue, H., Koyama, K., et al. 1984, PASJ, 36, 741 [NASA ADS] [Google Scholar]
 Mitsuda, K., Bautz, M., Inoue, H., et al. 2007, PASJ, 59, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mukherjee, D., Bhattacharya, D., & Mignone, A. 2013a, MNRAS, 430, 1976 [NASA ADS] [CrossRef] [Google Scholar]
 Mukherjee, D., Bhattacharya, D., & Mignone, A. 2013b, MNRAS, 435, 718 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, S., Ferrigno, C., Kühnel, M., et al. 2013, A&A, 551, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Poutanen, J. 2015, MNRAS, 447, 1847 [NASA ADS] [CrossRef] [Google Scholar]
 Nagel, W. 1980, ApJ, 236, 904 [NASA ADS] [CrossRef] [Google Scholar]
 Naik, S., Paul, B., & Ali, Z. 2011, ApJ, 737, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Nakajima, M., Mihara, T., Makishima, K., & Niko, H. 2006, ApJ, 646, 1125 [NASA ADS] [CrossRef] [Google Scholar]
 Negueruela, I., & Okazaki, A. T. 2001, A&A, 369, 108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neugebauer, P., Riffert, H., Herold, H., & Ruder, H. 1996, Phys. Rev. A, 54, 467 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Odaka, H., Khangulyan, D., Tanaka, Y. T., et al. 2013, ApJ, 767, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Odaka, H., Khangulyan, D., Tanaka, Y. T., et al. 2014, ApJ, 780, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Orlandini, M. 2004, Proc. Fourth AGILE Science Workshop, eds. M. Tavani, A. Pellizzoni, & S. Vercellone, Xray and Gammaray Astrophysics of Galactic Sources, 119 [Google Scholar]
 Orlandini, M. 2006, Adv. Space Res., 38, 2742 [NASA ADS] [CrossRef] [Google Scholar]
 Parmar, A. N., Martin, D. D. E., Bavdaz, M., et al. 1997, A&AS, 122, 309 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pavlov, G. G., & Panov, A. N. 1976, Sov. J. Exp. Theor. Phys., 44, 300 [NASA ADS] [Google Scholar]
 Pomraning, G. C. 1973, The equations of radiation hydrodynamics (Oxford: Pergamon Press) [Google Scholar]
 Pringle, J. E., & Rees, M. J. 1972, A&A, 21, 1 [NASA ADS] [Google Scholar]
 Psaltis, D., & Lamb, F. K. 1997, ApJ, 488, 881 (PL97) [NASA ADS] [CrossRef] [Google Scholar]
 Raichur, H., & Paul, B. 2008, MNRAS, 387, 439 [NASA ADS] [CrossRef] [Google Scholar]
 Rappaport, S., Clark, G. W., Cominsky, L., Li, F., & Joss, P. C. 1978, ApJ, 224, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Reynolds, A. P., Quaintrell, H., Still, M. D., et al. 1997, MNRAS, 288, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Riffert, H., Klingler, M., & Ruder, H. 1999, Phys. Rev. Lett., 82, 34321 (R99) [CrossRef] [Google Scholar]
 Santangelo, A., Segreto, A., Giarrusso, S., et al. 1999, ApJ, 523, L85 [NASA ADS] [CrossRef] [Google Scholar]
 Scott, D., Leahy, D., & Wilson, R. 2000, ApJ [Google Scholar]
 Shakura, N. I., Postnov, K. A., & Prokhorov, M. E. 1991, Sov. Astron. Lett., 17, 339 [NASA ADS] [Google Scholar]
 Staubert, R., Shakura, N. I., Postnov, K., et al. 2007, A&A, 465, L25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Staubert, R., Klochkov, D., Postnov, K., et al. 2009, A&A, 494, 1025 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Staubert, R., Klochkov, D., Vasco, D., et al. 2013, A&A, 550, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Staubert, R., Klochkov, D., Wilms, J., et al. 2014, A&A, 572, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suchy, S., Pottschmidt, K., Wilms, J., et al. 2008, ApJ, 675, 1487 [NASA ADS] [CrossRef] [Google Scholar]
 Sunyaev, R. A., & Titarchuk, L. G. 1980, A&A, 86, 121 [NASA ADS] [Google Scholar]
 Takahashi, T., Mitsuda, K., Kelley, R., et al. 2010, in SPIE Conf. Ser., 7732 [Google Scholar]
 Tananbaum, H., Gursky, H., Kellogg, E. M., et al. 1972, ApJ, 174, L143 [NASA ADS] [CrossRef] [Google Scholar]
 Thompson, T. W. J., & Rothschild, R. E. 2009, ApJ, 691, 1744 [NASA ADS] [CrossRef] [Google Scholar]
 Titarchuk, L. G. 1988, Sov. Astron. Lett., 14, 229 [NASA ADS] [Google Scholar]
 Titarchuk, L., Mastichiadis, A., & Kylafis, N. D. 1997, ApJ, 487, 834 [NASA ADS] [CrossRef] [Google Scholar]
 Truemper, J., Sacco, B., Pietsch, W., et al. 1977, Astronomische Gesellschaft and Deutsche Forschungsgemeinschaft, 42, 120 [NASA ADS] [Google Scholar]
 Tsygankov, S. S., Lutovinov, A. A., Churazov, E. M., & Sunyaev, R. A. 2006, MNRAS, 371, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Tsygankov, S. S., Lutovinov, A. A., Churazov, E. M., & Sunyaev, R. A. 2007, Astron. Lett., 33, 368 [NASA ADS] [CrossRef] [Google Scholar]
 van der Meer, A., Kaper, L., van Kerkwijk, M. H., Heemskerk, M. H. M., & van den Heuvel, E. P. J. 2007, A&A, 473, 523 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vasco, D., Staubert, R., Klochkov, D., et al. 2013, A&A, 550, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ventura, J. 1979, Phys. Rev. D, 19, 1684 [NASA ADS] [CrossRef] [Google Scholar]
 Walter, R., Lutovinov, A. A., Bozzo, E., & Tsygankov, S. S. 2015, A&ARv, 23, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Wasserman, I., & Shapiro, S. L. 1983, ApJ, 265, 1036 [NASA ADS] [CrossRef] [Google Scholar]
 Wheaton, W. A., Doty, J. P., Primini, F. A., et al. 1979, Nature, 282, 240 [NASA ADS] [CrossRef] [Google Scholar]
 White, N. E., Swank, J. H., & Holt, S. S. 1983, ApJ, 270, 711 [NASA ADS] [CrossRef] [Google Scholar]
 Whitlock, L., RousselDupre, D., & Priedhorsky, W. 1989, ApJ, 338, 381 [NASA ADS] [CrossRef] [Google Scholar]
 Wojdowski, P. S., Liedahl, D. A., Sako, M., Kahn, S. M., & Paerels, F. 2003, ApJ, 582, 959 [NASA ADS] [CrossRef] [Google Scholar]
 Zane, S., Ramsay, G., JimenezGarate, M. A., Willem den Herder, J., & Hailey, C. J. 2004, MNRAS, 350, 506 [NASA ADS] [CrossRef] [Google Scholar]
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_{e}c^{2}, 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 cosinedirection 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 X3. Also show are the bestfit compmag model and the residuals from the fit in units of σ (bottom panel). All bestfit parameters determined from the fit are given in Table 2. 

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

Open with DEXTER  
In the text 
Fig. 3 Same as Fig. 1 but for the combined Suzaku and NuStar observation of Her X1. 

Open with DEXTER  
In the text 
Fig. 4 From top to bottom: comptonized bremsstrahlung and cyclotron spectra of Cen X3, 4U 0115+63, and Her X1 obtained from the bestfit models given in Table 2. 

Open with DEXTER  
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 Xray 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). 

Open with DEXTER  
In the text 