A GENERALIZED ANALYTICAL MODEL FOR THERMAL AND BULK COMPTONIZATION IN ACCRETION-POWERED X-RAY PULSARS

1 INTRODUCTION

The X-ray emission from binary pulsars such as Her X-1, LMC X-4, Cen X-3, GX 304-1, and X Per is powered by the accretion of material from the “normal” companion onto the neutron star. The inflowing material is channeled by the strong magnetic field onto one or both of the neutron star’s magnetic poles. The X-ray luminosity, $L_{X}$ , of accretion-powered X-ray pulsars extends over a very wide range, from $L_{X}\sim 10^{34}\,{\rm ergs\,s}^{-1}$ for X Per up to $L_{X}\sim 10^{38}\,{\rm ergs\,s}^{-1}$ for LMC X-1. Modeling the accretion of material from the binary stellar companion onto the surface of the neutron star, and the formation of the associated X-ray spectrum, is a challenging physical problem that has not yet been fully solved. Physical simulation of these systems requires consideration of high-energy plasma processes in strong magnetic fields, strong shock physics involving both gas-dominated and radiation-dominated shocks, general relativistic light bending in the gravitational potential well of the neutron star, and complicated multi-dimensional radiative transfer. Many attempts have been made to relate the formation of the spectrum in X-ray pulsars to the physics of the accretion dynamics, based on either gas-mediated shocks (Langer & Rappaport, 1982), or Coulomb collisional stopping (e.g., Miller et al., 1989), but these have not demonstrated good agreement with X-ray pulsar spectral data (e.g., Meszaros & Nagel, 1985a, b).

More recently, the development of a new class of models that focuses on the effects of Comptonization in the accretion columns has significantly improved the agreement between the theoretical predictions and the observed phase-averaged X-ray spectra for a number of sources. The most comprehensive physical model currently available for the formation of the spectra in accretion-powered X-ray pulsars was developed by Becker & Wolff (2007, hereafter BW07), who studied the reprocessing of cyclotron, bremsstrahlung, and blackbody seed photons injected into a cylindrical accretion column located over one of the neutron star’s magnetic poles. The model is based on a detailed treatment of Comptonization, which is the energization of photons via electron scattering. This process can be broken into two components, namely thermal Comptonization (due to the stochastic velocity of the electrons), and bulk Comptonization (due to the accretion velocity of the electrons). Photons injected with a monochromatic energy distribution subsequently propagate throughout the energy space as a result of Comptonization, and they also diffuse through the physical space as they execute a random walk by scattering off the electrons, until they eventually escape through the walls or the top of the accretion column. Although the process is stochastic, in an average scattering event, net energy is transferred from the electrons to the photons, and the escape of the upscattered seed photons through the walls or top of the column carries away the kinetic energy of the accreting material in the form of X-rays, thereby allowing the gas to settle onto the star. This mechanism characteristically produces a power-law continuum in the energy range $\sim 1-30\,{\rm keV}$ , with a quasi-exponential cutoff (due to electron recoil) at higher energies. BW07 successfully applied their model to the interpretation of the phase-averaged X-ray spectra observed from Her X-1, Cen X-3, and LMC X-4, hence providing for the first time a firm theoretical foundation for the formation of the X-ray spectra emitted by luminous X-ray pulsars. Related work, reaching similar conclusions, was carried out by Farinelli et al. (2016), Wolff et al. (2016), West et al. (2017a, b), and Thalhammer et al. (2021).

The BW07 model has demonstrated success in applications to high-luminosity accretion-powered X-ray pulsars, with relatively flat power-law spectra, but it has proven more difficult to apply the model to low-luminosity sources such as X Per, which display steeper spectra. The problems with the low-luminosity sources stem from the fact that in the BW07 model, the accretion velocity is assumed to approach zero at the neutron star surface as a consequence of the strong deceleration occurring in an extended, smooth, radiation-dominated shock wave. This type of velocity profile is expected in high-luminosity “supercritical” sources such as Her X-1 (Becker et al., 2012), in which the radiation luminosity exceeds the effective Eddington limit for the gas, which is usually referred to as the critical luminosity, given by (e.g., Burnard et al., 1991)

L_{\rm crit}=L_{\rm E}\,\frac{\pi r_{0}^{2}}{4\pi R_{*}^{2}}\,\frac{\sigma_{{}_{\rm T}}}{\sigma_{||}}\ ,

(1)

where $r_{0}$ denotes the radius of the cylindrical accretion column, $R_{*}$ is the stellar radius, $\sigma_{{}_{\rm T}}$ is the Thomson cross section, $\sigma_{{}_{||}}$ denotes the electron scattering cross section for photons propagating parallel to the axis of the accretion column (which is also the magnetic field axis), and $L_{\rm E}$ represents the standard spherical Eddington limit for pure, fully-ionized hydrogen, defined by

L_{\rm E}\equiv\frac{4\pi GM_{*}m_{p}c}{\sigma_{{}_{\rm T}}}=1.26\times 10^{38}\left(\frac{M_{*}}{M_{\odot}}\right)\,{\rm ergs\ sec}^{-1}\ ,

(2)

with $m_{p}$ denoting the proton mass, $c$ the speed of light, and $M_{*}$ the stellar mass. The second factor on the right-hand side of Equation (1) represents the reduction in the critical luminosity due to the area of the accretion column, compared with the stellar surface area, and the third factor represents the increase in the critical luminosity, resulting from the fact that $\sigma_{||}\lesssim\sigma_{{}_{\rm T}}$ in a strong magnetic field, for energies significantly below the cyclotron energy (e.g., Ventura, 1979; Meszaros & Ventura, 1979).

The accretion dynamics for sources with X-ray luminosity $L_{X}\gtrsim L_{\rm crit}$ , such as Her X-1, are expected to be dominated by radiation pressure, with the gas decelerating essentially to rest at the stellar surface, after passing through a radiative, radiation-dominated standing shock (Becker & Wolff, 2005a, b). Conversely, for sources with $L_{X}\ll L_{\rm crit}$ , such as X Per, the role of radiation pressure is greatly reduced, and the gas may collide with the stellar surface with a substantial fraction of the local free-fall velocity. For sources with luminosity $L_{X}\sim L_{\rm crit}$ , the situation is more complex, and the accretion dynamics will be affected by additional details, such as the dependence of the electron scattering cross section on the photon energy and propagation direction (Becker et al., 2012; Mushtukov et al., 2015).

An interesting trend emerging from observations of accretion-powered X-ray pulsars over the past decades suggests that the index of the power-law continuum tends to reduce (i.e., the spectrum becomes flatter and harder) with increasing luminosity. In the case of Her X-1, this same trend is observed both in the long-term variability of the source luminosity, and also in the pulse-to-pulse variability (Klochkov et al., 2011). In the context of the BW07 model, this type of spectral variability with changes in the luminosity could be a natural consequence of changes in the magnitude of the radiation pressure, which ultimately controls the impact velocity of the material accreting onto the surface of the neutron star, and therefore determines the amount of $P$ d $V$ work done on the radiation field by the compressing gas. This effect is stronger in the high-luminosity (supercritical) sources, because of the enhanced compression, resulting in a flatter continuum. On the other hand, in the subcritical sources, the compression is weaker and the resulting X-ray spectrum is steeper and softer. We discuss this idea further in Section 10.

Motivated by the successes and limitations of the BW07 model, our goal in this paper is to develop a new, generalized model that expands upon the BW07 model to include a variety of new enhancements. Namely, the new model: (1) utilizes a conical geometry rather than the cylindrical geometry employed by BW07; (2) allows the computation of separate radiation components emitted through the walls and top of the column; (3) includes a new boundary condition at the stellar surface that allows for any impact velocity between free-fall and zero velocity; (4) incorporates a new velocity profile that smoothly merges with Newtonian free-fall far from the star; (5) utilizes a proper free-streaming boundary condition at the top of the accretion column; and (6) allows for the possibility of a standing shock at the column top.

We solve the steady-state radiation transport equation in a conical geometry, including the effects of thermal and bulk Comptonization, to obtain the analytical solution for the Green’s function, which represents the contribution to the observed steady-state photon spectrum resulting from the continual injection of monochromatic seed photons. By exploiting the linear structure of the mathematical problem, we show how the Green’s function can be convolved with an arbitrary source distribution to obtain the particular solution for the steady-state spectrum resulting from the continual injection of photons from any physical source mechanism. Examples include bremsstrahlung, blackbody, and cyclotron emission. The formalism also allows us to the calculate the separate radiation components emitted through the walls and top of the accretion column, which facilitates the computation of phase-dependent X-ray spectra, although we do not perform such calculations here.

We demonstrate that our new theoretical model is able to successfully reproduce the observed spectra for X-ray pulsars across five orders of magnitude in luminosity, from low-luminosity (subcritical) sources, up to relatively high-luminosity (supercritical) sources. We find that the spectra of the high-luminosity sources is dominated by Comptonized bremsstrahlung emission, and the spectra of the low-luminosity sources is dominated by Comptonized blackbody emission, powered by the residual kinetic energy of the flow at the stellar surface, as first suggested by Becker & Wolff (2005b). We illustrate the range of application of the new model by using it to qualitatively fit the X-ray continuum spectra for the supercritical source Her X-1 and the subcritical source X Per.

The remainder of the paper is organized as follows. In Section 2 we briefly review the nature of the primary radiation transport mechanisms in the accretion column with a focus on dynamical, thermal, and magnetic effects. The transport equation governing the formation of the radiation spectrum is introduced and analyzed in Section 3, and in Section 4 the exact analytical solution for the Green’s function describing the radiation distribution inside the accretion column is derived. The spectrum of the radiation escaping through the walls and top of the accretion column is developed in Section 5, and the physical constraints for the various model parameters are considered in Section 6. The nature of the source terms describing the injection of blackbody, cyclotron, and bremsstrahlung seed photons into the accretion column is discussed in Section 7. Emergent X-ray spectra are computed in Section 8, and the results are compared with the observational data for Her X-1 and X Per, which have widely differing luminosities. The self-consistency of the model is investigated in Section 9, and the implications of our work for the production of X-ray spectra in accretion-powered X-ray pulsars are discussed in Section 10.

2 RADIATIVE PROCESSES

The formation of the emergent X-ray continuum spectra in accretion-powered X-ray pulsars is a complex process that is powered fundamentally by the conversion of gravitational potential energy into kinetic energy of the accreting gas, which is transferred to the radiation field via electron scattering, and ultimately carried away by the energy of the escaping radiation. Accretion onto the surface of a neutron star differs qualitatively from black-hole accretion in the sense that the star obviously possesses a hard surface, and the accreting gas must therefore eventually decelerate to merge with the stellar crust. Conversely, in the case of black-hole accretion, no such solid solid barrier exists, although a centrifugal barrier may still develop in the flow due to a balance between centripetal and gravitational forces (e.g., Le & Becker, 2004). When an obstacle exists in the flow, whether due to a solid surface or due to a centrifugal barrier, shocks may form. In the case of neutron star accretion, the nature of the shock depends primarily on the accretion rate, which determines the luminosity of the escaping radiation field. In low-luminosity sources, the shock is expected to take the form of a discontinuous gas-mediated shock (Langer & Rappaport, 1982), and in high-luminosity sources, the shock is likely to be radiation-dominated, smooth, and radiative, meaning that the kinetic energy of the inflowing gas is radiated away in the shock itself.

The accretion flow in an X-ray pulsar is illustrated schematically in Figure 1. Physically, the accretion scenario corresponds to the flow of a mixture of gas and radiation inside a dipole-shaped magnetic “pipe” in which the fully-ionized gas is trapped by the strong magnetic field, but from which the radiation can escape via a three-dimensional random walk mediated by electron scattering. The primary mechanisms for the production of seed photons in X-ray pulsar accretion columns are cyclotron, bremsstrahlung, and blackbody emission. The seed photons are subsequently Compton scattered in energy due to collisions with hot electrons in the accreting gas, before escaping through the walls or top of the column. Bremsstrahlung emission provides a broad-band source of seed photons that are generated throughout the accretion column. Cyclotron emission also generates seed photons throughout the column, but rather than being a broad-band source, this process generates photons with an energy equal to the difference between the energy of the first excited Landau level and the ground state. Finally, blackbody seed photons are emitted from the surface of the “thermal mound” located close to the stellar surface, where the accreting gas gets sufficiently dense that thermodynamic equilibrium is achieved. Above the thermal mound, the opacity is dominated by electron scattering.

Refer to caption — Figure 1: Schematic depiction of gas in a conical accretion column accreting onto the magnetic pole of a neutron star. Seed photons created via blackbody (gold lines), cyclotron (red lines), and bremsstrahlung (cyan lines) emission are reprocessed via electron scattering and eventually escape through the walls and top of the column to form the observed X-ray spectrum (green lines).

2.1 Radiative Transfer in Magnetized Media

The strong ( $B\sim 10^{12}\,$ G) magnetic field channels the flow of gas from the accretion disk surrounding the neutron star onto one (or both) of the star’s magnetic poles. The presence of the strong magnetic field in X-ray pulsars has profound implications not only for the dynamics of the accretion flow, but also for the nature of the photon propagation through the accreting plasma. For example, vacuum polarization leads to birefringent behavior that produces two linearly polarized normal modes (Ventura, 1979; Nagel, 1980; Chanan et al., 1979). The electric field vector for the radiation is located in the plane formed by the photon propagation direction and the neutron star’s magnetic field for the ordinary mode. Conversely, for the extraordinary mode, the electric field vector of the radiation is pointed perpendicular to this plane. The nature of the photon-electron scattering process differs qualitatively for the two polarization modes, and it also depends critically on whether the photon energy, $\epsilon$ , exceeds the cyclotron energy, $\epsilon_{c}$ , defined by

\epsilon_{c}\equiv\frac{eBh}{2\pi m_{e}c}=11.57\left(\frac{B}{10^{12}\,{\rm G}}\right)\ {\rm keV}\ ,

(3)

where $c$ , $h$ , $m_{e}$ , and $e$ represent the speed of light, Planck’s constant, and the electron mass and charge, respectively.

2.1.1 Electron Scattering Cross Sections in Magnetized Plasma

Ventura (1979) derived expressions for the electron scattering cross sections for extraordinary and ordinary mode photons propagating in a magnetized plasma, neglecting the effects of vacuum polarization. The results he obtained for the plasma-only cross sections are valid provided the photon energy, $\epsilon$ , greatly exceeds the plasma energy, $\epsilon_{p}$ , given by

\epsilon_{p}=\frac{h\omega_{p}}{2\pi}=0.371\,\left(\frac{n_{e}}{10^{26}\,{\rm cm}^{-3}}\right)^{1/2}\,{\rm keV}\ ,

(4)

where $n_{e}$ denotes the electron number density, and the electron plasma frequency, $\omega_{p}$ , is defined by

\omega_{p}\equiv\left(\frac{4\pi n_{e}e^{2}}{m_{e}}\right)^{1/2}\ .

(5)

Our primary focus here is on the pulsars X Per and Her X-1, in which case the electron number density, $n_{e}$ , at the base of the accretion column is in the range $10^{18}\,{\rm cm}^{-3}\lesssim n_{e}\lesssim 10^{24}\,{\rm cm}^{-3}$ , with the lower limit corresponding to X Per and the upper limit to Her X-1. Substituting this range of number densities into Equation (4) leads to the conclusion that $\epsilon_{p}\ll 1\,$ keV for both of the sources. Hence $\epsilon\gg\epsilon_{p}$ in the X-ray pulsar application, and therefore the expressions derived by Ventura (1979) are applicable, provided the effects of vacuum polarization are negligible, which we discuss further in Section 2.1.2. The primary results derived by Ventura (1979) are summarized below.

In the absence of vacuum polarization effects, and assuming that $\epsilon\gg\epsilon_{p}$ , Ventura (1979) demonstrated that the pure-plasma electron scattering cross sections for extraordinary and ordinary mode photons, denoted by $\sigma_{{}_{P1}}$ and $\sigma_{{}_{P2}}$ , respectively, can be written as

\frac{\sigma_{{}_{P1}}}{\sigma_{{}_{\rm T}}}=\left[1+\alpha^{2}(\theta)\right]^{-1}\left\{\alpha^{2}(\theta)\sin^{2}\theta+\frac{1}{2}\left[\frac{1+\alpha(\theta)\cos\theta}{1+u^{1/2}}\right]^{2}+\frac{1}{2}\left[\frac{1-\alpha(\theta)\cos\theta}{1-u^{1/2}}\right]^{2}\right\}\ ,

(6)

and

\frac{\sigma_{{}_{P2}}}{\sigma_{{}_{\rm T}}}=\left[1+\alpha^{2}(\theta)\right]^{-1}\left\{\sin^{2}\theta+\frac{1}{2}\left[\frac{\cos\theta-\alpha(\theta)}{1+u^{1/2}}\right]^{2}+\frac{1}{2}\left[\frac{\cos\theta+\alpha(\theta)}{1-u^{1/2}}\right]^{2}\right\}\ ,

(7)

where $\theta$ is the angle between the magnetic field vector and the photon propagation direction, and the function $\alpha(\theta)$ is defined by

\alpha(\theta)\equiv\frac{-b}{1+(1+b^{2})^{1/2}}\ ,

(8)

with

b=2\,u^{-1/2}\,\frac{\cos\theta(1-v)}{\sin^{2}\theta}\ ,

(9)

and

u\equiv\left(\frac{\epsilon}{\epsilon_{c}}\right)^{-2}\ ,\qquad v\equiv\left(\frac{\epsilon}{\epsilon_{p}}\right)^{-2}\ .

(10)

We note that $v$ is essentially zero in the X-ray pulsar application since $\epsilon\gg\epsilon_{p}$ . In order to validate our implementation of the cross sections, in Figure 2 we plot $\sigma_{{}_{P1}}$ and $\sigma_{{}_{P2}}$ , evaluated using Equation (6) and (7), respectively. The results agree with Figure 2 from Ventura (1979). The extraordinary mode cross section exhibits a strong resonance at the cyclotron energy, as expected, while the ordinary mode cross section is non-resonant. Near the resonance, the extraordinary mode cross section greatly exceeds the Thomson value, whereas the ordinary mode cross section remains essentially sub-Thomson at all photon energies. We shall see below that the results for the scattering cross section are qualitatively different once the effects of vacuum polarization are considered.

2.1.2 Effects of Vacuum Polarization

The results obtained by Ventura (1979) and presented in Equations (6) and (7) are applicable when the effects of vacuum polarization are negligible, which corresponds to the energy range $\epsilon\lesssim\epsilon_{\rm vac}$ , where the vacuum polarization energy, $\epsilon_{\rm vac}$ , is given by (Meszaros & Ventura, 1979)

\epsilon_{\rm vac}\equiv\frac{h\omega_{\rm vac}}{2\pi}=1.32\left(\frac{n_{e}}{10^{20}\,{\rm cm}^{-3}}\right)^{1/2}\left(\frac{B}{10^{12}\,{\rm G}}\right)^{-1}\,{\rm keV}\ ,

(11)

and the vacuum polarization frequency, $\omega_{\rm vac}$ , is defined by

\omega_{\rm vac}\equiv\left(\frac{15\pi}{\alpha_{F}}\right)^{1/2}\left(\frac{B}{B_{c}}\right)^{-1}\omega_{p}\ .

(12)

Here, $\alpha_{F}$ and $\omega_{p}$ denote the fine-structure constant and the electron plasma frequency (Equation (5)), respectively, and $B_{c}$ represents the critical magnetic field strength, given by

B_{c}\equiv\frac{2\pi m_{e}^{2}c^{3}}{eh}=4.415\times 10^{13}\,{\rm G}\ ,

(13)

which is obtained by setting $\epsilon_{c}=m_{e}c^{2}$ (see Equation (3)). Vacuum polarization effects will strongly influence the electron scattering cross section when the photon energy $\epsilon\gtrsim\epsilon_{\rm vac}$ . In the X-ray pulsar accretion columns of interest here, $B\sim 10^{12}\,$ G, and the electron number density, $n_{e}$ , lies in the range $10^{18}\,{\rm cm}^{-3}\lesssim n_{e}\lesssim 10^{24}\,{\rm cm}^{-3}$ , where the lower and upper limits correspond to X Per and Her X-1, respectively. According to Equation (11), the corresponding range for the vacuum polarization energy is therefore $0.1\,{\rm keV}\lesssim\epsilon_{\rm vac}\lesssim 100\,{\rm keV}$ . Focusing on the case of Her X-1, with $\epsilon_{\rm vac}\sim 100\,{\rm keV}$ , it is clear that vacuum polarization effects will be unimportant for the X-ray continuum in this source. On the other hand, in the case of X Per, with $\epsilon_{\rm vac}\sim 0.1\,{\rm keV}$ , we note that vacuum polarization effects will be very important for the formation of the entire X-ray continuum. We summarize the results for the electron scattering cross sections including the effects of vacuum polarization below.

Meszaros & Ventura (1979) derived expressions for the electron scattering cross sections including the modifications due to vacuum polarization. The results they obtained for the cross sections for extraordinary and ordinary mode photons, denoted by $\sigma_{{}_{V1}}$ and $\sigma_{{}_{V2}}$ , respectively, can be written as

	$\displaystyle\frac{\sigma_{{}_{V1}}}{\sigma_{{}_{\rm T}}}=(1+K_{1}^{2})^{-1}$	$\displaystyle\bigg{\{}\frac{1}{2}(1-K_{1}\cos\theta)^{2}[(1-u^{1/2})^{2}+\gamma^{2}]^{-1}$		(14)
		$\displaystyle+\frac{1}{2}(1+K_{1}\cos\theta)^{2}[1+u^{1/2}]^{-2}+K_{1}^{2}\sin^{2}\theta\bigg{\}}\ ,$		(14)

and

	$\displaystyle\frac{\sigma_{{}_{V2}}}{\sigma_{{}_{\rm T}}}=(1+K_{2}^{2})^{-1}$	$\displaystyle\bigg{\{}\frac{1}{2}(1-K_{2}\cos\theta)^{2}[(1-u^{1/2})^{2}+\gamma^{2}]^{-1}$		(15)
		$\displaystyle+\frac{1}{2}(1+K_{2}\cos\theta)^{2}[1+u^{1/2}]^{-2}+K_{2}^{2}\sin^{2}\theta\bigg{\}}\ ,$		(15)

where the radiation damping constant $\gamma$ for photon frequency $\omega$ is defined by

\gamma\equiv\frac{2}{3}\frac{e^{2}\omega}{m_{e}c^{3}}\ ,

(16)

and the quantities $K_{1}$ and $K_{2}$ are computed using

K_{1}=K_{2}^{-1}=\alpha_{V}(\theta)\ ,

(17)

with

\alpha_{V}(\theta)\equiv\frac{-1}{{b_{V}+(1+b_{V}^{2})^{1/2}}}\ ,

(18)

and

b_{V}=\frac{u^{1/2}\sin^{2}\theta}{2\cos\theta(1-v)}\left[1+\frac{3\delta(1-u)}{vu}\right]\ .

(19)

The quantities $u$ and $v$ are defined in Equation (10), and the parameter $\delta$ is computed using

\delta\equiv\frac{2e^{2}}{45hc}\left(\frac{B}{B_{c}}\right)^{2}=5.16\times 10^{-5}\left(\frac{B}{B_{c}}\right)^{2}\ ,

(20)

where the critical field, $B_{c}$ , is defined in Equation (13).

In order to validate our implementation of the vacuum-modified cross sections, in Figure 3 we utilize Equations (14) and (15) to plot $\sigma_{{}_{V1}}$ and $\sigma_{{}_{V2}}$ using the parameter values $n_{e}=1.89\times 10^{21}\,{\rm cm}^{-3}$ and $B=4.41\times 10^{12}\,$ G, which are the values used by Meszaros & Ventura (1979) in their Figure 4. For these parameters, we obtain $\epsilon_{c}=51\,$ keV and $\epsilon_{\rm vac}=1.3\,$ keV, according to Equations (3) and (11), respectively. Hence in this case vacuum polarization will strongly modify the scattering cross sections for photons with energy $\epsilon\gtrsim 1.3\,$ keV, which comprises the entire X-ray band. This behavior is clearly seen in Figure 3, especially for energies close to the cyclotron resonance, where we note that the vacuum-modified extraordinary mode cross section (solid lines) displays a rapid increase as a function of the propagation angle $\theta$ , from sub-Thomson values for $\theta\sim 0$ to highly super-Thomson values for $\theta\sim 90^{\circ}$ . This behavior differs qualitatively from that of the plasma-only extraordinary-mode cross section (dot-dashed lines), which does not display a strong dependence on $\theta$ for any value of the photon energy.

In Figure 4 we use Equations (14) and (15) to plot the vacuum-modified electron scattering cross sections for the extraordinary and ordinary mode photons as functions of the photon energy ratio, $\epsilon/\epsilon_{c}$ , for the same values of $n_{e}$ and $B$ used in Figure 3. Due to the effect of vacuum polarization, the extraordinary mode cross section is now a very weak function of the propagation angle $\theta$ , for $\theta\gtrsim 20^{\circ}$ , and the ordinary mode now displays a resonance at the cyclotron energy, which was absent in the plasma-only cross section plotted in Figure 2. It is interesting to note that the value of the electron number density used in Figures 3 and 4 lies between the values expected at the base of the accretion columns in Her X-1 and X Per, and therefore the behaviors displayed in these two figures are likely to be intermediate between the two sources. We will consider the cross section values in more detail in Section 9.

2.1.3 Simplified Scattering Cross Sections

The detailed energy and angular dependences of the electron scattering cross sections for the extraordinary and ordinary polarization modes, including the effects of vacuum polarization, are given by Equations (14) and (15), respectively. Unfortunately, it is not possible to include the full complexity of these expressions into the model developed here, and we will therefore follow BW07 and Wang & Frank (1981) by introducing a set of approximate scattering cross sections, averaged over the photon energy and the two polarization modes. The cross sections for photons propagating parallel and perpendicular to the radial direction are denoted by $\sigma_{{}_{||}}$ and $\sigma_{\perp}$ , respectively, and the angle-averaged cross section is denoted by $\overline{\sigma}$ . The magnitudes of these cross sections are fairly well understood for relatively luminous sources such as Her X-1, in which case one usually finds that $\sigma_{\perp}\sim\sigma_{{}_{\rm T}}$ and $\sigma_{\perp}\gg\overline{\sigma}\gg\sigma_{{}_{||}}$ (e.g., Wolff et al., 2016). The hierarchy of these cross section values stems from the fact that $\epsilon<\epsilon_{c}<\epsilon_{\rm vac}$ in Her X-1, and therefore the super-Thomson cross section values near the cyclotron resonance are avoided (see Figures 2 and 3). Conversely, in the case of low-luminosity sources such as X Per, the situation is reversed, and we generally find that $\epsilon_{\rm vac}<\epsilon<\epsilon_{c}$ . The change in the energy hierarchy leads to a greatly increased contribution from the cyclotron resonance, which results in cross section values $\sigma_{\perp}\gg\sigma_{{}_{\rm T}}$ and $\sigma_{{}_{||}}\lesssim\overline{\sigma}\ll\sigma_{\perp}$ . This issue is further discussed in Section 9.

3 RADIATIVE TRANSFER

The exact solution for the radiation spectrum inside a cylindrical X-ray pulsar accretion column was obtained by BW07 under the assumption that the flow velocity of the accreting gas approaches zero at the stellar surface, as a result of passage through a radiative, radiation-dominated shock wave. This assumption leads to a theoretical model that produces X-ray spectra that agree quite well with the relatively flat power-law spectra observed from luminous X-ray pulsars such as Her X-1 and LMC X-4. However, attempts to fit the same model to low-luminosity sources such as X Per, with steeper spectra, yield parameter values that are difficult to justify physically. We suspect that the problems with the low-luminosity sources stem from a need to generalize the boundary conditions utilized in the model. This observation has motivated a comprehensive reexamination of the boundary conditions, which in turn led to a complete reformulation of the model that not only generalizes the boundary conditions, but also incorporates a more realistic conical geometry for the accretion column. In addition, the new model also includes a more realistic velocity profile, and the capability to separately compute the radiation components emitted through the walls and the top of the accretion column. The transport equation introduced below includes the effects of special relativity up to first order in $\upsilon/c$ , where $\upsilon$ is the accretion velocity. However, the effects of general relativity are not included, as this would greatly complicate the model, and the resulting corrections are not likely to be important at the level of approximation employed here. We present the core elements of the new model in this section.

3.1 Transport Equation

The time-dependent transport equation governing the propagation, scattering, and escape of radiation in a pulsar accretion column with an arbitrary geometry can be written in the vector form

	$\displaystyle\frac{\partial f}{\partial t}+\vec{\upsilon}\cdot\vec{\nabla}f$	$\displaystyle=$	$\displaystyle(\vec{\nabla}\cdot\vec{\upsilon})\,\frac{\epsilon}{3}\,\frac{\partial f}{\partial\epsilon}+\vec{\nabla}\cdot\left(\frac{c}{3n_{e}\sigma_{{}_{\|\|}}}\,\vec{\nabla}\cdot f\right)$		(21)
		$\displaystyle+$	$\displaystyle\frac{n_{e}\overline{\sigma}c}{m_{e}c^{2}}\frac{1}{\epsilon^{2}}\frac{\partial}{\partial\epsilon}\left[\epsilon^{4}\left(f+kT_{e}\,\frac{\partial f}{\partial\epsilon}\right)\right]+Q\ ,$		(21)

where $\vec{\upsilon}$ is the accretion velocity field, $f(\vec{r},\epsilon,t)$ is the photon distribution function, $\vec{r}$ is the spatial vector, $t$ is time, $\epsilon$ is the photon energy, $n_{e}$ is the electron number density, $T_{e}$ is the electron temperature, $k$ is Boltzmann’s constant, $\sigma_{{}_{||}}$ represents the electron scattering cross section for photons propagating in the radial direction, and $\overline{\sigma}$ denotes the angle-averaged scattering cross section. The distribution function $f$ is related to the occupation number $\bar{n}$ via $f=8\pi\,\bar{n}/(c^{3}h^{3})$ , and the quantity $\epsilon^{2}f(\vec{r},\epsilon,t)\,d\epsilon$ gives the number density of photons in the energy range between $\epsilon$ and $\epsilon+d\epsilon$ . Hence the total photon number density, $n_{r}$ , and energy density, $U_{r}$ , are computed using the integrals

n_{r}(\vec{r},t)=\int_{0}^{\infty}\epsilon^{2}\,f(\vec{r},\epsilon,t)\,d\epsilon\ ,\qquad U_{r}(\vec{r},t)=\int_{0}^{\infty}\epsilon^{3}\,f(\vec{r},\epsilon,t)\,d\epsilon\ .

(22)

The left-hand side of Equation (21) denotes the comoving time derivative of $f$ , and the terms on the right-hand side represent first-order Fermi energization (“bulk Comptonization”), spatial diffusion, thermal Comptonization, and the radiation source term, $Q$ , respectively. We note that equation (21) is valid in the region of the accretion column above the thermal mound, where the opacity is dominated by electron scattering (see Section 6.2).

BW07 solved Equation (21) under the assumption of a cylindrical accretion column geometry. Here, we adopt a hollow conical geometry, in which the solid angle of the column, $\Omega$ , is given by

\Omega=2\pi(\cos\Theta_{2}-\cos\Theta_{1})\ ,

(23)

where $\Theta_{1}$ and $\Theta_{2}$ denote the opening angles for the outer and inner walls of the column, respectively. The conical shape is a more reasonable approximation of the true dipole geometry of the column, especially in the lower region close to the stellar surface. In the conical accretion column geometry, the accretion rate, $\dot{M}$ , is related to the electron number density, $n_{e}$ , the radius measured from the center of the neutron star, $R$ , and the radial accretion velocity, $\upsilon<0$ , via

\dot{M}=\Omega R^{2}n_{e}m_{p}|\upsilon|\ ,

(24)

where we have assumed that the accreting gas is composed of pure, fully-ionized hydrogen.

In order to treat a variety of photon injection scenarios, such as bremsstrahlung, cyclotron, and blackbody emission, it is convenient to compute the Green’s function, $f_{{}_{\rm G}}$ , which represents the photon distribution inside the accretion column, resulting from the continual injection of $\dot{N}_{0}$ monochromatic seed photons per unit time, each with energy $\epsilon_{0}$ , from a source located at radius $R_{0}$ . In a conical geometry, the equation satisfied by $f_{{}_{\rm G}}$ can be written as

	$\displaystyle\frac{\partial f_{{}_{\rm G}}}{\partial t}+\upsilon\,\frac{\partial f_{{}_{\rm G}}}{\partial R}$	$\displaystyle=$	$\displaystyle\frac{1}{R^{2}}\frac{\partial}{\partial R}\left(R^{2}\,\upsilon\right)\,\frac{\epsilon}{3}\,\frac{\partial f_{{}_{\rm G}}}{\partial\epsilon}+\frac{1}{R^{2}}\frac{\partial}{\partial R}\left(\frac{R^{2}\,c}{3n_{e}\sigma_{{}_{\|\|}}}\,\frac{\partial f_{{}_{\rm G}}}{\partial R}\right)-\frac{f_{{}_{\rm G}}}{t_{\rm esc}}$		(25)
		$\displaystyle+$	$\displaystyle\frac{n_{e}\overline{\sigma}c}{m_{e}c^{2}}\frac{1}{\epsilon^{2}}\frac{\partial}{\partial\epsilon}\left[\epsilon^{4}\left(f_{{}_{\rm G}}+kT_{e}\,\frac{\partial f_{{}_{\rm G}}}{\partial\epsilon}\right)\right]+\frac{\dot{N}_{0}\,\delta(\epsilon-\epsilon_{0})\,\delta(R-R_{0})}{\Omega R_{0}^{2}\,\epsilon_{0}^{2}}\ ,$		(25)

where $t_{\rm esc}$ represents the mean time for photons to escape through the walls of the accretion column. The escape of radiation through the top of the accretion column is implemented using a free-streaming boundary condition as discussed in Section 4.1. Equation (25) is obtained by adopting a conical geometry in Equation (21) and then averaging over the azimuthal direction in the column. Hence the only spatial variable appearing in Equation (25) is the central radius $R$ . We are interested in solving the steady-state version of Equation (25) and therefore we will set $\partial f_{{}_{\rm G}}/\partial t=0$ .

The escape of radiation through the walls of the conical accretion column is modeled using an escape-probability formalism, based on the mean escape timescale, $t_{\rm esc}$ , appearing in Equation (25). The value of $t_{\rm esc}$ is equal to the average time required for photons to diffuse through the walls of the column, which is computed using the expressions

t_{\rm esc}(R)=\frac{r_{\perp}}{\upsilon_{\rm diff}}\ ,\qquad\upsilon_{\rm diff}(R)=\frac{c}{\tau_{\perp}}\ ,

(26)

where $\upsilon_{\rm diff}$ represents the radiation diffusion velocity. The perpendicular scattering optical thickness, $\tau_{\perp}$ , and the perpendicular radius of the accretion column, $r_{\perp}$ , are given by

\tau_{\perp}(R)=n_{e}\,\sigma_{\perp}\,r_{\perp}\ ,\qquad r_{\perp}(R)=b\,R\ ,

(27)

where the constant $b$ is given by

b=\tan\Theta_{1}-\tan\Theta_{2}\ .

(28)

The diffusion approximation employed in Equation (26) is valid provided the perpendicular scattering optical thickness $\tau_{\perp}>1$ , and it is important to confirm that this is the case in our applications. By combining relations, Equation (26) for the mean escape timescale can be rewritten as

t_{\rm esc}(R)=\frac{b^{2}\sigma_{\perp}\dot{M}}{\Omega\,m_{p}\,c\,|\upsilon|}\ .

(29)

Following BW07 and Lyubarskii & Syunyaev (1982), we shall assume that the electron distribution is isothermal, with constant temperature $T_{e}$ , throughout the emission region. This is an acceptable approximation, since the electron temperature is expected to be regulated in a fairly small interval via thermal Comptonization (e.g., Sunyaev & Titarchuk, 1980; West et al., 2017a, b). Under the assumption of a constant electron temperature, it is convenient to work in terms of the dimensionless energy variable, $\chi$ , defined by

\chi(\epsilon)\equiv\frac{\epsilon}{kT_{e}}\ .

(30)

It is also useful to transform the spatial variable from the radius, $R$ , to the scattering optical depth, $\tau$ , measured in the radial direction, starting from the stellar surface. The differential $d\tau$ is related to $dR$ via

d\tau=n_{e}(R)\,\sigma_{{}_{||}}\,dR\ ,

(31)

which can be rewritten in terms of the dimensionless radius, $y$ , defined by

y(R)\equiv\frac{R}{R_{*}}\ ,

(32)

to obtain

d\tau=n_{e}(y)\,\sigma_{{}_{||}}R_{*}\,dy\ .

(33)

Integration of Equation (33) yields

\tau(y)=\int_{1}^{y}n_{e}(y^{\prime})\,\sigma_{{}_{||}}R_{*}\,dy^{\prime}\ ,

(34)

so that $\tau=0$ at the stellar surface, where $y=1$ .

Making the change of variable from $(R,\epsilon)$ to $(\tau,\chi)$ in Equation (25) and using Equation (29) to substitute for the escape timescale $t_{\rm esc}$ , we find after some algebra that the transport equation satisfied by the steady-state Green’s function can be written in the form

	$\displaystyle u\,\frac{\partial f_{{}_{\rm G}}}{\partial\tau}$	$\displaystyle=$	$\displaystyle\frac{1}{y^{2}}\,\frac{d}{d\tau}(y^{2}u)\,\frac{\chi}{3}\,\frac{\partial f_{{}_{\rm G}}}{\partial\chi}+\frac{1}{3y^{2}}\,\frac{\partial}{\partial\tau}\left(y^{2}\,\frac{\partial f_{{}_{\rm G}}}{\partial\tau}\right)-\xi^{2}u^{2}y^{2}f_{{}_{\rm G}}$		(35)
		$\displaystyle+$	$\displaystyle\frac{\overline{\sigma}}{\sigma_{{}_{\|\|}}}\,\frac{kT_{e}}{m_{e}c^{2}}\frac{1}{\chi^{2}}\frac{\partial}{\partial\chi}\left[\chi^{4}\left(f_{{}_{\rm G}}+\frac{\partial f_{{}_{\rm G}}}{\partial\chi}\right)\right]+\frac{\dot{N}_{0}\,\delta(\chi-\chi_{0})\,\delta(\tau-\tau_{0})}{\Omega\,c\,R_{0}^{2}\chi_{0}^{2}(kT_{e})^{3}}\ ,$		(35)

where $\tau_{0}$ and $\chi_{0}$ denote the injection optical depth and the dimensionless injection energy, respectively, and we have introduced the dimensionless accretion velocity, $u<0$ , defined by

u(\tau)\equiv\frac{\upsilon(\tau)}{c}\ .

(36)

The dimensionless escape constant, $\xi$ , introduced in Equation (35), parameterizes the rate of diffusion of photons through the walls of the conical accretion column, and is defined by

\xi\equiv\frac{\Omega\,R_{*}m_{p}c}{b\,\dot{M}(\sigma_{{}_{||}}\sigma_{\perp})^{1/2}}\ ,

(37)

which can be rewritten as

\xi=\frac{\Omega}{4\pi b}\left(\frac{\dot{M}}{\dot{M}_{\rm E}}\right)^{-1}\left(\frac{\sigma_{\perp}}{\sigma_{{}_{\rm T}}}\right)^{-1/2}\left(\frac{\sigma_{{}_{||}}}{\sigma_{{}_{\rm T}}}\right)^{-1/2}\left(\frac{R_{*}}{R_{g}}\right)\ ,

(38)

where $R_{g}=GM_{*}/c^{2}$ denotes the gravitational radius of the neutron star, and the Eddington accretion rate is defined by

\dot{M}_{\rm E}\equiv\frac{4\pi GM_{*}m_{p}}{\sigma_{{}_{\rm T}}c}=1.40\times 10^{17}\,\left(\frac{M_{*}}{M_{\odot}}\right)\ {\rm g\,s}^{-1}\ .

(39)

3.2 Separability

It is important to identify the physical situations that result in the separability of the transport equation (Equation (35)), since separability allows one to obtain the exact analytical solution for the steady-state Green’s function $f_{{}_{\rm G}}(\tau,\chi)$ . Equation (35) is separable if each of the terms containing energy operators has the same dependence on the scattering optical depth, $\tau$ . In this case, separability can be accomplished if the accretion velocity, $u$ , satisfies the equation

\frac{1}{y^{2}}\frac{d}{d\tau}(y^{2}u)=-\alpha\ ,

(40)

where $\alpha$ is a positive constant. By utilizing Equations (24), (33), and (34) to transform from $d\tau$ to $dy$ in Equation (40) and rearranging the resulting expression, we can show that

\frac{d}{dy}(y^{2}u)^{2}=\frac{2\alpha\dot{M}\sigma_{{}_{||}}}{\Omega R_{*}m_{p}c}\,y^{2}\ ,

(41)

which can be integrated to obtain

(y^{2}u)^{2}=\left(\frac{2R_{g}}{R_{*}}\right)[k_{\infty}^{2}(y^{3}-1)+k_{0}^{2}]\ ,

(42)

where $k_{0}$ is a constant of integration, and the parameter $\alpha$ is related to $k_{\infty}$ via

\alpha=\frac{3k_{\infty}^{2}\Omega\,R_{g}\,m_{p}c}{\sigma_{{}_{||}}\dot{M}}\ .

(43)

The expression for the dimensionless constant $\alpha$ can also be written as

\alpha=\frac{3k_{\infty}^{2}\Omega}{4\pi}\left(\frac{\dot{M}}{\dot{M}_{\rm E}}\right)^{-1}\left(\frac{\sigma_{{}_{||}}}{\sigma_{{}_{\rm T}}}\right)^{-1}\ ,

(44)

where the Eddington accretion rate $\dot{M}_{\rm E}$ is defined in Equation (39). Based on Equation (42), we find that the general expression for the accretion velocity profile, $u(y)$ , resulting from the separability condition can be written in the form

u(y)=\frac{\upsilon(y)}{c}=-\left(\frac{2R_{g}}{R_{*}}\right)^{1/2}\,\left[\frac{k_{\infty}^{2}(y^{3}-1)+k_{0}^{2}}{y^{4}}\right]^{1/2}\ .

(45)

It is important to discuss the physical interpretation of the dimensionless constants $k_{0}$ and $k_{\infty}$ introduced in Equation (42), which are connected with the asymptotic behaviors of the velocity profile function $u(y)$ as $y\to\infty$ or $y\to 1$ . We note that at the stellar surface ( $y=1$ ), Equation (45) yields

u(y)\to-k_{0}\left(\frac{2GM_{*}}{R_{*}c^{2}}\right)^{1/2}\equiv\frac{\upsilon_{*}}{c}\ ,\qquad y\to 1\ ,

(46)

where $\upsilon_{*}$ is the impact velocity. The constant $k_{0}$ is equal to the ratio of the impact velocity divided by the Newtonian free-fall velocity, and we therefore refer to $k_{0}$ as the “impact velocity constant.” Next we note that as $y\to\infty$ , the asymptotic behavior of Equation (45) is given by

u(y)\to-k_{\infty}\left(\frac{2GM_{*}}{R\,c^{2}}\right)^{1/2}\ ,\qquad y\to\infty\ ,

(47)

which establishes that the constant $k_{\infty}$ is equal to the ratio of the accretion velocity divided by the Newtonian free-fall velocity in the upstream region, far from the neutron star. Hence we refer to $k_{\infty}$ as the “velocity at infinity constant.” When treating high-luminosity (supercritical) pulsars, such as Her X-1, we set $k_{\infty}=1$ so that the velocity of the accreting gas merges smoothly with the Newtonian free-fall velocity profile as $y\to\infty$ . On the other hand, setting $k_{\infty}\sim 0.25$ , for example, results in a sub-free-fall velocity profile above the star. This may be appropriate when treating low-luminosity (subcritical) sources such as X Per, in which a collisionless shock is thought to be located at the top of the accretion column, so that the downstream (post-shock) velocity is reduced from the free-fall value by a factor of 4 (Langer & Rappaport, 1982).

The $P$ d $V$ work done on the radiation field via bulk Comptonization in the accretion column depends sensitively on the amount of compression occurring near the stellar surface, where the gas decelerates and impacts against the star. Lower impact velocities lead to higher levels of compression, and enhanced bulk Comptonization, resulting in flatter power-law continuum spectra (BW07). Since the parameter $k_{0}$ expresses the impact velocity relative to the free-fall velocity at the stellar surface, we can expect that the value of $k_{0}$ will be connected with the power-law slope of the X-ray continuum spectrum escaping through the walls and the top of the accretion column. High-luminosity sources such as Her X-1, with relatively flat continuum spectra, are expected to have values of $k_{0}$ close to zero, which maximizes the compression occurring at the base of the accretion column. Conversely, low-luminosity sources sources such as X Per, which display relatively steep power-law continuum spectra, are expected to have non-zero values of $k_{0}$ .

In Figure 5 we plot several examples of the input velocity function, $u(y)$ , computed using Equation (45). Note that a wide variety of accretion scenarios can be modeled, depending on the values selected for the impact velocity constant, $k_{0}$ , and the velocity at infinity constant, $k_{\infty}$ . For example, setting $k_{0}=0$ results in a gradual settling of the gas onto the neutron star surface, appropriate for high-luminosity sources such as Her X-1, in which radiation pressures dominates the accretion dynamics. On the other hand, setting $k_{0}\sim 0.4$ may provide a reasonable description of the flow velocity near the stellar surface in low-luminosity sources such as X Per, in which radiation pressure is probably insufficient to decelerate the gas, and instead the gas impacts on the stellar surface with a high residual velocity, equal to $0.64\,k_{0}\,c$ . In this situation, the final merger with the stellar crust occurs within the final few cm above the neutron star surface via Coulombic deceleration (e.g., Sokolova-Lapa et al., 2021). It is also interesting to examine the dependence of the input velocity profile (Equation (45)) on the parameter $k_{\infty}$ . The cases with $k_{\infty}=1$ plotted in Figure 5 merge smoothly with the Newtonian free-fall profile as $y\to\infty$ , and the cases with $k_{\infty}=0.25$ approximate the velocity profile in the post-shock region below a standing discontinuous shock located at the top of the accretion column.

3.3 Optical Depth Variation

It is useful to develop an expression for the variation of the electron scattering optical depth, $\tau$ , measured from the stellar surface, for photons propagating in the radial direction, as a function of the radius $R$ . We will work in terms of the dimensionless radius $y\equiv R/R_{*}$ introduced in Equation (32). Starting with Equation (33) for the differential optical depth, $d\tau$ , we can use Equation (24) to substitute for the electron number density, $n_{e}$ , to obtain

d\tau=\frac{\dot{M}\sigma_{{}_{||}}}{\Omega R_{*}y^{2}m_{p}|\upsilon(y)|}\,dy\ .

(48)

Next, we can utilize Equation (43) to eliminate $\dot{M}$ in Equation (48), which yields

d\tau=\frac{3k_{\infty}^{2}R_{g}}{\alpha R_{*}y^{2}|u(y)|}\,dy\ ,

(49)

where $u=\upsilon/c$ according to Equation (36), and $\alpha$ is given by Equation (44). Equation (49) can be integrated with respect to $y$ to obtain

\tau(y)=\int_{1}^{y}\frac{3k_{\infty}^{2}R_{g}}{\alpha R_{*}y^{\prime 2}|u(y^{\prime})|}\,dy^{\prime}\ .

(50)

We can use Equation (45) to substitute for the velocity profile, $u(y)$ , in Equation (50), which yields

\tau(y)=\frac{3}{\alpha\sqrt{2}}\sqrt{\frac{R_{g}}{R_{*}}}\int_{1}^{y}\frac{k_{\infty}^{2}}{\sqrt{k_{\infty}^{2}(y^{\prime 3}-1)+k_{0}^{2}}}\,dy^{\prime}\ .

(51)

The elliptic integral in Equation (51) can be evaluated analytically with the assistance of Equations (15.2.5) and (15.3.5) from Abramowitz & Stegun (1970), giving the closed-form result

\tau(y)=\frac{3\sqrt{2}\,k_{\infty}}{\alpha}\left(\frac{R_{g}}{R_{*}}\right)^{1/2}\left[G(k_{0},k_{\infty},1)-G(k_{0},k_{\infty},y)\right]\ ,

(52)

where the function $G$ is defined by

G(k_{0},k_{\infty},y)\equiv\,_{2}F_{1}\left(\frac{1}{6},\frac{1}{2};\frac{7}{6};\frac{1-k_{0}^{2}/k_{\infty}^{2}}{y^{3}}\right)\frac{1}{\sqrt{y}}\ ,

(53)

and ${}_{2}F_{1}$ denotes the hypergeometric function (Abramowitz & Stegun, 1970). Note that according to Equation (52), $\tau=0$ at the stellar surface ( $y=1$ ) as required. In Figure 6 we plot the variation of the optical depth $\tau(y)$ for a few representative values of the parameters $\alpha$ and $k_{0}$ . In each of the cases plotted in Figure 6, we have set $k_{\infty}=1$ , so that the flow velocity merges smoothly with the Newtonian free-fall profile as $y\to\infty$ .

4 SOLUTION FOR THE GREEN’S FUNCTION

Adopting the velocity profile given by Equation (45), we can utilize Equation (40) to rewrite the transport equation (Equation (35)) in the form

	$\displaystyle u\,\frac{\partial f_{{}_{\rm G}}}{\partial\tau}$	$\displaystyle=$	$\displaystyle-\alpha\,\frac{\chi}{3}\,\frac{\partial f_{{}_{\rm G}}}{\partial\chi}+\frac{1}{3y^{2}}\,\frac{\partial}{\partial\tau}\left(y^{2}\,\frac{\partial f_{{}_{\rm G}}}{\partial\tau}\right)-\xi^{2}u^{2}y^{2}f_{{}_{\rm G}}$		(54)
		$\displaystyle+$	$\displaystyle\frac{\overline{\sigma}}{\sigma_{{}_{\|\|}}}\,\frac{kT_{e}}{m_{e}c^{2}}\frac{1}{\chi^{2}}\frac{\partial}{\partial\chi}\left[\chi^{4}\left(f_{{}_{\rm G}}+\frac{\partial f_{{}_{\rm G}}}{\partial\chi}\right)\right]+\frac{\dot{N}_{0}\,\delta(\chi-\chi_{0})\,\delta(\tau-\tau_{0})}{\Omega\,c\,R_{0}^{2}\chi_{0}^{2}(kT_{e})^{3}}\ .$		(54)

This equation is separable if we restrict attention to the case with $\chi\neq\chi_{0}$ , so that the source term vanishes. In this situation, we can reorganize the transport equation to obtain

	$\displaystyle\alpha\,\frac{\chi}{3}\,\frac{\partial f_{{}_{\rm G}}}{\partial\chi}$	$\displaystyle-$	$\displaystyle\frac{\overline{\sigma}}{\sigma_{{}_{\|\|}}}\,\frac{kT_{e}}{m_{e}c^{2}}\frac{1}{\chi^{2}}\frac{\partial}{\partial\chi}\left[\chi^{4}\left(f_{{}_{\rm G}}+\frac{\partial f_{{}_{\rm G}}}{\partial\chi}\right)\right]$		(55)
		$\displaystyle=$	$\displaystyle\frac{1}{3y^{2}}\,\frac{\partial}{\partial\tau}\left(y^{2}\,\frac{\partial f_{{}_{\rm G}}}{\partial\tau}\right)-u\,\frac{\partial f_{{}_{\rm G}}}{\partial\tau}-\xi^{2}u^{2}y^{2}f_{{}_{\rm G}}\ ,$		(55)

which has been rearranged so that all of the energy operators are all on the left-hand side and all of the spatial operators are on the right-hand side. Following BW07, Equation (55) can be separated in energy and space using the functions

f_{\lambda}(\tau,\chi)\equiv g(\lambda,\tau)\ h(\lambda,\chi)\ ,

(56)

where $\lambda$ is the separation constant and $g$ and $h$ denote the spatial and energy separation functions, respectively. Setting $f_{{}_{\rm G}}=f_{\lambda}=gh$ in Equation (55) yields

	$\displaystyle\frac{1}{h}\,\bigg{\{}\alpha\,\frac{\chi}{3}\,\frac{dh}{d\chi}$	$\displaystyle-$	$\displaystyle\frac{\overline{\sigma}}{\sigma_{{}_{\|\|}}}\,\frac{kT_{e}}{m_{e}c^{2}}\frac{1}{\chi^{2}}\frac{d}{d\chi}\left[\chi^{4}\left(h+\frac{dh}{d\chi}\right)\right]\bigg{\}}$		(57)
		$\displaystyle=$	$\displaystyle\frac{1}{g}\,\bigg{[}\frac{1}{3y^{2}}\,\frac{d}{d\tau}\left(y^{2}\,\frac{dg}{d\tau}\right)-u\,\frac{dg}{d\tau}-\xi^{2}u^{2}y^{2}g\bigg{]}=-\frac{\alpha}{3}\lambda\ ,$		(57)

from which we conclude that the functions $g$ and $h$ satisfy the differential equations

\frac{1}{3y^{2}}\,\frac{d}{d\tau}\left(y^{2}\,\frac{dg}{d\tau}\right)-u\,\frac{dg}{d\tau}+\left(\frac{\alpha\lambda}{3}-\xi^{2}u^{2}y^{2}\right)g=0\ ,

(58)

and

\frac{1}{\chi^{2}}\frac{d}{d\chi}\left[\chi^{4}\left(h+\frac{dh}{d\chi}\right)\right]-\psi\,\chi\,\frac{dh}{d\chi}-\psi\,\lambda\,h=0\ ,

(59)

respectively, where the constant $\psi$ is defined by

\psi\equiv\frac{\alpha}{3}\,\frac{\sigma_{{}_{||}}}{\overline{\sigma}}\,\frac{m_{e}c^{2}}{kT_{e}}\ ,

(60)

which is equivalent to Equation (38) from BW07. By combining Equations (44) and (60), we can rewrite the dimensionless constant $\psi$ as

\psi=\frac{k_{\infty}^{2}\Omega}{4\pi}\left(\frac{\dot{M}}{\dot{M}_{\rm E}}\right)^{-1}\left(\frac{\overline{\sigma}}{\sigma_{{}_{\rm T}}}\right)^{-1}\,\left(\frac{kT_{e}}{m_{e}c^{2}}\right)^{-1}\ ,

(61)

where the Eddington accretion rate $\dot{M}_{\rm E}$ is defined in Equation (39).

4.1 Spatial Boundary Conditions

The spatial separation functions $g(\lambda,\tau)$ , must satisfy Equation (58) in combination with suitable boundary conditions. There are a variety of boundary conditions that can be employed in the application to X-ray pulsar accretion columns, which are based on placing various constraints on the radiation flux at the stellar surface and at the top of the accretion column. We discuss the relevant formulations here, and motivate the specific choice we make in carrying out the derivation of the Green’s function $f_{{}_{\rm G}}$ obtained in this paper.

The spatial boundary conditions are based fundamentally on the behavior of the specific radiation flux, $F$ , also called the “streaming function,” which is related to the radiation distribution function, $f$ , via the expression (e.g., Becker, 1992)

F(R,\epsilon)\equiv-\frac{\epsilon\upsilon}{3}\,\frac{\partial f}{\partial\epsilon}-\kappa\,\frac{\partial f}{\partial R}\quad\propto\quad{\rm s^{-1}\,cm^{-2}\,erg^{-3}}\ ,

(62)

where the spatial diffusion coefficient, $\kappa$ , is given by

\kappa(R)=\frac{c}{3\,\sigma_{{}_{||}}n_{e}(R)}\ .

(63)

The first and second terms on the right-hand side of Equation (62) correspond to advection and spatial diffusion, respectively. The corresponding total number flux of the radiation, $F_{\#}$ , is computed from $F$ using the energy integration

F_{\#}(R)=\int_{0}^{\infty}\epsilon^{2}F(R,\epsilon)\,d\epsilon\quad\propto\quad{\rm s^{-1}\,cm^{-2}}\ ,

(64)

which yields

F_{\#}(R)=\upsilon\,n_{r}-\kappa\,\frac{\partial n_{r}}{\partial R}\ ,

(65)

where $n_{r}$ is the radiation number density defined in Equation (22). Likewise, the total radiation energy flux, $F_{\rm en}$ , is computed using the integral

F_{\rm en}(R)=\int_{0}^{\infty}\epsilon^{3}F(R,\epsilon)\,d\epsilon\quad\propto\quad{\rm ergs\,s^{-1}\,cm^{-2}}\ ,

(66)

from which we obtain

F_{\rm en}(R)=\frac{4}{3}\,\upsilon\,U_{r}-\kappa\,\frac{\partial U_{r}}{\partial R}\ ,

(67)

where $U_{r}$ is the radiation energy density, evaluated using Equation (22). At the top of the accretion column, we will utilize a free-streaming boundary condition to ensure that spatial diffusion inside the column makes a proper transition when the accreting gas becomes optically thin to electron scattering. At the lower surface of the accretion column, where the gas merges with the neutron star crust, there are three primary scenarios for imposing boundary conditions on the radiation flux, as discussed below.

4.1.1 Free-Streaming Upper Boundary Condition

At the upper surface of the accretion column, the gas becomes optically thin to electron scattering, and therefore the radiation transport makes a transition from spatial diffusion (a three-dimensional random walk) to radial free streaming of the escaping photons at the speed of light. The free-streaming boundary condition applicable at the top of the accretion column can be written as

-\kappa\,\frac{\partial f}{\partial R}=cf\ ,\qquad R\to R_{\rm top}\ ,

(68)

where $R_{\rm top}$ denotes the radius at the top of the column, which is a free parameter in our model. By employing Equations (31) and (63), we can transform the spatial variable using

-\kappa\,\frac{\partial}{\partial R}=-\frac{c}{3}\,\frac{\partial}{\partial\tau}\ ,

(69)

which allows us to rewrite the free-streaming upper boundary condition as

-\frac{1}{3}\,\frac{\partial f}{\partial\tau}=f\ ,\qquad\tau\to\tau_{\rm top}\ .

(70)

Here, $\tau_{\rm top}$ denotes the scattering optical depth at the top of the accretion column, which is related to $R_{\rm top}$ via (see Equation (52))

\tau_{\rm top}=\tau(y_{\rm top})\ ,

(71)

where $y_{\rm top}\equiv R_{\rm top}/R_{*}$ is the dimensionless radius at the top of the accretion column (see Equation (32)).

4.1.2 Zero Diffusion Flux Lower Boundary Condition

There are three possible prescriptions for setting the radiation boundary condition imposed at the lower surface of the accretion column, where the gas merges with the crust of the neutron star. The simplest boundary condition is to set the diffusive flux equal to zero at the stellar surface, i.e.,

-\kappa\,\frac{\partial f}{\partial R}=0\ ,\qquad R\to R_{*}\ ,

(72)

which can be expressed in terms of the optical depth as

-\frac{1}{3}\,\frac{\partial f}{\partial\tau}=0\ ,\qquad\tau\to 0\ .

(73)

This boundary condition is applicable in situations in which the flow smoothly decelerates to rest at the base of the accretion column, as expected in luminous sources such as Her X-1, in which the gas passes through a radiative, radiation-dominated shock before setting onto the stellar surface (BW07). However, in low-luminosity sources such as X Per, this boundary condition may not be applicable, since the gas is expected to strike the stellar surface with a substantial residual velocity due to the lack of sufficient radiation pressure to accomplish smooth deceleration. In this case, the final merger occurs via Coulombic deceleration (Sokolova-Lapa et al., 2021).

4.1.3 Zero Number Flux Lower Boundary Condition

The second possibility is that the net flux of the photon number vanishes at the lower boundary of the accretion column, which is the photon-conservation scenario. In this case, we can use see Equation (65) to write the radiation boundary condition at the surface of the neutron star as

F_{\#}(R)=\upsilon\,n_{r}-\kappa\,\frac{\partial n_{r}}{\partial R}=0\ ,\qquad R\to R_{*}\ ,

(74)

or, in terms of the optical depth,

F_{\#}(\tau)=\upsilon\,n_{r}-\frac{c}{3}\,\frac{\partial n_{r}}{\partial\tau}=0\ ,\qquad\tau\to 0\ .

(75)

The advantage of this boundary condition is that photon conservation is guaranteed at the stellar surface. However, unfortunately, this condition does not guarantee the conservation of the photon energy flux at the surface of the neutron star. That problem leads to the third possibility for the spatial boundary condition at the stellar surface.

4.1.4 Zero Energy Flux Lower Boundary Condition

The third alternative for the spatial boundary condition at the lower boundary of the accretion column is based on the requirement that the net flux of the photon energy vanishes at the stellar surface, which is an energy-conservation principle. In this scenario, we can use Equation (67) to write the radiation boundary condition at the surface of the neutron star in the form

F_{\rm en}(R)=\frac{4}{3}\,\upsilon\,U_{r}-\kappa\,\frac{\partial U_{r}}{\partial R}=0\ ,\qquad R\to R_{*}\ ,

(76)

or, in terms of the optical depth,

F_{\rm en}(\tau)=\frac{4}{3}\,\upsilon\,U_{r}-\frac{c}{3}\,\frac{\partial U_{r}}{\partial\tau}=0\ ,\qquad\tau\to 0\ .

(77)

This condition ensures that there is no net photon energy flux either into or out of the star at the stellar surface, which establishes energy conservation at the base of the accretion column. In practical terms, this means that the luminosity of the radiated X-rays is strictly connected with processes occurring inside the accretion column, and there is no anomalous source of energy emitted from the interior of the neutron star, so that the emergent X-ray luminosity, $L_{X}$ , is related to the accretion rate, $\dot{M}$ , via

L_{X}=\frac{GM_{*}\dot{M}}{R_{*}}\ .

(78)

In our view, this is the overarching requirement of any physical model attempting to explain the formation of the X-ray continuum spectrum in an X-ray pulsar accretion column. Hence we will utilize Equation (77) as the fundamental spatial boundary condition for the radiation field at the stellar surface as we develop our solution for the emitted X-ray spectrum.

4.2 Eigenvalues and Spatial Eigenfunctions

The solution for the Green’s function, $f_{{}_{\rm G}}$ , can be expressed as the infinite series

f_{{}_{\rm G}}(\tau_{0},\chi_{0},\tau,\chi)=\sum_{n=0}^{\infty}\ C_{n}\,g_{n}(\tau)\,h_{n}(\chi)\ ,

(79)

where the expansion coefficients are denoted by $C_{n}$ , and the spatial eigenfunctions, $g_{n}(\tau)$ , are defined by

g_{n}(\tau)\equiv g(\lambda_{n},\tau)\ ,

(80)

with $\lambda_{n}$ denoting an eigenvalue of the separation constant $\lambda$ . Equation (79) is valid provided the spatial eigenfunctions $g_{n}(\tau)$ form an orthogonal set, which we will establish in Section 4.3. The eigenvalues $\lambda_{n}$ are determined by requiring that the spatial separation functions $g_{n}$ satisfy appropriate physical boundary conditions at the surface of the neutron star ( $\tau=0$ ), and at the top of the accretion column ( $\tau=\tau_{\rm top}$ ). We discuss the development of the upper and lower boundary conditions for the function $g_{n}$ below.

4.2.1 Upper Boundary Condition

At the top of the accretion column, located at radius $R=R_{\rm top}$ and scattering optical depth $\tau=\tau_{\rm top}$ , the radiation distribution must satisfy the free-streaming boundary condition given by Equation (70), which can be rewritten in terms of the Green’s function $f_{{}_{\rm G}}$ as

\frac{1}{3}\,\frac{\partial f_{{}_{\rm G}}}{\partial\tau}+f_{{}_{\rm G}}=0\ ,\qquad\tau\to\tau_{\rm top}\ .

(81)

We can use this expression to derive a related boundary condition for the spatial separation function $g_{n}$ . By substituting Equation (79) into Equation (81) and interchanging the order of summation and differentiation, we obtain

\sum_{n=0}^{\infty}\ C_{n}\,h_{n}(\chi)\left(\frac{1}{3}\,\frac{dg_{n}}{d\tau}+g_{n}\right)=0\ ,\qquad\tau\to\tau_{\rm top}\ .

(82)

This result implies that at the top of the accretion column, the spatial separation functions $g_{n}$ must satisfy the boundary condition

\frac{1}{3}\,\frac{dg_{n}}{\partial\tau}+g_{n}=0\ ,\qquad\tau\to\tau_{\rm top}\ .

(83)

The satisfaction of Equation (83) by the spatial separation functions $g_{n}$ is a necessary and sufficient condition to ensure that the Green’s function, $f_{{}_{\rm G}}$ , satisfies the free-streaming boundary condition at the top of the accretion column expressed by Equation (81).

4.2.2 Lower Boundary Condition

We can use the vanishing of the net photon energy flux at the stellar surface (expressed by Equation (77)) as the starting point for developing the appropriate physical boundary condition applicable to the spatial separation function $g_{n}$ at $\tau=0$ . By operating on Equation (79) with $\int_{0}^{\infty}\epsilon^{3}\,d\epsilon$ and integrating term-by-term, we can show that the solution for the radiation energy density, $U_{r}$ (Equation (22)), is given by the series

U_{r}(\tau)=(kT_{e})^{4}\sum_{n=0}^{\infty}\ C_{n}\,g_{n}(\tau)\,\int_{0}^{\infty}\chi^{3}h_{n}(\chi)\,d\chi\ ,

(84)

where we have transformed from $\epsilon$ to $\chi$ using Equation (30). Next, we can substitute Equation (84) into Equation (77) to show that the radiation energy flux is given by

F_{\rm en}(\tau)=(kT_{e})^{4}\sum_{n=0}^{\infty}\ C_{n}\,\left(\frac{4}{3}\,\upsilon\,g_{n}-\frac{c}{3}\,\frac{dg_{n}}{d\tau}\right)\,\int_{0}^{\infty}\chi^{3}h_{n}(\chi)\,d\chi\ .

(85)

According to Equation (77), the radiation energy flux vanishes at the stellar surface, and therefore $F_{\rm en}(\tau)\to 0$ as $\tau\to 0$ . In combination with Equation (85), this condition implies that the spatial separation functions $g_{n}$ must satisfy the boundary condition

\frac{4}{3}\,\upsilon\,g_{n}-\frac{c}{3}\,\frac{dg_{n}}{d\tau}=0\ ,\qquad\tau\to 0\ .

(86)

The satisfaction of Equation (86) by the spatial separation functions $g$ is a necessary and sufficient condition to guarantee that the Green’s function $f_{{}_{\rm G}}$ complies with the zero-energy-flux boundary condition at the bottom of the accretion column expressed by Equation (77).

4.2.3 Spatial Eigenfunctions

The global eigenfunctions $g_{n}(\tau)$ are those functions that satisfy the fundamental differential equation given by Equation (58), in combination with the upper and lower boundary conditions, expressed by Equations (83) and Equation (86), respectively. Equation (86) ensures that the radiation energy flux vanishes at the base of the accretion column, and Equation (83) ensures that the radiation diffusion flux correctly transitions to the proper free-streaming form at the top of the accretion column. This combination of relations is only satisfied when the separation constant $\lambda$ equals one of the discrete eigenvalues $\lambda_{n}$ .

The general procedure required to solve for the eigenvalues $\lambda_{n}$ and the associated eigenfunctions $g_{n}(\tau)$ involves the bi-directional integration of Equation (58). The first step is to choose a candidate value for $\lambda$ , which we wish to investigate in order to determine whether or not it is an eigenvalue. The second step is to integrate numerically Equation (58) starting at the stellar surface, $\tau=0$ , and utilizing the lower boundary condition given by Equation (86). The solution obtained in this step is referred to as $g_{\rm lower}(\lambda,\tau)$ . The third step is to integrate numerically Equation (58) starting at the top of the accretion column, $\tau=\tau_{\rm top}$ , and utilizing the upper boundary condition expressed by Equation (83). The solution obtained in this step is referred to as $g_{\rm upper}(\lambda,\tau)$ . Because $g_{\rm lower}(\lambda,\tau)$ and $g_{\rm upper}(\lambda,\tau)$ are continuous functions of $\tau$ we can confirm that the candidate value of $\lambda$ is an eigenvalue by establishing that these two functions are linearly dependent, which requires evaluation of the Wronskian of the two functions, ${\cal W}$ , defined by

{\cal W}(\lambda,\tau)\equiv g_{\rm lower}(\lambda,\tau)\,g_{\rm upper}^{\prime}(\lambda,\tau)-g_{\rm upper}(\lambda,\tau)\,g_{\rm lower}^{\prime}(\lambda,\tau)\ ,

(87)

where primes denote differentiation with respect to $\tau$ . The two functions $g_{\rm lower}$ and $g_{\rm upper}$ are linearly dependent if the Wronskian ${\cal W}$ vanishes. Hence, the eigenvalues $\lambda_{n}$ are the roots of the equation

{\cal W}(\lambda,\tau_{*})=0\ ,

(88)

where $\tau_{*}$ is an arbitrary value of $\tau$ within the computational domain, so that $\tau_{\rm top}>\tau_{*}>0$ . When the separation constant $\lambda$ is equal to one of the eigenvalues $\lambda_{n}$ , then the resulting functions $g_{\rm upper}(\tau)$ and $g_{\rm lower}(\tau)$ describe the same global eigenfunction, $g_{n}(\tau)$ , which is a smooth continuous function throughout the computational domain.

4.3 Eigenfunction Orthogonality

In order to evaluate the Green’s function, $f_{{}_{\rm G}}$ , using the series expansion in Equation (79), it is essential that we establish the orthogonality of the spatial eigenfunctions $g_{n}(\tau)$ . We begin by noting that the general Sturm-Liouville operator can be written as

\frac{d}{d\tau}\left[S(\tau)\frac{dg_{n}}{d\tau}\right]+\alpha\lambda_{n}\omega(\tau)g_{n}(\tau)+V(\tau)g_{n}(\tau)=0\ ,

(89)

where $\omega(\tau)$ is the weight function and $S(\tau)$ and $V(\tau)$ are auxiliary functions. Expansion and reorganization of Equation (89) yields the equivalent differential equation

g_{n}^{\prime\prime}(\tau)+\frac{S^{\prime}(\tau)}{S(\tau)}\,g_{n}^{\prime}(\tau)+\alpha\lambda_{n}\frac{\omega(\tau)}{S(\tau)}\,g_{n}(\tau)+\frac{V(\tau)}{S(\tau)}\,g_{n}(\tau)=0\ ,

(90)

where primes denote differentiation with respect to $\tau$ . Next we note that the fundamental differential equation (Equation (58)) satisfied by the spatial separation functions $g(\tau)$ can be rewritten in the form

g_{n}^{\prime\prime}(\tau)+\left[\frac{2y^{\prime}(\tau)}{y(\tau)}-3u(\tau)\right]g_{n}^{\prime}(\tau)+\alpha\lambda_{n}\,g_{n}(\tau)-3\xi^{2}u^{2}(\tau)y^{2}(\tau)\,g_{n}(\tau)=0\ .

(91)

Comparison of Equations (90) and (91) yields the identifications

\omega(\tau)=S(\tau)=y^{2}(\tau)\exp\left[-3\int_{0}^{\tau}u(\tau^{\prime})\,d\tau^{\prime}\right]\ ,

(92)

and

V(\tau)=-3\xi^{2}u^{2}(\tau)y^{2}(\tau)S(\tau)\ ,

(93)

where $\omega(\tau)$ is the weight function.

We can make further progress by invoking the relationship between the scattering optical depth $\tau$ and the dimensionless radius $y$ . Based on Equation (49), we can write

\frac{d\tau}{dy}=\frac{3k_{\infty}^{2}R_{g}}{\alpha R_{*}y^{2}|u(y)|}\ ,

(94)

which can be utilized to transform the variable of integration in Equation (92) from $\tau$ to $y$ . After simplification, the result obtained for the weight function is

\omega(\tau)=S(\tau)=y^{2}(\tau)\exp\left\{\frac{9R_{g}k_{\infty}^{2}}{\alpha R_{*}}\left[1-\frac{1}{y(\tau)}\right]\right\}\ .

(95)

Likewise, Equation (93) now becomes

V(\tau)=-3\,\xi^{2}u^{2}(\tau)\,y^{4}(\tau)\exp\left\{\frac{9R_{g}k_{\infty}^{2}}{\alpha R_{*}}\left[1-\frac{1}{y(\tau)}\right]\right\}\ .

(96)

We are now in position to establish the orthogonality of the family of spatial eigenfunctions, $g_{n}(\tau)$ , using a standard analysis procedure based on the Sturm-Liouville form of the transport equation (Equation (89)). The first step is to replace $g_{n}$ with $g_{m}$ in Equation (89) to obtain

\frac{d}{d\tau}\left[S(\tau)\frac{dg_{m}}{d\tau}\right]+\alpha\lambda_{m}\omega(\tau)g_{m}(\tau)+V(\tau)g_{m}(\tau)=0\ .

(97)

Next we multiply Equation (89) by $g_{m}$ and Equation (97) by $g_{n}$ and subtract the first from the second, which yields

g_{n}(\tau)\frac{d}{d\tau}\left[S(\tau)\frac{dg_{m}}{d\tau}\right]-g_{m}(\tau)\frac{d}{d\tau}\left[S(\tau)\frac{dg_{n}}{d\tau}\right]=(\lambda_{n}-\lambda_{m})\,\alpha\,\omega(\tau)\,g_{n}(\tau)\,g_{m}(\tau)\ .

(98)

This relation can be integrated by parts from $\tau=0$ to $\tau=\tau_{\rm top}$ to obtain

S(\tau)\left[g_{n}(\tau)\frac{dg_{m}}{d\tau}-g_{m}(\tau)\frac{dg_{n}}{d\tau}\right]\bigg{|}_{0}^{\tau_{\rm top}}=(\lambda_{n}-\lambda_{m})\int_{0}^{\tau_{\rm top}}\alpha\,\omega(\tau)\,g_{n}(\tau)\,g_{m}(\tau)\,d\tau\ .

(99)

Based on the boundary conditions for $g(\tau)$ given by Equations (83) and (86), we conclude that in the two limits $\tau\to 0$ and $\tau\to\tau_{\rm top}$ , the left-hand side of Equation (99) vanishes, which therefore establishes the orthogonality relation

\int_{0}^{\tau_{\rm top}}\omega(\tau)\,g_{n}(\tau)\,g_{m}(\tau)\,d\tau=0,\qquad n\neq m\ .

(100)

This important result confirms that the family of spatial eigenfunctions $g_{n}(\tau)$ form an orthogonal set relative to the weight function $\omega(\tau)$ defined in Equation (95). This is an essential property for carrying out the series expansion for the Green’s function $f_{{}_{\rm G}}$ in Equation (79).

4.4 Energy Eigenfunctions

The energy separation functions, $h(\lambda,\chi)$ , are the solutions to Equation (59) that satisfy appropriate physical boundary conditions in the energy space. As $\chi\to 0$ , we require that $h$ not increase faster than $\chi^{-3}$ in order to ensure that the Green’s function possesses a finite total photon number density (see Equation (22)). Conversely, as $\chi\to\infty$ , we require that $h$ decrease more rapidly than $\chi^{-4}$ since the Green’s function must also contain a finite total photon energy density. We also require that $h$ be continuous at the injection energy $\chi=\chi_{0}$ in order to avoid an infinite diffusive flux in the energy space. The fundamental solution to Equation (59) that satisfies the various boundary and continuity conditions can be written as (BW07)

h(\lambda,\chi)=\begin{cases}\chi^{\kappa-4}\,e^{-\chi/2}\,W_{\kappa,\mu}(\chi_{0})\,M_{\kappa,\mu}(\chi)\ ,&\chi\leq\chi_{0}\ ,\\ \chi^{\kappa-4}\,e^{-\chi/2}\,M_{\kappa,\mu}(\chi_{0})\,W_{\kappa,\mu}(\chi)\ ,&\chi\geq\chi_{0}\ ,\\ \end{cases}

(101)

where $M_{\kappa,\mu}$ and $W_{\kappa,\mu}$ denote Whittaker functions (Abramowitz & Stegun, 1970), and we have made the definitions

\kappa\equiv\frac{1}{2}\,(\psi+4)\ ,\ \ \ \ \ \mu\equiv\frac{1}{2}\left[(3-\psi)^{2}+4\,\psi\lambda\right]^{1/2}\ ,

(102)

with the parameter $\psi$ defined in Equation (60).

When the separation constant $\lambda$ is equal to one of the eigenvalues $\lambda_{n}$ , then $h(\lambda,\chi)$ is an energy eigenfunction, denoted by $h_{n}(\chi)$ , which can be written in the compact form

h_{n}(\chi)\equiv h(\lambda_{n},\chi)=\chi^{\kappa-4}\,e^{-\chi/2}\,M_{\kappa,\mu}(\chi_{\rm min})\,W_{\kappa,\mu}(\chi_{\rm max})\ ,

(103)

where

\chi_{\rm min}\equiv\min(\chi,\chi_{0})\ ,\ \ \ \ \ \chi_{\rm max}\equiv\max(\chi,\chi_{0})\ .

(104)

We note that each of the eigenvalues $\lambda_{n}$ results in a different value for $\mu$ by setting $\lambda=\lambda_{n}$ in Equation (102).

4.5 Green’s Function Expansion

Equation (79) gives the series solution for the Green’s function, $f_{{}_{\rm G}}$ , as a function of the injection optical depth, $\tau_{0}$ , the injection dimensionless energy, $\chi_{0}$ , the evaluation optical depth, $\tau$ , and the evaluation dimensionless energy, $\chi$ . The computation of the expansion coefficients, $C_{n}$ , in Equation (79) is accomplished by exploiting the orthogonality of the spatial eigenfunctions, $g_{n}(\tau)$ , along with the derivative jump condition

\lim_{\varepsilon\to 0}\ \bigg{\{}\frac{\partial f_{{}_{\rm G}}}{\partial\chi}\Bigg{|}_{\chi=\chi_{0}+\varepsilon}-\frac{\partial f_{{}_{\rm G}}}{\partial\chi}\Bigg{|}_{\chi=\chi_{0}-\varepsilon}\ \bigg{\}}=-\,\frac{3\,\psi\dot{N}_{0}\,\delta(\tau-\tau_{0})}{\alpha\,\Omega\,c\,R_{*}^{2}\,y_{0}^{2}\,\chi_{0}^{4}\,(kT_{e})^{3}}\ ,

(105)

which is obtained by integrating Equation (54) with respect to $\chi$ in a very small range surrounding the injection energy $\chi_{0}$ , and we have introduced the dimensionless injection radius, $y_{0}\equiv R_{0}/R_{*}$ . The parameter $\psi$ appearing in Equation (105) is defined in Equation (60). Combining Equations (79) and (105) yields

\lim_{\varepsilon\to 0}\ \sum_{n=0}^{\infty}\ C_{n}\,g_{n}(\tau)\,\left[h_{n}^{\prime}(\chi_{0}+\varepsilon)-h_{n}^{\prime}(\chi_{0}-\varepsilon)\right]=-\,\frac{3\,\psi\dot{N}_{0}\,\delta(\tau-\tau_{0})}{\alpha\,\Omega\,c\,R_{*}^{2}\,y_{0}^{2}\,\chi_{0}^{4}\,(kT_{e})^{3}}\ ,

(106)

where primes denote differentiation with respect to $\chi$ . By using Equation (103) to substitute for $h_{n}(\chi)$ , we find that

\sum_{n=0}^{\infty}\ C_{n}\,g_{n}(\tau)\,\chi_{0}^{\kappa-4}\,e^{-\chi_{0}/2}\mathfrak{W}(\chi_{0})=-\,\frac{3\,\psi\dot{N}_{0}\,\delta(\tau-\tau_{0})}{\alpha\,\Omega\,c\,R_{*}^{2}\,y_{0}^{2}\,\chi_{0}^{4}\,(kT_{e})^{3}}\ ,

(107)

where $\mathfrak{W}$ denotes the Wronskian of the Whittaker functions, defined by

\mathfrak{W}(\chi_{0})\equiv M_{\kappa,\mu}(\chi_{0})\,W^{\prime}_{\kappa,\mu}(\chi_{0})-W_{\kappa,\mu}(\chi_{0})\,M^{\prime}_{\kappa,\mu}(\chi_{0})\ .

(108)

The Wronskian can be evaluated analytically using Equation (55) from BW07, which yields

\mathfrak{W}(\chi_{0})=-\,\frac{\Gamma(1+2\mu)}{\Gamma(\mu-\kappa+1/2)}\ .

(109)

Using Equation (109) to substitute for $\mathfrak{W}(\chi_{0})$ in Equation (107), and rearranging factors, we find that

\sum_{n=0}^{\infty}\ \frac{\Gamma(1+2\mu)\,C_{n}\,g_{n}(\tau)}{\Gamma(\mu-\kappa+1/2)}=\,\frac{3\,\psi\dot{N}_{0}e^{\chi_{0}/2}\,\delta(\tau-\tau_{0})}{\alpha\,\Omega\,c\,R_{*}^{2}\,y_{0}^{2}\,\chi_{0}^{\kappa}\,(kT_{e})^{3}}\ ,

(110)

where $\mu$ is a function of $\lambda=\lambda_{n}$ according to Equation (102). We can derive an expression for the expansion coefficients $C_{n}$ by exploiting the orthogonality of the spatial eigenfunctions $g_{n}(\tau)$ . Multiplying both sides of Equation (110) by the weight function $\omega(\tau)$ and the eigenfunction $g_{m}(\tau)$ and integrating with respect to $\tau$ from $\tau=0$ to $\tau=\tau_{\rm top}$ , we find that, according to the orthogonality relation in Equation (100), only the $n=m$ term in the sum survives, and we therefore obtain

C_{n}=\frac{3\,\psi\dot{N}_{0}e^{\chi_{0}/2}\,\omega(\tau_{0})}{\alpha\,\Omega\,c\,R_{*}^{2}\,y_{0}^{2}\,\chi_{0}^{\kappa}\,(kT_{e})^{3}}\,\frac{\Gamma(\mu-\kappa+1/2)\,g_{n}(\tau_{0})}{{\cal I}_{n}\,\Gamma(1+2\mu)}\ ,

(111)

where we have introduced the quadratic normalization integrals, ${\cal I}_{n}$ , defined by

{\cal I}_{n}\equiv\int_{0}^{\tau_{\rm top}}g_{n}^{2}(\tau)\omega(\tau)\,d\tau\ .

(112)

The final closed-form solution for the Green’s function is obtained by combining Equations (79), (103), and (111), which yields

	$\displaystyle f_{{}_{\rm G}}(\tau_{0},\tau,\chi_{0},\chi)$	$\displaystyle=$	$\displaystyle\frac{3\,\psi\dot{N}_{0}\,\omega(\tau_{0})\,\chi^{\kappa-4}\,e^{(\chi_{0}-\chi)/2}}{\alpha\,\Omega\,c\,R_{*}^{2}\,y_{0}^{2}\,\chi_{0}^{\kappa}\,(kT_{e})^{3}}\sum_{n=0}^{\infty}\ \frac{\Gamma(\mu-\kappa+1/2)}{{\cal I}_{n}\,\Gamma(1+2\mu)}$		(113)
		$\displaystyle\times$	$\displaystyle g_{n}(\tau_{0})\,g_{n}(\tau)\,M_{\kappa,\mu}(\chi_{\rm min})\,W_{\kappa,\mu}(\chi_{\rm max})\ ,$		(113)

where the parameters $\kappa$ and $\mu$ are computed using Equation (102), and $\chi_{\rm min}$ and $\chi_{\rm max}$ are defined in Equation (104). The Green’s function can also be expressed directly in terms of the photon energy $\epsilon$ by writing

	$\displaystyle f_{{}_{\rm G}}(\tau_{0},\tau,\epsilon_{0},\epsilon)$	$\displaystyle=$	$\displaystyle\frac{3\,\psi\dot{N}_{0}\,\omega(\tau_{0})\,kT_{e}\,\epsilon^{\kappa-4}\,e^{(\epsilon_{0}-\epsilon)/(2kT_{e})}}{\alpha\,\Omega\,c\,R_{*}^{2}\,y_{0}^{2}\,\epsilon_{0}^{\kappa}}\sum_{n=0}^{\infty}\ \frac{\Gamma(\mu-\kappa+1/2)}{{\cal I}_{n}\,\Gamma(1+2\mu)}$		(114)
		$\displaystyle\times$	$\displaystyle g_{n}(\tau_{0})\,g_{n}(\tau)\,M_{\kappa,\mu}\left(\frac{\epsilon_{\rm min}}{kT_{e}}\right)W_{\kappa,\mu}\left(\frac{\epsilon_{\rm max}}{kT_{e}}\right)\ ,$		(114)

where

\epsilon_{\rm min}\equiv\min(\epsilon,\epsilon_{0})\ ,\ \ \ \ \ \epsilon_{\rm max}\equiv\max(\epsilon,\epsilon_{0})\ .

(115)

Equations (113) and (114) provide the exact solution for the steady-state Green’s function, $f_{{}_{\rm G}}$ , which represents the radiation spectrum inside the accretion column, for selected values of the photon energy $\epsilon$ and the scattering optical depth $\tau$ , resulting from the continual injection of $\dot{N}_{0}$ monochromatic seed photons per unit time, each with energy $\epsilon_{0}$ , from a source located at optical depth $\tau_{0}$ . The series expansions in Equations (113) and (114) converge rapidly, and we generally obtain three significant figures of accuracy in our calculations of $f_{{}_{\rm G}}$ using the first 5-10 terms in the series.

4.6 Radiation Distribution for Arbitrary Source Term

The fundamental transport equation (Equation (25)) is linear, and therefore the exact solution obtained for the Green’s function $f_{{}_{\rm G}}$ in Section 4.5 can be used to compute the particular solution, $f$ , representing the radiation distribution inside the accretion column resulting from an arbitrary source term $Q$ in the fundamental transport equation (Equation (21)). The particular solution for the radiation distribution, $f$ , associated with an arbitrary source function $Q$ is given by the integral convolution (Becker, 2003)

f(R,\epsilon)=\int_{R_{*}}^{R_{\rm top}}\int_{0}^{\infty}\frac{1}{\dot{N}_{0}}\,f_{{}_{\rm G}}(R_{0},R,\epsilon_{0},\epsilon)\,\Omega R_{0}^{2}\,Q(R_{0},\epsilon_{0})\,d\epsilon_{0}\,dR_{0}\ ,

(116)

where the Green’s function, $f_{{}_{\rm G}}$ , is computed using Equation (114), and the source function $Q$ is normalized so that the quantity $\epsilon_{0}^{2}\,\Omega R_{0}^{2}\,Q(R_{0},\epsilon_{0})\,d\epsilon_{0}\,dR_{0}$ gives the number of seed photons injected per unit time with energy between $\epsilon_{0}$ and $\epsilon_{0}+d\epsilon_{0}$ from the radial zone between $R_{0}$ and $R_{0}+dR_{0}$ . Examples of source functions relevant for accretion-powered X-ray pulsars include bremsstrahlung, cyclotron, and blackbody emission.

Equation (114) expresses the Green’s function, $f_{{}_{\rm G}}$ , in terms of the scattering optical depth, rather than the radius, and therefore it is convenient to transform the spatial variable of integration in Equation (116) from $R_{0}$ to $\tau_{0}$ using the expression (cf. Equation (49))

dR_{0}=R_{*}\,dy_{0}=\frac{\alpha R_{*}^{2}y_{0}^{2}|u|}{3k_{\infty}^{2}R_{g}}\,d\tau_{0}\ ,

(117)

which yields for the particular solution

f(\tau,\epsilon)=\int_{0}^{\tau_{\rm top}}\int_{0}^{\infty}\frac{1}{\dot{N}_{0}}\,f_{{}_{\rm G}}(\tau_{0},\tau,\epsilon_{0},\epsilon)\,Q(\tau_{0},\epsilon_{0})\,\frac{\alpha\Omega R_{*}^{4}\,y_{0}^{4}|u|}{3k_{\infty}^{2}R_{g}}\,d\epsilon_{0}\,d\tau_{0}\ ,

(118)

where $y_{0}$ is related to $\tau_{0}$ via Equation (52). Equation (118) gives the particular solution for the radiation distribution inside the accretion column, $f$ , resulting from the continual emission of photons from a source with distribution $Q$ . In Section 7 we will use Equation (118) to develop particular solutions for the radiation distribution $f$ inside the accretion column resulting from a variety of physical emission mechanisms.

5 EMITTED RADIATION SPECTRUM

Our transport equation formalism includes both an escape-probability formalism to model the diffusion of radiation through the column walls (Equation (29)), as well as a free-streaming boundary condition to model the escape of radiation through the top of the accretion column (Equation (68)). Since these two processes are included in our formalism, we can therefore self-consistently compute the radiation emission components due to the escape of photons through either the walls or the top of the column. The capability to compute these two emission components separately facilitates more detailed simulations of X-ray pulsar phase-averaged spectra than are possible using the formalism of BW07, who only considered the escape of radiation through the column walls. Phase-averaged spectra evaluated using the new two-component model are presented in Section 8. Here we discuss the specific procedures used to compute the two separate emission components.

5.1 Green’s Function for Column Wall Spectrum

In the escape-probability approach employed here, the Green’s function describing the number spectrum of the photons escaping through the walls of the conical accretion column is given by

\dot{N}_{\epsilon}^{\rm G}(R_{0},R,\epsilon_{0},\epsilon)\equiv\frac{\Omega\,R^{2}\epsilon^{2}}{t_{\rm esc}(\tau)}\,f_{{}_{\rm G}}(R_{0},R,\epsilon_{0},\epsilon)\ ,

(119)

where the escape timescale, $t_{\rm esc}$ , is evaluated using Equation (29). The quantity $\dot{N}_{\epsilon}^{\rm G}\,dR\,d\epsilon$ represents the number of photons escaping per unit time, with energy between $\epsilon$ and $\epsilon+d\epsilon$ , from the disk-shaped volume between radii $R$ and $R+dR$ . We can combine Equations (29), (37), and (43) to show that the escape timescale, $t_{\rm esc}$ , can also be written in the form

t_{\rm esc}(R)=\frac{\alpha R_{*}^{2}}{3k_{\infty}^{2}\xi^{2}R_{g}|\upsilon|}\ ,

(120)

which can be used to substitute for $t_{\rm esc}$ in Equation (119), yielding

\dot{N}_{\epsilon}^{\rm G}(R_{0},R,\epsilon_{0},\epsilon)=\Omega\,R^{2}\epsilon^{2}\,\frac{3k_{\infty}^{2}\xi^{2}R_{g}|\upsilon|}{\alpha R_{*}^{2}}\,f_{{}_{\rm G}}(R_{0},R,\epsilon_{0},\epsilon)\ .

(121)

Transforming the spatial variable from the radius $R$ to the scattering optical depth $\tau$ , we can combine Equations (114) and (121) to express the exact solution for Green’s function describing the photon number spectrum escaping through the column walls as

	$\displaystyle\dot{N}_{\epsilon}^{\rm G}(\tau_{0},\tau,\epsilon_{0},\epsilon)$	$\displaystyle=$	$\displaystyle\frac{9k_{\infty}^{2}y^{2}\xi^{2}R_{g}\|\upsilon\|\psi\dot{N}_{0}\,\omega(\tau_{0})\,kT_{e}\,\epsilon^{\kappa-2}\,e^{(\epsilon_{0}-\epsilon)/(2kT_{e})}}{\alpha^{2}cR_{*}^{2}\,y_{0}^{2}\,\epsilon_{0}^{\kappa}}$		(122)
		$\displaystyle\times$	$\displaystyle\sum_{n=0}^{\infty}\ \frac{\Gamma(\mu-\kappa+1/2)}{{\cal I}_{n}\,\Gamma(1+2\mu)}g_{n}(\tau_{0})\,g_{n}(\tau)\,M_{\kappa,\mu}\left(\frac{\epsilon_{\rm min}}{kT_{e}}\right)W_{\kappa,\mu}\left(\frac{\epsilon_{\rm max}}{kT_{e}}\right)\ ,$		(122)

where $\epsilon_{\rm min}$ and $\epsilon_{\rm max}$ are defined in Equations (115), and the ${\cal I}_{n}$ integrals are computed using Equation (112).

In many situations, we do not need to consider the detailed distribution of photons escaping through the column walls as a function of $R$ , and instead it is sufficient to consider the column-integrated Green’s function, $\Phi_{\epsilon}^{\rm G}$ , which describes the number distribution of the radiation escaping through the walls of the entire conical accretion column, from the stellar surface at $R=R_{*}$ up to the top of the column at $R=R_{\rm top}$ . The function $\Phi_{\epsilon}^{\rm G}$ is related to the distribution $\dot{N}_{\epsilon}^{\rm G}$ via

\Phi_{\epsilon}^{\rm G}(R_{0},\epsilon_{0},\epsilon)\equiv\int_{R_{*}}^{R_{\rm top}}\dot{N}_{\epsilon}^{\rm G}(R_{0},R,\epsilon_{0},\epsilon)\,dR\ ,

(123)

and the quantity $\Phi_{\epsilon}^{\rm G}\,d\epsilon$ gives the number of photons escaping through the column walls per unit time with energy between $\epsilon$ and $\epsilon+d\epsilon$ . Since Equation (122) gives $\dot{N}_{\epsilon}^{\rm G}$ as a function of the scattering optical depth $\tau$ , it is convenient to transform the variable of integration in Equation (123) from $R$ to $\tau$ using the expression (cf. Equation (117))

dR=R_{*}\,dy=\frac{\alpha R_{*}^{2}y^{2}|u|}{3k_{\infty}^{2}R_{g}}\,d\tau\ ,

(124)

which yields

\Phi_{\epsilon}^{\rm G}(\tau_{0},\epsilon_{0},\epsilon)\equiv\int_{0}^{\tau_{\rm top}}\dot{N}_{\epsilon}^{\rm G}(\tau_{0},\tau,\epsilon_{0},\epsilon)\,\frac{\alpha R_{*}^{2}y^{2}|u|}{3k_{\infty}^{2}R_{g}}\,d\tau\ ,

(125)

where $y$ is related to $\tau$ via Equation (52). Next, we employ Equation (122) to substitute for $\dot{N}_{\epsilon}^{\rm G}$ in Equation (125) and interchange the order of summation and integration to obtain

	$\displaystyle\Phi_{\epsilon}^{\rm G}(\tau_{0},\epsilon_{0},\epsilon)$	$\displaystyle=$	$\displaystyle\frac{6\psi\xi^{2}\dot{N}_{0}R_{g}\omega(\tau_{0})\,kT_{e}\,\epsilon^{\kappa-2}\,e^{(\epsilon_{0}-\epsilon)/(2kT_{e})}}{\alpha y_{0}^{2}R_{*}\epsilon_{0}^{\kappa}}$		(126)
		$\displaystyle\times$	$\displaystyle\sum_{n=0}^{\infty}\ \frac{\Gamma(\mu-\kappa+1/2)}{{\cal I}_{n}\,\Gamma(1+2\mu)}g_{n}(\tau_{0})\,I_{n}\,M_{\kappa,\mu}\left(\frac{\epsilon_{\rm min}}{kT_{e}}\right)W_{\kappa,\mu}\left(\frac{\epsilon_{\rm max}}{kT_{e}}\right)\ ,$		(126)

where the quantity $I_{n}$ denotes the spatial integral

I_{n}\equiv\left(\frac{R_{*}}{2R_{g}}\right)\int_{0}^{\tau_{\rm top}}g_{n}(\tau)\,y^{4}(\tau)\,u^{2}(\tau)\,d\tau\ ,

(127)

and the ${\cal I}_{n}$ integrals are defined in Equation (112). Equation (127) can be simplified by using Equation (45) to substitute for the dimensionless velocity, $u(\tau)$ , which yields

I_{n}=\int_{0}^{\tau_{\rm top}}g_{n}(\tau)\,[k_{0}^{2}+k_{\infty}^{2}(y^{3}-1)]\,d\tau\ .

(128)

Once the set of spatial eigenfunctions $g_{n}(\tau)$ has been obtained, the integrals in Equation (128) can be evaluated numerically.

5.2 Wall Spectrum for a General Source

Once the particular solution $f$ is determined using Equation (116), we can show by analogy with Equation (119) that the corresponding height-dependent number distribution describing the radiation spectrum escaping through the column walls for an arbitrary source term $Q$ can be computed using

\dot{N}_{\epsilon}(R,\epsilon)\equiv\frac{\Omega\,R^{2}\,\epsilon^{2}}{t_{\rm esc}(R)}\,f(R,\epsilon)\ .

(129)

By combining Equations (116), (119), and (129), we find that the particular solution for the emitted photon spectrum can be written as

\dot{N}_{\epsilon}(R,\epsilon)=\int_{R_{*}}^{R_{\rm top}}\int_{0}^{\infty}\frac{\dot{N}_{\epsilon}^{\rm G}(R_{0},R,\epsilon_{0},\epsilon)}{\dot{N}_{0}}\ \epsilon_{0}^{2}\,\Omega R_{0}^{2}\,Q(R_{0},\epsilon_{0})\,d\epsilon_{0}\,dR_{0}\ ,

(130)

where the quantity $\dot{N}_{\epsilon}\,dR\,d\epsilon$ represents the number of photons emitted per unit time with energy between $\epsilon$ and $\epsilon+d\epsilon$ from the disk-shaped volume between radii $R$ and $R+dR$ . It is convenient to transform the spatial variable of integration in Equation (130) from $R_{0}$ to $\tau_{0}$ using Equation (117), because $\dot{N}_{\epsilon}^{\rm G}$ is evaluated as a function of $(\tau_{0},\tau,\epsilon_{0},\epsilon)$ according to Equation (122). The change of variables yields for the height-dependent escaping photon number distribution the general expression

\dot{N}_{\epsilon}(R,\epsilon)=\int_{0}^{\tau_{\rm top}}\int_{0}^{\infty}\frac{\dot{N}_{\epsilon}^{\rm G}(\tau_{0},\tau,\epsilon_{0},\epsilon)}{\dot{N}_{0}}\ \epsilon_{0}^{2}\,Q(\tau_{0},\epsilon_{0})\,\frac{\alpha\Omega R_{*}^{4}\,y_{0}^{4}|u|}{3k_{\infty}^{2}R_{g}}\,d\epsilon_{0}\,d\tau_{0}\ ,

(131)

where the relationship between $y_{0}$ and $\tau_{0}$ is given by Equation (52). Equation (131) allows us to calculate the number distribution of the photons escaping through the walls of the accretion column at any radius $R$ as a function of the photon energy $\epsilon$ for an arbitrary source distribution $Q$ .

By analogy with Equation (123), we can also obtain the column-integrated photon spectrum escaping through the walls of the conical accretion column due to a general photon source $Q$ using the integral

\Phi_{\epsilon}(\epsilon)\equiv\int_{R_{*}}^{R_{\rm top}}\dot{N}_{\epsilon}(R,\epsilon)\,dR\ ,

(132)

where the quantity $\Phi_{\epsilon}\,d\epsilon$ gives the number of photons escaping through the column walls per unit time with energy between $\epsilon$ and $\epsilon+d\epsilon$ due to the photon source function $Q$ . By combining Equations (123), (130), and (132), we can show that the particular solution for the column-integrated photon number spectrum escaping through the column walls for an arbitrary source function $Q$ is given by

\Phi_{\epsilon}(\epsilon)=\int_{R_{*}}^{R_{\rm top}}\int_{0}^{\infty}\frac{\Phi_{\epsilon}^{\rm G}(R_{0},\epsilon_{0},\epsilon)}{\dot{N}_{0}}\ \epsilon_{0}^{2}\,\Omega R_{0}^{2}\,Q(R_{0},\epsilon_{0})\,d\epsilon_{0}\,dR_{0}\ .

(133)

Equation (126) gives $\Phi_{\epsilon}^{\rm G}$ as a function of $(\tau_{0},\epsilon_{0},\epsilon)$ , and therefore it is convenient to transform the spatial variable of integration in Equation (133) from $R_{0}$ to $\tau_{0}$ using Equation (117), which yields

\Phi_{\epsilon}(\epsilon)=\int_{0}^{\tau_{\rm top}}\int_{0}^{\infty}\frac{\Phi_{\epsilon}^{\rm G}(\tau_{0},\epsilon_{0},\epsilon)}{\dot{N}_{0}}\ \epsilon_{0}^{2}\,Q(\tau_{0},\epsilon_{0})\,\frac{\alpha\Omega R_{*}^{4}\,y_{0}^{4}|u|}{3k_{\infty}^{2}R_{g}}\,d\epsilon_{0}\,d\tau_{0}\ ,

(134)

where $y_{0}$ is related to $\tau_{0}$ via Equation (52). Equation (134) allows us to calculate the column-integrated photon number spectrum escaping through the walls of the accretion column, $\Phi_{\epsilon}$ , as a function of the photon energy, $\epsilon$ , for an arbitrary photon source function $Q$ . In Section 7 we will use Equation (134) to derive the particular solution for $\Phi_{\epsilon}$ associated with each of the primary photon emission mechanisms thought to be important in the formation of the observed spectra from accretion-powered X-ray pulsars.

5.3 Green’s Function for Column Top Spectrum

The escape of photons out of the top of the conical accretion column is mediated by the process of spatial diffusion, in which photons execute a random walk through the plasma by scattering off electrons. At the top of the column, the gas becomes optically thin, and the escape of radiation through the upper surface of the column is therefore described by the free-streaming boundary condition expressed by Equation (68). It follows that the Green’s function describing the number distribution of the photons diffusing through the upper surface of the accretion column, denoted b $\dot{\cal N}^{\,\rm G}_{\epsilon}$ , is given by

\dot{\cal N}^{\,\rm G}_{\epsilon}(\tau_{0},\epsilon_{0},\epsilon)=c\,\Omega\,R^{2}_{\rm top}\,\epsilon^{2}f_{{}_{\rm G}}(\tau_{0},\tau_{\rm top},\epsilon_{0},\epsilon)\ ,

(135)

where $R_{\rm top}$ is the radius at the column top, $\tau_{\rm top}$ denotes the scattering optical depth from the stellar surface to the column top, and the function $f_{{}_{\rm G}}$ is evaluated using Equation (114). The quantity $\dot{\cal N}^{\,\rm G}_{\epsilon}\,d\epsilon$ gives the number of photons with energy between $\epsilon$ and $\epsilon+d\epsilon$ escaping per unit time through the top of the accretion column. By combining Equation (114) and (135), we can show that the closed-form expression for $\dot{\cal N}^{\,\rm G}_{\epsilon}$ is given by

	$\displaystyle\dot{\cal N}^{\,\rm G}_{\epsilon}(\tau_{0},\epsilon_{0},\epsilon)$	$\displaystyle=$	$\displaystyle\frac{3\,\psi\dot{N}_{0}\,\omega(\tau_{0})\,kT_{e}\,R^{2}_{\rm top}\,\epsilon^{\kappa-2}\,e^{(\epsilon_{0}-\epsilon)/(2kT_{e})}}{\alpha R_{*}^{2}\,y_{0}^{2}\,\epsilon_{0}^{\kappa}}\sum_{n=0}^{\infty}\ \frac{\Gamma(\mu-\kappa+1/2)}{{\cal I}_{n}\,\Gamma(1+2\mu)}$		(136)
		$\displaystyle\times$	$\displaystyle g_{n}(\tau_{0})\,g_{n}(\tau_{\rm top})\,M_{\kappa,\mu}\left(\frac{\epsilon_{\rm min}}{kT_{e}}\right)W_{\kappa,\mu}\left(\frac{\epsilon_{\rm max}}{kT_{e}}\right)\ ,$		(136)

where $\epsilon_{\rm min}$ and $\epsilon_{\rm max}$ are evaluated using Equations (115). Equation (136) facilitates the computation of the Green’s function for the number distribution of the radiation escaping through the column top, resulting from the continual injection of seed photons with energy $\epsilon_{0}$ at optical depth $\tau_{0}$ . In the next section we will consider the convolution of the Green’s function with an arbitrary source term.

5.4 Column Top Spectrum for a General Source

In our applications to accretion-powered X-ray pulsars, we must consider the injection of seed photons via the processes of cyclotron, bremsstrahlung, and blackbody emission. By analogy with Equation (130) for the wall emission, we can show that the particular solution for the number distribution of the photons escaping through the top of the accretion column, denoted by $\dot{\cal N}_{\epsilon}$ , is given by the convolution

\dot{\cal N}_{\epsilon}(\epsilon)=\int_{R_{*}}^{R_{\rm top}}\int_{0}^{\infty}\frac{\dot{\cal N}^{\,\rm G}_{\epsilon}(R_{0},\epsilon_{0},\epsilon)}{\dot{N}_{0}}\ \epsilon_{0}^{2}\,\Omega R_{0}^{2}\,Q(R_{0},\epsilon_{0})\,d\epsilon_{0}\,dR_{0}\ ,

(137)

where $\dot{\cal N}^{\,\rm G}_{\epsilon}$ is evaluated using Equation (136). Since Equation (136) gives $\dot{\cal N}^{\,\rm G}_{\epsilon}$ as a function of optical depth rather than radius, it is convenient to transform the spatial variable of integration in Equation (137) from $R_{0}$ to $\tau_{0}$ using Equation (117), which yields (cf. Equation (118))

\dot{\cal N}_{\epsilon}(\epsilon)=\int_{0}^{\tau_{\rm top}}\int_{0}^{\infty}\frac{\dot{\cal N}^{\,\rm G}_{\epsilon}(\tau_{0},\epsilon_{0},\epsilon)}{\dot{N}_{0}}\ \epsilon_{0}^{2}\,Q(\tau_{0},\epsilon_{0})\,\frac{\alpha\Omega R_{*}^{4}\,y_{0}^{4}|u|}{3k_{\infty}^{2}R_{g}}\,d\epsilon_{0}\,d\tau_{0}\ ,

(138)

where $y_{0}$ and $\tau_{0}$ are related via Equation (52), and $\dot{\cal N}^{\,\rm G}_{\epsilon}$ is computed using Equation (136). The quantity $\dot{\cal N}_{\epsilon}\,d\epsilon$ gives the number of photons with energy between $\epsilon$ and $\epsilon+d\epsilon$ escaping through the top of the accretion column per unit time due to the photon source function $Q$ . Hence we can use Equation (138) to compute the number distribution of the photons escaping through the top of the conical accretion column for any desired source function $Q$ .

6 MODEL PARAMETERS AND CONSTRAINTS

If we adopt canonical values for the neutron star mass $M_{*}$ and radius $R_{*}$ , then the remaining fundamental physical free parameters in our model comprise the accretion rate, $\dot{M}$ , the electron temperature, $T_{e}$ , the magnetic field strength, $B$ , the dimensionless radius at the column top, $y_{\rm top}\equiv R_{\rm top}/R_{*}$ , the electron scattering cross sections, $\sigma_{\perp}$ , $\sigma_{{}_{||}}$ , and $\overline{\sigma}$ , the dynamical constants $k_{0}$ and $k_{\infty}$ , and the outer and inner column wall angles, $\Theta_{1}$ and $\Theta_{2}$ , respectively. The theory also includes three dimensionless auxiliary parameters that are related to the physical parameters mentioned above, namely $\xi$ (Equation (38)), $\alpha$ (Equation (44)), and $\psi$ (Equation (61)). In this section, we explore the relationships between these three dimensionless quantities and the fundamental physical free parameters. We also consider the thermodynamic structure of the accretion column and establish a framework for determining the location of the thermal mound.

6.1 Scattering Cross Sections

The detailed energy and angular dependences of the vacuum-modified electron scattering cross sections for the extraordinary and ordinary polarization modes are given by Equations (14) and (15), respectively. Since it is not possible to include the full complexity of these expressions into the model developed here, we have followed BW07 and Wang & Frank (1981) by introducing a set of approximate scattering cross sections, averaged over the photon energy and the two polarization modes. The cross sections for photons propagating parallel and perpendicular to the radial direction are denoted by $\sigma_{{}_{||}}$ and $\sigma_{\perp}$ , respectively, and the angle-averaged cross section is denoted by $\overline{\sigma}$ . In our applications to specific sources, we shall generally treat the dimensionless constants $\alpha$ , $\xi$ , and $\psi$ as free parameters, and derive the set of associated scattering cross sections $\sigma_{\perp}$ , $\sigma_{{}_{||}}$ , and $\overline{\sigma}$ using Equations (38), (44), and (61). The procedure is described below.

First, we can derive an explicit expression for $\sigma_{{}_{||}}$ by rearranging Equation (44) to obtain

\frac{\sigma_{{}_{||}}}{\sigma_{{}_{\rm T}}}=\frac{3k_{\infty}^{2}}{\alpha}\left(\frac{\Omega}{4\pi}\right)\left(\frac{\dot{M}}{\dot{M}_{\rm E}}\right)^{-1}\ ,

(139)

where the solid angle $\Omega$ is a function of the outer and inner column wall angles $\Theta_{1}$ and $\Theta_{2}$ via Equation (23). Next, we can combine Equations (38) and (139) to obtain an expression for $\sigma_{\perp}$ , which yields

\frac{\sigma_{\perp}}{\sigma_{{}_{\rm T}}}=\frac{\alpha}{3k_{\infty}^{2}b^{2}\xi^{2}}\left(\frac{\Omega}{4\pi}\right)\left(\frac{\dot{M}}{\dot{M}_{\rm E}}\right)^{-1}\left(\frac{R_{*}}{R_{g}}\right)^{2}\ ,

(140)

where the parameter $b$ is a function of $\Theta_{1}$ and $\Theta_{2}$ via Equation (28). Finally, we can obtain an explicit expression for $\overline{\sigma}$ by rearranging Equation (61) to obtain

\frac{\overline{\sigma}}{\sigma_{{}_{\rm T}}}=\frac{k_{\infty}^{2}}{\psi}\left(\frac{\Omega}{4\pi}\right)\left(\frac{\dot{M}}{\dot{M}_{\rm E}}\right)^{-1}\left(\frac{kT_{e}}{m_{e}c^{2}}\right)^{-1}\ .

(141)

For given values of the free parameters $\dot{M}$ , $T_{e}$ , $\Theta_{1}$ , $\Theta_{2}$ , $k_{\infty}$ , $\alpha$ , $\xi$ , and $\psi$ , we can compute the set of scattering cross sections $\sigma_{{}_{||}}$ , $\sigma_{\perp}$ , and $\overline{\sigma}$ using Equations (139), (140), and (141), respectively, and this is the procedure used to determine the cross sections in our model. In Section 9, we will confirm that the values obtained for the cross sections in our applications to specific pulsars are reasonably consistent with the results obtained using the detailed scattering cross sections given by Equations (14) and (15).

6.2 Thermal Mound Properties

The “thermal mound” in an X-ray pulsar accretion column is the effective surface for photon creation and destruction, which occurs mainly via free-free emission and absorption. Inside the thermal mound, we expect thermodynamic equilibrium to prevail, and above the mound, the opacity is dominated by electron scattering. We can therefore define the top of the thermal mound as the radius at which the free-free absorption optical thickness across the column is equal to unity, so that

\tau^{\rm ff}\equiv r_{\perp}\,\alpha_{\rm R}^{\rm ff}=1\ ,

(142)

where $r_{\perp}$ is the radius of the accretion column, measured perpendicular to the column axis, and the Rosseland mean absorption coefficient for the free-free process in pure, fully-ionized hydrogen gas, is given in cgs units by (Rybicki & Lightman, 1979)

\alpha_{\rm R}^{\rm ff}=A_{0}\ T_{\rm th}^{-7/2}\ \rho_{\rm th}^{2}\ ,

(143)

where

A_{0}\equiv 6.10\times 10^{22}\ ,

(144)

and $T_{\rm th}$ and $\rho_{\rm th}$ denote the temperature and density at the top of the thermal mound, respectively. Using Equation (27) to substitute for $r_{\perp}$ in Equation (142) yields

b\,R_{\rm th}\,\alpha_{\rm R}^{\rm ff}=1\ ,

(145)

where the geometrical constant $b$ is computed using Equation (28), and $R_{\rm th}$ denotes the radius at the top of the thermal mound.

The temperature of the gas in the thermal mound, $T_{\rm th}$ , is computed using the energy conservation relation (see Equation (91) from BW07)

\sigma_{{}_{\rm SB}}\,T_{\rm th}^{4}=\frac{1}{2}\,J_{\rm th}\,\upsilon_{\rm th}^{2}\ ,

(146)

where $\upsilon_{\rm th}\equiv\upsilon(R_{\rm th})$ denotes the accretion velocity at the top of the thermal mound, $\sigma_{\rm SB}$ is the Stephan-Boltzmann constant, and the mass flux at the top of the thermal mound, $J_{\rm th}$ , is given by

J_{\rm th}=\rho_{\rm th}\,|\upsilon_{\rm th}|\ .

(147)

Using Equation (147) to eliminate $\rho_{\rm th}$ in Equation (143) and combining the result with Equation (145) yields the condition

\upsilon^{2}_{\rm th}(R_{\rm th}bA_{0})^{-8/15}J_{\rm th}^{-3/5}(2\sigma_{\rm SB})^{-7/15}=1\ ,

(148)

which is satisfied at the top of the thermal mound. In the conical accretion flow treated here, with solid angle $\Omega$ , the mass flux at the top of the thermal mound, $J_{\rm th}$ , is related to the accretion rate, $\dot{M}$ , via

J_{\rm th}=\frac{\dot{M}}{\Omega R_{\rm th}^{2}}\ .

(149)

The accretion velocity at the top of the thermal mound, $\upsilon_{\rm th}$ , is related to the thermal mound radius, $R_{\rm th}$ , via Equation (45), which yields

\upsilon_{\rm th}^{2}\equiv\upsilon^{2}(y_{\rm th})=c^{2}\left(\frac{2R_{g}}{R_{*}}\right)\,\frac{k_{\infty}^{2}(y_{\rm th}^{3}-1)+k_{0}^{2}}{y_{\rm th}^{4}}\ ,

(150)

where the dimensionless thermal mound radius, $y_{\rm th}$ , is related to $R_{\rm th}$ via

y_{\rm th}\equiv\frac{R_{\rm th}}{R_{*}}\ .

(151)

Combining Equations (148) and (149) and utilizing Equation (150) to substitute for the accretion velocity, $\upsilon_{\rm th}$ , yields an equation satisfied by $y_{\rm th}$ . The result obtained is

y_{\rm th}^{-10/3}[(y_{\rm th}^{3}-1)k_{\infty}^{2}+k_{0}^{2}]\left(\frac{2R_{g}}{R_{*}}\right)c^{2}(R_{*}bA_{0})^{-8/15}\left(\frac{\dot{M}}{R_{*}^{2}\Omega}\right)^{-3/5}(2\sigma_{\rm SB})^{-7/15}=1\ .

(152)

The dimensionless thermal mound radius, $y_{\rm th}$ , is the root of this equation, if one exists for $y_{\rm th}>1$ . If no such root exists, then the thermal mound is a thin “hot spot” on the stellar surface, and therefore $y_{\rm th}=1$ . The associated thermal mound altitude, $z_{\rm th}$ , is defined by

z_{\rm th}\equiv(y_{\rm th}-1)R_{*}\ .

(153)

Once the value of $y_{\rm th}$ has been determined, we can compute the associated accretion velocity, $\upsilon_{\rm th}$ , using Equation (150). Thereafter, the thermal mound density, $\rho_{\rm th}$ , can be evaluated by combining Equations (147), (149), and (151) to obtain

\rho_{\rm th}=\frac{\dot{M}}{\Omega\,R_{*}^{2}\,y_{\rm th}^{2}\,\upsilon_{\rm th}}\ .

(154)

Next, we can compute the thermal mound temperature, $T_{\rm th}$ , by combining Equations (146) and (147), which yields

T_{\rm th}=\left(\frac{\rho_{\rm th}\,\upsilon_{\rm th}^{3}}{2\,\sigma_{\rm SB}}\right)^{1/4}\ .

(155)

The electron scattering optical depth measured from the surface of the neutron star to the top of the thermal mound is defined by

\tau_{\rm th}\equiv\tau(y_{\rm th})\ ,

(156)

which is evaluated using Equation (52). We will use $\tau_{\rm th}$ as the lower limit for the spatial integrations performed in Section 7 when we calculate the emergent spectra resulting from the Comptonization of bremsstrahlung and cyclotron seed photons.

7 PHOTON SOURCES AND ASSOCIATED SPECTRA

The X-ray spectrum observed from an accretion-powered X-ray pulsar is generated via the thermal and bulk Comptonization of seed photons injected into the accreting gas via the cyclotron, blackbody, and bremsstrahlung emission processes (BW07), which are each represented by different source functions $Q$ in the fundamental transport equation (Equation (21)). Because the transport equation is linear, the emergent spectrum can be treated as a superposition of the individual spectra resulting from the reprocessing of the various seed photon populations. For a given source function $Q$ , the general expression for the height-dependent photon number distribution escaping through the column walls, $\dot{N}_{\epsilon}$ , is given by the integral convolution in Equation (131). Likewise, the corresponding particular solution for the column-integrated photon number distribution, $\Phi_{\epsilon}$ , escaping through the walls of the accretion column, is computed using the integral convolution given by Equation (134). Finally, the particular solution for the photon number distribution, $\Psi_{\epsilon}$ , escaping through the top of the accretion column is given by the integral convolution in Equation (138).

7.1 Seed Photon Source Functions

The primary sources of seed photons in the context of accretion-powered X-ray pulsars are blackbody, cyclotron, and bremsstrahlung emission (Arons et al., 1987). Blackbody emission occurs at the top of the “thermal mound,” which is the “photosphere” for the accretion column. The top of the thermal mound is expected to be located near the stellar surface, at the base of the accretion column, where the gas achieves sufficient optical thickness to thermal absorption. Blackbody emission is a broad-band source of seed photons, but it is localized in space, since it is concentrated at the top of the thermal mound. On the other hand, cyclotron emission occurs throughout the accretion column as the result of electron-ion collisions, which can excite electrons to the $n=1$ Landau level. Radiative deexcitation back to the ground state results in the emission of cyclotron photons. Finally, bremsstrahlung emission also occurs throughout the accretion column, but in contrast to cyclotron emission, bremsstrahlung is a broad-band process that generates a roughly uniform energy distribution up to a high-energy exponential turnover at a photon energy comparable to the electron thermal energy. Photons produced via any of these three mechanisms are subsequently Comptonized before escaping from the accretion column.

The radiation source function $Q$ appearing in the transport Equation (21) is related to the photon emissivity, $\dot{n}_{\epsilon}$ , via

\epsilon^{2}Q(R,\epsilon)=\dot{n}_{\epsilon}(R,\epsilon)\ ,

(157)

where $\dot{n}_{\epsilon}(R,\epsilon)\,d\epsilon$ expresses the number of photons injected into the accretion column per unit time per unit volume at radius $R$ with energy between $\epsilon$ and $\epsilon+d\epsilon$ . In this section we will use Equation (157) to derive the source functions $Q$ for bremsstrahlung, cyclotron, and blackbody emission based on the associated emissivity functions $\dot{n}_{\epsilon}$ .

7.2 Cyclotron Radiation

Cyclotron photons are produced when electrons are excited to the $n=1$ Landau level as a result of collisions with protons. The excitations are rapidly followed by a radiative decay back to the ground state ( $n=0$ ). Estimates indicate that in X-ray pulsars, radiative deexcitation dominates over collisional deexcitation, and therefore the production of cyclotron photons tends to cool the gas (e.g., Nagel, 1980). The cyclotron source function is essentially a monochromatic function centered at the cyclotron energy $\epsilon_{c}$ (Equation (3)), although thermal effects can lead to significant broadening of the monochromatic source, and additional broadening can also occur as a result of magnetic field gradients. However, the dominant broadening mechanism is thermal Comptonization, which tends to drive any source function towards a Wien spectrum (Rybicki & Lightman, 1979).

Assuming that the accreting gas is composed of pure, fully-ionized hydrogen, we can use Equations (7) and (11) from Arons et al. (1987) or Equation (114) from BW07 to obtain for the cyclotron photon emissivity in cgs units

\dot{n}_{\epsilon}^{\rm cyc}=2.10\times 10^{36}\,\rho^{2}\,B_{12}^{-3/2}\,H\left(\frac{\epsilon_{c}}{kT_{e}}\right)\,e^{-\epsilon_{c}/kT_{e}}\,\delta(\epsilon-\epsilon_{c})\ ,

(158)

where $\rho=m_{p}n_{e}$ denotes the mass density, the cyclotron energy $\epsilon_{c}$ is given by Equation (3), and the function $H$ is defined by

H\left(\frac{\epsilon_{c}}{kT_{e}}\right)\equiv\begin{cases}0.41\ ,&\epsilon_{c}/kT_{e}>7.5\ ,\\ 0.15\,\sqrt{\epsilon_{c}/kT_{e}}\ ,&\epsilon_{c}/kT_{e}<7.5\ .\\ \end{cases}

(159)

The associated source function, $Q^{\rm cyc}$ , is given by the expression

Q^{\rm cyc}(R,\epsilon)=2.10\times 10^{36}\,\rho^{2}\,B_{12}^{-7/2}\,H\left(\frac{\epsilon_{c}}{kT_{e}}\right)\,\epsilon^{-2}\,e^{-\epsilon_{c}/kT_{e}}\,\delta(\epsilon-\epsilon_{c})\ ,

(160)

obtained by combining Equations (157) and (158). The cyclotron source function $Q^{\rm cyc}(R,\epsilon)$ is operative between the thermal mound radius, $R_{\rm th}$ , and the top of the accretion column, $R_{\rm top}$ , and it is concentrated near the stellar surface due to the appearance of the factor $\rho^{2}$ , which reflects the two-body nature of the excitation process. Because of the concentration of the emission process near the stellar surface, we can assume a roughly constant value for the magnetic field strength in the cyclotron emission region.

The particular solution for the column-integrated photon number spectrum resulting from the escape of reprocessed cyclotron seed photons through the walls of the accretion column, denoted by $\Phi^{\rm cyc}_{\epsilon}(\epsilon)$ , is obtained by using Equation (160) to substitute for $Q$ in Equation (134). Due to the delta-function dependence on the photon energy $\epsilon$ , the energy integration is trivial. By integrating over $\tau_{0}$ term-by-term, after some algebra we find that the particular solution for the column-integrated escaping photon number spectrum due to Comptonized cyclotron emission is given in cgs units by

	$\displaystyle\Phi^{\rm cyc}_{\epsilon}(\epsilon)$	$\displaystyle=$	$\displaystyle\frac{1.15\times 10^{-12}\dot{M}^{2}\,T_{e}\,\psi\,\xi^{2}\,H(\chi_{c})\,\epsilon^{\kappa-2}}{k_{\infty}^{2}\sqrt{R_{g}R_{*}}\,\Omega\,e^{(\epsilon+\epsilon_{c})/(2kT_{e})}\,\epsilon_{c}^{\kappa+3/2}}\sum_{n=0}^{\infty}\frac{\Gamma(\mu-\kappa+1/2)}{\Gamma(1+2\mu)}\frac{I_{n}\,J_{n}}{{\cal I}_{n}}$
		$\displaystyle\times$	$\displaystyle M_{\kappa,\mu}[\min(\chi,\chi_{c})]\,W_{\kappa,\mu}[\max(\chi,\chi_{c})]\ ,\phantom{\left(\frac{\epsilon_{\rm min}}{kT_{e}}\right)}$		(161)

where

\chi\equiv\frac{\epsilon}{kT_{e}}\ ,\qquad\chi_{c}\equiv\frac{\epsilon_{c}}{kT_{e}}\ ,

(162)

and the quantity $J_{n}$ represents the spatial integral

J_{n}\equiv\int_{\tau_{\rm th}}^{\tau_{\rm top}}\frac{\omega(\tau_{0})\,g_{n}(\tau_{0})}{\sqrt{(y_{0}^{3}-1)k_{\infty}^{2}+k_{0}^{2}}}\,d\tau_{0}\ .

(163)

The weight function $\omega(\tau)$ appearing in Equation (163) is defined in Equation (95), and the quantities $\tau_{\rm top}$ (Equation (71)) and $\tau_{\rm th}$ (Equation (156)) represent the scattering optical depths at the top of the accretion column and at the surface of the thermal mound, respectively. The integrals ${\cal I}_{n}$ and $I_{n}$ in Equation (161) are defined in Equations (112) and (127), respectively. The quantity $\Phi^{\rm cyc}_{\epsilon}(\epsilon)$ gives the number of photons escaping through the walls of the accretion column per unit time in the energy range between $\epsilon$ and $\epsilon+d\epsilon$ .

We can also obtain the particular solution for the photon number spectrum resulting from the escape of Comptonized cyclotron emission through the top of the accretion column, denoted by $\dot{\cal N}_{\epsilon}^{\rm cyc}(\epsilon)$ , by using Equation (160) to substitute for $Q$ in Equation (138). After simplification, the result obtained in cgs units is

	$\displaystyle\dot{\cal N}_{\epsilon}^{\rm cyc}(\epsilon)$	$\displaystyle=$	$\displaystyle\frac{5.76\times 10^{-13}\dot{M}^{2}\,T_{e}\,y_{\rm top}^{2}\,\psi\sqrt{R_{g}R_{*}}\,H(\chi_{c})\,\epsilon^{\kappa-2}}{\Omega R_{g}^{2}\,k_{\infty}^{2}\,e^{(\epsilon+\epsilon_{c})/(2kT_{e})}\,\epsilon_{c}^{\kappa+3/2}}\sum_{n=0}^{\infty}\frac{\Gamma(\mu-\kappa+1/2)}{\Gamma(1+2\mu)}\frac{J_{n}\,g_{n}(\tau_{\rm top})}{{\cal I}_{n}}$
		$\displaystyle\times$	$\displaystyle M_{\kappa,\mu}[\min(\chi,\chi_{c})]\,W_{\kappa,\mu}[\max(\chi,\chi_{c})]\ ,\phantom{\left(\frac{\epsilon_{\rm min}}{kT_{e}}\right)}$		(164)

where $\chi$ and $\chi_{c}$ are defined in Equation (162). The quantity $\dot{\cal N}_{\epsilon}^{\rm cyc}(\epsilon)$ gives the number of photons escaping through the top of the accretion column per unit time in the energy range between $\epsilon$ and $\epsilon+d\epsilon$ .

7.3 Blackbody Radiation

At the base of the accretion column, the gas becomes sufficiently dense that thermodynamic equilibrium is established, in the region referred to as the thermal mound. The upper surface of the thermal mound radiates a blackbody distribution, with surface energy flux equal to $\pi$ times the Planck intensity function, $B_{\epsilon}(\epsilon)$ , given by (Rybicki & Lightman, 1979)

B_{\epsilon}(\epsilon)=\frac{2\,\epsilon^{3}}{c^{2}h^{3}}\,\frac{1}{e^{\epsilon/kT_{\rm th}}-1}\ ,

(165)

where $T_{\rm th}$ is the temperature of the gas in the thermal mound. Note that the function $B_{\epsilon}$ denotes the blackbody “energy intensity,” with units ${\rm ergs\ s^{-1}\,ster^{-1}\,cm^{-2}\,erg^{-1}}$ . Following Becker & Wolff (2005b), and BW07, we define the thermal mound surface emission function, $S(\epsilon)$ , using

\epsilon^{3}S(\epsilon)\,d\epsilon\equiv\pi\,B_{\epsilon}(\epsilon)\,d\epsilon\ .

(166)

According to Equation (166), the quantity $\epsilon^{3}S(\epsilon)\,d\epsilon$ gives the energy emitted per unit time per unit area from the upper surface of the thermal mound in the energy range between $\epsilon$ and $\epsilon+d\epsilon$ .

The function $S$ is related to the source term $Q$ appearing in the transport Equation (21) via

Q^{\rm bb}(R,\epsilon)\equiv S(\epsilon)\,\delta(R-R_{\rm th})\ ,

(167)

which can be combined with Equations (165) and (166) to obtain

Q^{\rm bb}(R,\epsilon)=\frac{2\pi}{c^{2}h^{3}}\,\frac{\delta(R-R_{\rm th})}{e^{\epsilon/kT_{\rm th}}-1}\ .

(168)

This result for the blackbody source term is localized in physical space, but it possesses a distributed (broad-band) energy dependence, which is the exact opposite of the cyclotron source given by Equation (160).

We can obtain the particular solution for the column-integrated photon number spectrum resulting from the escape of Comptonized blackbody radiation through the walls of the accretion column by using Equation (168) to substitute for $Q$ in Equation (133). The resulting spatial integration is trivial, and only the energy integral remains, which reduces to

\Phi^{\rm bb}_{\epsilon}(\epsilon)=\frac{2\pi\,\Omega R_{\rm th}^{2}}{c^{2}h^{3}}\int_{0}^{\infty}\frac{\Phi_{\epsilon}^{\rm G}(R_{\rm th},\epsilon_{0},\epsilon)}{\dot{N}_{0}}\ \frac{\epsilon_{0}^{2}}{e^{\epsilon_{0}/kT_{\rm th}}-1}\ d\epsilon_{0}\ ,

(169)

where $\Phi_{\epsilon}^{\rm G}$ is computed using Equation (126). Likewise, we can obtain the particular solution for the photon number spectrum resulting from the escape of Comptonized blackbody photons through the top of the accretion column by combining Equations (137) and (168), which yields

\dot{\cal N}_{\epsilon}^{\rm bb}(\epsilon)=\frac{2\pi\,\Omega R_{\rm th}^{2}}{c^{2}h^{3}}\int_{0}^{\infty}\frac{\dot{\cal N}^{\,\rm G}_{\epsilon}(R_{\rm th},\epsilon_{0},\epsilon)}{\dot{N}_{0}}\ \,\frac{\epsilon_{0}^{2}}{e^{\epsilon_{0}/kT_{\rm th}}-1}\ d\epsilon_{0}\ ,

(170)

where $\dot{\cal N}^{\,\rm G}_{\epsilon}$ is evaluated using Equation (136).

7.4 Bremsstrahlung Radiation

Bremsstrahlung is the dominant contributor to the seed photon distribution in luminous accretion-powered X-ray pulsars such as Her X-1, LMC X-4, and Cen X-3 (e.g., Becker & Wolff, 2007; Wolff et al., 2016; West et al., 2017a, b). Following BW07, we can use Equation (5.14b) from Rybicki & Lightman (1979) to write the bremsstrahlung (free-free) photon production rate per unit volume per unit energy in cgs units as

\dot{n}_{\epsilon}^{\rm ff}=3.68\times 10^{36}\,\rho^{2}\,T_{e}^{-1/2}\,\epsilon^{-1}\,e^{-\epsilon/kT_{e}}\ ,

(171)

for a plasma composed of pure, fully-ionized hydrogen. Equation (171) can be combined with Equation (157) to obtain the free-free source function

Q^{\rm ff}(R,\epsilon)=3.68\times 10^{36}\,\rho^{2}\,T_{e}^{-1/2}\,\epsilon^{-3}\,e^{-\epsilon/kT_{e}}\ ,\ \ \ \ \ R_{\rm th}<R<R_{\rm top}\ ,

(172)

where $R_{\rm th}$ and $R_{\rm top}$ denote the radii at the surface of the thermal mound and at the top of the accretion column, respectively. Equation (172) gives the spectrum of the injected, optically-thin bremsstrahlung radiation for photons with energy $\epsilon>\epsilon_{\rm abs}$ , where $\epsilon_{\rm abs}$ denotes the self-absorption cutoff, below which the plasma is optically thick and therefore the radiation is thermalized. The value of $\epsilon_{\rm abs}$ is related to the gas density $\rho$ and the electron temperature $T_{e}$ via Equation (127) from BW07, which gives

\frac{\epsilon_{\rm abs}}{kT_{e}}=6.08\times 10^{12}\,T_{e}^{-7/4}\,\rho^{1/2}\ .

(173)

In our applications, we set $\epsilon_{\rm abs}/kT_{e}=0.05$ , which is a reasonable approximation for the calculations performed here, since the emergent spectrum has a weak dependence on this parameter.

The bremsstrahlung source term given by Equation (172) has non-trivial dependences on both the radius $R$ and the photon energy $\epsilon$ and hence it is more complex to implement computationally than either the cyclotron or blackbody sources given by Equations (160) and (168), respectively. Using Equation (172) to substitute for $Q$ in Equation (133), we find that, after performing the spatial and energy integrals and simplifying, the result obtained for the column-integrated photon number spectrum resulting from the escape of Comptonized bremsstrahlung through the walls of the accretion column can be written in cgs units as

\Phi^{\rm ff}_{\epsilon}(\epsilon)=\frac{6.80\times 10^{7}\,\dot{M}^{2}\,\psi\,\xi^{2}\,\epsilon^{\kappa-2}\,e^{-\epsilon/(2kT_{e})}}{k_{\infty}^{2}\sqrt{R_{g}R_{*}}\,\Omega\,(kT_{e})^{\kappa-1/2}}\ \sum_{n=0}^{\infty}\ \frac{\Gamma(\mu-\kappa+1/2)}{\Gamma(1+2\mu)}\ \frac{I_{n}\,J_{n}}{{\cal I}_{n}}\,B_{n}\ ,

(174)

where ${\cal I}_{n}$ , $I_{n}$ , and $J_{n}$ are computed using Equations (112), (128), and (163), respectively. The quantity $B_{n}$ introduced in Equation (174) denotes the energy integral

B_{n}\equiv\int_{\chi_{\rm abs}}^{\infty}\chi_{0}^{-1-\kappa}\,e^{-\chi_{0}/2}\,M_{\kappa,\mu}[\min(\chi,\chi_{0})]\,W_{\kappa,\mu}[\max(\chi,\chi_{0})]\,d\chi_{0}\ ,

(175)

where $\chi\equiv\epsilon/kT_{e}$ , $\chi_{0}\equiv\epsilon_{0}/kT_{e}$ , and $\chi_{\rm abs}\equiv\epsilon_{\rm abs}/kT_{e}$ .

We can also derive the particular solution for the photon number spectrum escaping through the top of the accretion column as a result of Comptonized bremsstrahlung emission by combining Equations (138) and (172). After some algebra, the result obtained in cgs units is

	$\displaystyle\dot{\cal N}_{\epsilon}^{\rm ff}(\epsilon)$	$\displaystyle=$	$\displaystyle\frac{3.40\times 10^{7}\,\dot{M}^{2}\,y_{\rm top}^{2}\,\psi\,\epsilon^{\kappa-2}\,e^{-\epsilon/(2kT_{e})}\sqrt{R_{g}R_{*}}}{k_{\infty}^{2}R_{g}^{2}\,\Omega\,(kT_{e})^{\kappa-1/2}}$		(176)
		$\displaystyle\times$	$\displaystyle\ \sum_{n=0}^{\infty}\ \frac{\Gamma(\mu-\kappa+1/2)}{\Gamma(1+2\mu)}\ \frac{g_{n}(\tau_{\rm top})\,J_{n}}{{\cal I}_{n}}\,B_{n}\ .$		(176)

8 ASTROPHYSICAL APPLICATIONS

The results from the previous sections in the paper can be combined to compute the spectrum emitted by an X-ray pulsar accretion column due to the thermal and bulk Comptonization of cyclotron, bremsstrahlung, and blackbody seed photons. Because the theoretical model developed here computes separate emission components for photons escaping through either the walls or the top of the accretion column, we are able to accomplish more detailed calculations of the approximate phase-averaged spectrum than were presented by BW07. In this section we develop the formalism for calculating the contributions of the wall and top emission components to the observed X-ray spectrum, and we apply the new model to Her X-1 and X Per.

The theoretical column-integrated photon count rate spectrum observed at Earth due to the escape of Comptonized photons through the column walls, denoted by $F^{\rm wall}_{\epsilon}(\epsilon)$ , is computed using the expression

F^{\rm wall}_{\epsilon}(\epsilon)\equiv\frac{\Phi_{\epsilon}^{\rm tot}(\epsilon)}{4\pi D^{2}}\ ,

(177)

where

\Phi_{\epsilon}^{\rm tot}(\epsilon)\equiv\left[\Phi^{\rm cyc}_{\epsilon}(\epsilon)+\Phi^{\rm bb}_{\epsilon}(\epsilon)+\Phi^{\rm ff}_{\epsilon}(\epsilon)\right]A_{c}(\epsilon)

(178)

denotes the total photon number spectrum escaping through the column walls, $D$ is the distance to the source, and the cyclotron absorption feature commonly seen in X-ray pulsar spectra is represented by the function $A_{c}(\epsilon)$ , defined using the standard Gaussian form (e.g., Heindl & Chakrabarty, 1999; Orlandini et al., 1998; Soong et al., 1990)

A_{c}(\epsilon)\equiv 1-\frac{d_{c}}{\sigma_{c}\sqrt{2\pi}}\,e^{-(\epsilon-\epsilon_{c})^{2}/(2\sigma^{2}_{c})}\ .

(179)

The quantities on the right-hand side of Equation (178) express the contributions to the escaping column-wall spectrum due to Comptonized cyclotron emission ( $\Phi^{\rm cyc}_{\epsilon}$ ; Equation (161)), blackbody emission ( $\Phi^{\rm bb}_{\epsilon}$ ; Equation (169)), and bremsstrahlung emission ( $\Phi^{\rm ff}_{\epsilon}$ ; Equation (174)), respectively.

We can also use the theoretical formalism developed here to compute the observed spectrum due to the escape of Comptonized photons through the top of the accretion column. The total column-top spectrum observed at Earth, denoted by $F^{\rm top}_{\epsilon}(\epsilon)$ , is calculated using the expression

F^{\rm top}_{\epsilon}(\epsilon)\equiv\frac{\dot{\cal N}_{\epsilon}(\epsilon)}{4\pi D^{2}}\ ,

(180)

where

\dot{\cal N}^{\rm top}_{\epsilon}(\epsilon)\equiv\left[\dot{\cal N}^{\rm cyc}_{\epsilon}(\epsilon)+\dot{\cal N}^{\rm bb}_{\epsilon}(\epsilon)+\dot{\cal N}^{\rm ff}_{\epsilon}(\epsilon)\right]A_{c}(\epsilon)\ ,

(181)

where $A_{c}(\epsilon)$ is defined in Equation (179). The quantities on the right-hand side of Equation (181) denote the contributions to the column-top spectrum due to Comptonized cyclotron emission ( $\dot{\cal N}^{\rm cyc}_{\epsilon}$ ; Equation (164)), blackbody emission ( $\dot{\cal N}^{\rm bb}_{\epsilon}$ ; Equation (170)), and bremsstrahlung emission ( $\dot{\cal N}^{\rm ff}_{\epsilon}$ ; Equation (176)), respectively.

Following West et al. (2017b), we can use our model to obtain an approximation of the observed phase-averaged X-ray spectrum by adding the contributions emitted through the walls of the accretion column (Equation (177)) and through the column top (Equation (180)). The approximate phase-averaged spectrum, denoted by $F^{\rm total}_{\epsilon}(\epsilon)$ , is therefore given by the sum

F^{\rm total}_{\epsilon}(\epsilon)\equiv F^{\rm wall}_{\epsilon}(\epsilon)+F^{\rm top}_{\epsilon}(\epsilon)\ .

(182)

We will use this approach to compute theoretical predictions for the observed phase-averaged X-ray spectra of Her X-1 and X Per. These two sources were selected because they span a very wide range in X-ray luminosity, and because good quality X-ray spectra are available for both sources in the published literature. The spectral results presented here are considered example calculations, rather than detailed quantitative fits, which will require further software development and are beyond the scope of this paper. In each case we adopt canonical values for the stellar mass and radius, with $M_{*}=1.4\,M_{\odot}$ and $R_{*}=10\,$ km, respectively. The accretion rate, $\dot{M}$ , is set using observational estimates of the X-ray luminosity, combined with Equation (78), and the distance to the source, $D$ , is adopted from published results. The remaining fundamental theoretical free parameters that must be set in order to evaluate the X-ray spectrum using our model are $\alpha$ , $\xi$ , $\psi$ , $B$ , $\Theta_{1}$ , $\Theta_{2}$ , $k_{0}$ , $k_{\infty}$ , $y_{\rm top}$ , and $T_{e}$ , which are listed in Table 1 for each of the models considered here. These fundamental free parameters are varied in order to achieve a satisfactory qualitative fit to the observed phase-averaged X-ray spectrum for a given source. In addition, the free parameters are also constrained via analysis of the thermodynamic and hydrodynamic structure of the accretion column, as discussed in Section 9.

Once the fundamental free parameters $\alpha$ , $\xi$ , $\psi$ , $B$ , $\Theta_{1}$ , $\Theta_{2}$ , $k_{0}$ , $k_{\infty}$ , $y_{\rm top}$ , and $T_{e}$ are determined via the qualitative fitting process, a number of additional auxiliary parameters are computed, with values reported in Table 2. These include the perpendicular electron scattering cross section, $\sigma_{\perp}$ , the parallel electron scattering cross section, $\sigma_{{}_{||}}$ , the angle-averaged scattering cross section, $\overline{\sigma}$ , the thermal mound temperature, $T_{\rm th}$ , the thermal mound altitude, $z_{\rm th}$ , the thermal mound scattering optical depth, $\tau_{\rm th}$ , and the optical depth at the top of the column, $\tau_{\rm top}$ . We also include the radius of the magnetic cap at the base of the accretion column, defined by

r_{0}\equiv R_{*}\,\Theta_{1}\ .

(183)

For the two sources treated here, we find that the best results are obtained by setting the inner column wall angle, $\Theta_{2}$ , equal to zero, so that the accretion columns are both completely filled cones, with no hollow interior region. In addition to the requirement of an acceptable qualitative fit to the observed X-ray spectrum, we also impose additional constraints related to the dynamical and thermal structure of the accretion column (see Section 9).

The accretion dynamics in Her X-1 is expected to be dominated by radiation pressure (BW07), in which case the accretion column contains a smooth, radiative, radiation-dominated shock. In this source, most of the kinetic energy of the accreting gas is radiated away through the walls of the accretion column, and the gas settles onto the stellar surface with approximately zero velocity. Conversely, in the case of X Per, it is thought that a strong, discontinuous, gas-mediated standing shock is located at the top of the accretion column (Langer & Rappaport, 1982). As the gas passes through the shock, its velocity drops by a factor of 4 below the incident free-fall value, and in the downstream region below the shock, the gas may reaccelerate until it collides with the surface of the star. The gas may approach the stellar surface at a large fraction of the speed of light, with the final deceleration occurring via Coulomb collisions as the gas merges with the stellar crust.

Once all of the theoretical parameters are specified for a given source, the resulting X-ray count-rate spectrum is computed using Equations (177) and (180), which give the column-wall and the column-top spectra, respectively. The spectral results presented here were computed using the first 5-10 terms in the analytical expansions, which we find yields at least three decimal digits of accuracy in the resulting X-ray spectra. Photon conservation is a necessary condition for the steady-state model considered here, and therefore the integrity of the calculation is checked by confirming that the total number of photons escaping through the column walls and column top per unit time is equal to the number of seed photons injected per unit time. This check was carried out independently for each photon source mechanism considered here (bremsstrahlung, blackbody, and cyclotron).

8.1 Her X-1

In Figure 7 we plot the theoretical count-rate spectrum for Her X-1, computed using Equations (177), (180), and (182), and compare it with the phase-averaged NuSTAR spectrum reported by Wolff et al. (2016), which is based on an XSPEC analysis of an observation by Fürst et al. (2013). The deconvolved data from the NuSTAR FPMA and FPMB modules are indicated by the magenta and cyan lines, respectively. Results are presented for the total theoretical phase-averaged spectrum, as well as for the individual contributions to the observed flux due to the Comptonization of cyclotron, blackbody, and bremsstrahlung seed photons. The spectrum for each mechanism is further separated into column-wall and column-top emission components.

The input values for the fundamental theory parameters used to compute the spectral results for Her X-1 correspond to model 1, with $\dot{M}=1.90\times 10^{17}\,{\rm g\,s^{-1}}$ , $T_{e}=5.50\times 10^{7}\,$ K, $B=3.26\times 10^{12}\,$ G, $D=6.6\,$ kpc, $\alpha=0.35$ , $\xi=1.14$ , $\psi=1.42$ , $k_{0}=0$ , $k_{\infty}=1$ , $\Theta_{1}=0.315^{\circ}$ , $\Theta_{2}=0$ , and $y_{\rm top}=2.15$ , as indicated in Table 1. The associated values for the various auxiliary parameters are listed in Table 2. In particular, the values obtained for the scattering cross sections are $\sigma_{\perp}/\sigma_{{}_{\rm T}}=0.537$ , $\sigma_{{}_{||}}/\sigma_{{}_{\rm T}}=6.68\times 10^{-5}$ , and $\overline{\sigma}/\sigma_{{}_{\rm T}}=5.90\times 10^{-4}$ . Table 2 also includes values for the scattering optical depth from the stellar surface to the top of the thermal mound, $\tau_{\rm th}=0.08$ , the scattering optical depth from the surface to the top of the column, $\tau_{\rm top}=2.91$ , the distance from the stellar surface to the top of the thermal mound, $z_{\rm th}=692\,$ cm, the radius of the magnetic polar cap at the base of the flow, $r_{0}=55\,$ m, the impact velocity at the stellar surface, $\upsilon_{*}=0$ , the accretion velocity at the top of the thermal mound, $\upsilon_{\rm th}=-0.03\,$ c, and the derived temperature of the thermal mound, $T_{\rm th}=6.07\times 10^{7}\,$ K, which is the temperature used to compute the blackbody seed photon distribution. The velocity profile for the Her X-1 accretion column is computed by substituting the parameter values $k_{0}=0$ and $k_{\infty}=1$ into Equation (45), and the resulting velocity profile is plotted in Figure 8.

It is clear from the spectral components for the Her X-1 theoretical model plotted in Figure 7 that the phase-averaged X-ray spectrum for this source is dominated by Comptonized bremsstrahlung emission, radiated primarily through the column walls rather than through the top of the column. However, the cyclotron emission components make a significant contribution to the observed spectrum in the region around the cyclotron absorption feature, at photon energy $\epsilon_{c}\sim 38\,$ keV. The power-law shape of the continuum, combined with the lack of a strong Wien peak in the spectrum, indicates that thermal Comptonization is unsaturated in this source. We note that the blackbody emission components make a negligible contribution to the observed X-ray spectrum for Her X-1.

The values obtained for the electron scattering cross sections in the application to Her X-1, listed in Table 2, satisfy the relation $\sigma_{\perp}\sim\sigma_{{}_{\rm T}}\gg\overline{\sigma}\gg\sigma_{{}_{||}}$ , which indicates that the cross sections are not significantly influenced by the effects of vacuum polarization, as discussed in Section 2. This is reasonable, since the electron number density near the base of the accretion column in Her X-1 is $n_{e}\sim 10^{24}\,{\rm cm}^{-3}$ , which yields for the vacuum energy $\epsilon_{\rm vac}\sim 100\,$ keV (see Equation (11)). Recalling that vacuum polarization effects are unimportant for photon energies $\epsilon\lesssim\epsilon_{\rm vac}$ , it is clear that the X-ray continuum in Her X-1 is unaffected by this process. In addition to the effect of vacuum polarization, one must also consider the distinction between the scattering cross sections experienced by ordinary and extraordinary mode photons propagating in the strong magnetic field. According to Figure 2, the extraordinary mode photons will experience a resonance in the scattering cross section at photon energy $\epsilon\sim\epsilon_{c}$ . Hence, if extraordinary mode photons are dominant in Her X-1, then we may expect to see super-Thomson cross section values for photons energies close to $\epsilon_{c}$ . However, the dominance of bremsstrahlung emission (which is non-resonant) in the spectrum of Her X-1 implies that ordinary mode photons dominate the seed photon population in this source. Furthermore, detailed consideration of the cross sections for mode-change scattering indicates that the ordinary mode photons produced via bremsstrahlung will primarily experience “mode-coherent” scattering, and will therefore remain as ordinary mode photons. This argument supports the cross section relation $\sigma_{\perp}\sim\sigma_{{}_{\rm T}}\gg\overline{\sigma}\gg\sigma_{{}_{||}}$ which we obtain in our application of the model to Her X-1. This issue will be discussed in more detail in Section 9.

The relatively low value obtained for the thermal mound velocity, $\upsilon_{\rm th}=-0.03\,$ c, in Her X-1 indicates that radiation pressure has nearly removed all of the kinetic energy from the gas by the time it reaches the top of the thermal mound. All of the remaining kinetic energy is used to power the blackbody “hotspot” at the base of the column, which generates the blackbody seed photons, as the gas continues to decelerate due to radiation pressure. By the time the gas reaches the stellar surface, the impact velocity is $\upsilon_{*}=0$ , as implied by the value $k_{0}=0$ listed in Table 1 (see Equations (45) and (46)). The bremsstrahlung and cyclotron seed photons are generated throughout the entire accretion column, but their production is concentrated near the base of the flow since the density reaches its maximum value there. While bremsstrahlung and cyclotron emission both make important contributions to the observed X-ray spectrum in Her X-1, the blackbody component makes a negligible contribution because very little kinetic energy is left by the time the gas enters the thermal mound.

In the case of Her X-1, the spectrum plotted in Figure 7 also includes a Gaussian Fe K $\alpha$ emission line profile computed using the function

F^{\rm K}_{\epsilon}(\epsilon)\equiv\frac{d_{\rm K}}{\sigma_{\rm K}\sqrt{2\pi}}\,e^{-(\epsilon-\epsilon_{\rm K})^{2}/(2\sigma^{2}_{\rm K})}\ ,

(184)

where $d_{\rm K}$ is the total iron line photon flux measured at Earth. It is important to note that in our approach, the iron line is simply added to the final spectrum rather than being subject to Comptonization inside the accretion column. A more sophisticated approach is beyond the scope of this paper due to the complexity of the spectral formation process for the iron line photons. The associated auxiliary parameter value we used to model the iron line for Her X-1 are $\epsilon_{\rm K}=6.55\,$ keV, $d_{\rm K}=4.00\times 10^{-3}\,{\rm s^{-1}\,cm^{-2}}$ , and $\sigma_{\rm K}=0.30\,$ keV, and the parameter values used to treat cyclotron absorption via Equation (179) are are $\epsilon_{c}=37.7\,$ keV, $d_{c}=18.0\,$ keV, and $\sigma_{c}=10.0\,$ keV. We note that most of the parameter values obtained here are very similar to those reported by BW07 and by Wolff et al. (2016) in their simulations of the phase-averaged X-ray spectrum of Her X-1.

8.2 X Per

In Figure 9 we compare the theoretical count-rate spectrum computed using Equations (177), (180), and (182) with the deconvolved, phase-averaged RXTE spectrum for X Per, which is based on an XSPEC analysis of an archival observation taken in 1998 July and reported by Delgado-Martí et al. (2001). The RXTE spectrum considered here is the same one analyzed by Becker & Wolff (2005a, b) using their bulk Comptonization model. The theoretical results plotted in Figure 9 include the total phase-averaged spectrum, and the individual column-wall and column-top contributions to the observed spectrum due to the Comptonization of cyclotron, blackbody, and bremsstrahlung seed photons.

The magnetic field strength in X Per has been the subject of intense debate for several decades, and it has been estimated using a variety of techniques. For example, Ghosh & Lamb (1979) used their model for torque-driven pulsar spin-down to estimate the surface magnetic field strength in X Per as $B=4.8\times 10^{12}\,$ G, assuming a canonical neutron star with radius $R_{*}=10\,$ km. Di Salvo et al. (1998) used an analysis of BeppoSAX observations of X Per to obtain the similar estimate $B=5.6\times 10^{12}\,$ G, and Coburn et al. (2001) used the detection of a putative cyclotron absorption feature in the RXTE spectrum of X Per to conclude that the surface magnetic field strength is $B=3.3\times 10^{12}\,$ G. Maitra et al. (2017) used an analysis of Suzaku data to obtain the similar value $B=3.4\times 10^{12}\,$ G. On the other hand, Yatabe et al. (2018) have argued that application of the Ghosh & Lamb (1979) theory to the interpretation of MAXI data yields a best fit that indicates a much higher estimate for the magnetic field strength in the range $B\sim 10^{13-14}\,$ G. Given the uncertainty in the data and the lack of consensus regarding the strength of the magnetic field in this source, we have chosen to adopt the representative value $B=3.0\times 10^{12}\,$ G.

The theoretical spectrum for X Per plotted in Figure 9 is based on model 2, with fundamental theory parameters given by $\dot{M}=4.16\times 10^{14}\,{\rm g\,s^{-1}}$ , $T_{e}=1.20\times 10^{8}\,$ K, $B=3.00\times 10^{12}\,$ G, $D=0.725\,$ kpc, $\alpha=0.08$ , $\xi=1.20$ , $\psi=0.75$ , $k_{0}=0.3$ , $k_{\infty}=0.244$ , $\Theta_{1}=5.536^{\circ}$ , $\Theta_{2}=0$ , and $y_{\rm top}=2.20$ , as listed in Table 1. The associated auxiliary parameters are listed in Table 2. In particular, the values obtained for the scattering cross sections are $\sigma_{\perp}/\sigma_{{}_{\rm T}}=848$ , $\sigma_{{}_{||}}/\sigma_{{}_{\rm T}}=2.46$ , and $\overline{\sigma}/\sigma_{{}_{\rm T}}=4.31$ . These cross section values are radically different from those obtained in the case of Her X-1, which reflects the distinctly different nature of the conditions in the accreting plasma in the two sources. This is further discussed below. Table 2 also includes the scattering optical depth at the top of the thermal mound, $\tau_{\rm th}=0$ , the scattering optical depth at the top of the accretion column, $\tau_{\rm top}=1.75$ , the distance between the stellar surface and the top of the thermal mound, $z_{\rm th}=0$ , the magnetic polar cap radius, $r_{0}=966\,$ m, the impact velocity at the stellar surface, $\upsilon_{*}=-0.19\,c$ , the accretion velocity at the top of the thermal mound, $\upsilon_{\rm th}=-0.19\,$ c, and the thermal mound temperature, $T_{\rm th}=8.04\times 10^{6}\,$ K, which is used to compute the blackbody seed photon distribution.

We note that in the case of X Per, the thermal mound is located directly on the stellar surface, because there is insufficient absorption optical depth to form an extended region of thermal equilibrium. Hence the thermal mound velocity, $\upsilon_{\rm th}$ , is equal to the impact velocity onto the stellar surface, $\upsilon_{*}$ . The theoretical spectrum plotted in Figure 9 also includes interstellar absorption, with a hydrogen column density $N_{\rm H}=3.0\times 10^{21}\,{\rm cm}^{-2}$ (Becker & Wolff, 2005a, b; La Palombara & Mereghetti, 2007). However, the model for X Per does not include either the cyclotron absorption feature (Equation (179)) or the iron emission line (Equation (184)), since there is no strong observational evidence for either of these features in the X-ray spectrum for this source. The associated velocity profile for the X Per accretion column is obtained by substituting the parameter values $k_{0}=0.3$ and $k_{\infty}=0.244$ into Equation (45), and is plotted in Figure 10.

Comparison of the various spectral components plotted in Figure 9 with the RXTE data for X Per clearly indicates that the X-ray spectrum is dominated by Comptonized blackbody radiation, escaping through the top of the accretion column. The seed blackbody photons for this component are generated in the “hot spot” at the base of the accretion column, where the gas collides at high speed with the surface of the star. According to Equations (45) and (46), the surface velocity is determined by the value of the parameter $k_{0}$ , with $k_{0}=0.3$ in the case of X Per. The resulting large impact velocity, $\upsilon_{*}=-0.19\,c$ , reflects the inability of radiation pressure to decelerate the gas once it passes through the standing, discontinuous shock at the top of the accretion column. Instead, the gas gradually reaccelerates in the post-shock region (see Figure 10). The kinetic energy associated with the large impact velocity powers the blackbody seed emission, making it a much more effective source of seed photons for the Comptonization process than either bremsstrahlung or cyclotron. These behaviors are qualitatively different from those observed in the case of Her X-1, and stem from the fact that the luminosity of X Per is roughly three orders of magnitude lower, which profoundly alters the accretion dynamics and the properties of the accreting plasma.

The values obtained for the electron scattering cross sections in the application of the model to X Per are listed in Table 2. In contrast to the results obtained in the case of Her X-1, we find that in the case of X Per, the cross sections are very strongly influenced by the effects of vacuum polarization (see Section 2). This is expected, because the electron number density near the base of the accretion column in X Per is $n_{e}\sim 10^{18}\,{\rm cm}^{-3}$ , which yields for the vacuum energy $\epsilon_{\rm vac}\sim 0.1\,$ keV (see Equation (11)). Since vacuum polarization is important for photons with energy $\epsilon\gtrsim\epsilon_{\rm vac}$ , we conclude that the entire X-ray continuum will be strongly modified by this effect in X Per. Further insight into this issue can be obtained by examining the distinction between the scattering cross sections for the ordinary and extraordinary mode photons. According to Figures 3 and 4, both the ordinary and extraordinary mode photons with energy $\epsilon\sim\epsilon_{c}$ in X Per will experience super-Thomson scattering cross sections due to vacuum polarization. However, the angular dependence of the cross sections for the two modes is quite different. In particular, the cross section for the extraordinary mode photons increases as a function of the propagation angle $\theta$ , whereas it decreases with increasing $\theta$ for the ordinary mode. Keeping in mind that photons propagating perpendicular to the magnetic field see the perpendicular scattering cross section, $\sigma_{\perp}$ , we conclude that if extraordinary mode photons dominate in X Per, then we would expect to see super-Thomson cross sections, with $\sigma_{\perp}\gg\sigma_{{}_{||}}\gtrsim\sigma_{{}_{\rm T}}$ , which is indeed the behavior observed in our model (see Table 2).

This naturally leads to the question of whether there is any physical reason that extraordinary mode photons would actually dominate in the X Per accretion column. We argue in Section 9 that “mode pumping” of blackbody seed photons from the ordinary mode to the extraordinary mode will tend to reinforce the dominance of extraordinary mode photons in the X Per accretion column, which supports the super-Thomson cross section results we obtain when applying our model to that source. We also note the surprising similarity between the values obtained for $\sigma_{{}_{||}}$ and $\overline{\sigma}$ in the case of X Per, and suggest that this effect is due to a strong anisotropy in the radiation field inside the accretion column. The anisotropy would result in the beaming of most of the X-rays along the column axis, which is consistent with the dominance of the column-top emission over the column-wall emission in X Per.

9 MODEL SELF-CONSISTENCY

The fundamental physical parameters in the theoretical model developed here are the stellar mass, set using $M_{*}=1.4\,M_{\odot}$ , the stellar radius, set to the canonical value $R_{*}=10\,$ km, the accretion rate, $\dot{M}$ , set using observational estimates of the X-ray luminosity, combined with Equation (78), and the source distance, $D$ , set using published estimates. The remaining free parameters in the theory are $\alpha$ , $\xi$ , $\psi$ , $B$ , $\Theta_{1}$ , $\Theta_{2}$ , $k_{0}$ , $k_{\infty}$ , $y_{\rm top}$ , and $T_{e}$ . These parameters are varied in order to achieve satisfactory qualitative fits between the theoretical model and the observed X-ray spectra of Her X-1 and X Per, with the results displayed in Figures 7 and 9, respectively.

In addition to the need to fit the X-ray data, the theoretical free parameters can also be constrained by considering the dynamical and thermal structure of the accretion column. This is an important issue, since the velocity profile used in our model, given by Equation (45), is a mathematical result obtained via application of the separation of variables technique to the solution of the transport equation (Equation (25)). Hence the velocity profile used for each of the sources treated here must be physically validated based on the values selected for the constants $k_{0}$ and $k_{\infty}$ . We examine this issue below by performing an analysis of the hydrodynamic and thermodynamic structure of the accretion columns in Her X-1 and X Per obtained using our model.

We must also critically examine the results we have obtained for the electron scattering cross sections $\sigma_{\perp}$ , $\sigma_{{}_{||}}$ , and $\overline{\sigma}$ when using our model to analyze the X-ray spectra of Her X-1 and X Per. The results obtained for the cross sections are radically different for the two sources (see Table 2). In the application to Her X-1, we obtain $\sigma_{\perp}/\sigma_{{}_{\rm T}}=0.537$ , $\sigma_{{}_{||}}/\sigma_{{}_{\rm T}}=6.68\times 10^{-5}$ , and $\overline{\sigma}/\sigma_{{}_{\rm T}}=5.90\times 10^{-4}$ , and in the case of X Per, we obtain $\sigma_{\perp}/\sigma_{{}_{\rm T}}=848$ , $\sigma_{{}_{||}}/\sigma_{{}_{\rm T}}=2.46$ , and $\overline{\sigma}/\sigma_{{}_{\rm T}}=4.31$ . We argue below that these apparently contradictory results can be understood as a consequence of the effect of vacuum polarization in X Per, which strongly alters the entire X-ray continuum in that source, whereas this effect is unimportant in the case of Her X-1.

9.1 Column Structure

The dynamical results for the velocity profiles in Her X-1 and X Per, plotted in Figures 8 and 10, respectively, are based on the substitution of specific values for the parameters $k_{0}$ and $k_{\infty}$ into the velocity function given by Equation (45). The values for $k_{0}$ and $k_{\infty}$ are obtained as part of the spectral fitting process for the two sources. While the resulting spectral fits show good qualitative agreement with the X-ray spectra plotted in Figures 7 and 9, it is important to investigate whether the velocity profiles utilized in the models for the two sources satisfy the physical hydrodynamic and thermodynamic conservation equations.

9.1.1 Thermodynamic Structure

In order to validate the input dynamical profiles for Her X-1 and X Per, we need to analyze the structure of the accretion columns by solving a set of physical conservation equations. The results of this calculation will be a “hydrodynamical” velocity profile for each source that can be compared with the input profiles used in our model, which are computed using Equation (45). The thermodynamic and hydrodynamic structures are coupled since the pressure of the radiation field and the gas may each make important contributions. We shall begin by considering the thermodynamic structure of the accretion column by examining the variation of the electron and proton temperatures. In the case of Her X-1, we would expect these two temperatures to track closely, since Coulomb coupling will be very efficient in the dense plasma. However, in the case of X Per, the situation may be quite different. In particular, Langer & Rappaport (1982) have shown that the plasma in the accretion column in X Per is so tenuous that a two-temperature plasma may result, in which the proton temperature, $T_{p}$ , is much higher than the electron temperature, $T_{e}$ , due to inefficient Coulomb coupling between the two species.

In our model, the electrons are assumed to form an isothermal distribution, with a temperature, $T_{e}$ , that is regulated mainly by thermal Comptonization (Lyubarskii & Syunyaev, 1982; Sunyaev & Titarchuk, 1980; West et al., 2017a, b). The isothermal assumption is likely to be approximately satisfied in both the Her X-1 and X Per accretion columns. An analysis of the variation of the proton temperature, $T_{p}$ , requires the solution to an energy equation for the protons, and we can expect to obtain very different results in the cases of Her X-1 and X Per due to the vastly different plasma densities in the two sources. We can write the proton energy equation as

\upsilon\frac{dU_{p}}{dR}=\gamma_{p}\,\frac{U_{p}}{\rho}\,\upsilon\frac{d\rho}{dR}+\dot{U}_{p}\ ,

(185)

where $R$ is the radius from the center of the star, $\rho=m_{p}\,n_{e}$ is the mass density, $\gamma_{p}=5/3$ denotes the proton adiabatic index, and the proton energy density, $U_{p}$ , is related to the proton temperature, $T_{p}$ , via

U_{p}=\frac{1}{\gamma_{p}-1}\,\frac{\rho\,kT_{p}}{m_{p}}\ .

(186)

The first term on the right-hand side of Equation (185) represents adiabatic compression, and the second term, $\dot{U}_{p}$ , denotes the heating rate for the protons via Coulomb coupling with the electrons. Based on Equation (C1) from Langer & Rappaport (1982), we can write

\dot{U}_{p}=-\beta_{e}\,R_{0}\,n_{e}^{2}\,\frac{m_{e}}{m_{p}}\left(\frac{8}{\pi}\right)^{1/2}\frac{T_{p}-T_{e}}{T_{\rm eff}}\left(\frac{m_{e}c^{2}}{kT_{\rm eff}}\right)^{1/2}\ln\Lambda_{\rm QM}\ ,

(187)

where the effective temperature, $T_{\rm eff}$ , is given by

T_{\rm eff}\equiv T_{e}+\frac{m_{e}}{m_{p}}\,T_{p}\ ,

(188)

the constant $R_{0}$ is defined by

R_{0}\equiv\frac{3}{4}\sigma_{{}_{\rm T}}m_{e}c^{3}\ ,

(189)

and $\Lambda_{\rm QM}$ is the Coulomb logarithm, computed using

\Lambda_{\rm QM}=5.41+\frac{1}{4}\,\ln\left(\frac{kT_{\rm eff}}{20\,{\rm keV}}\,\frac{B}{10^{12}\,{\rm G}}\,\frac{10^{20}\,{\rm cm^{-3}}}{n_{e}}\right)\ .

(190)

Note that the heating rate for the protons, $\dot{U}_{p}$ , is negative if $T_{p}>T_{e}$ . The quantity $\beta_{e}$ appearing in Equation (187) is an efficiency factor for Coulomb cooling. In a collisional plasma, such as expected in Her X-1, we set $\beta_{e}=1$ , but in a collisionless plasma, such as expected in X Per, the efficiency of the Coulomb cooling may be reduced significantly below unity, and we therefore set $\beta_{e}=0.1$ in this situation. Physically, this reflects the fact that the mean free path for proton-electron collisions is much larger than the system size in a collisionless plasma, and the proton distribution may therefore deviate significantly from an isotropic Maxwellian (Spitzer, 1962; Kirk & Galloway, 1982).

We can transform Equation (185) into an ordinary differential equation for the proton temperature, $T_{p}$ , by combining it with Equation (186), which yields

\frac{d\ln T_{p}}{dR}=-(\gamma_{p}-1)\frac{d\ln(R^{2}\,|\upsilon|)}{dR}+\frac{1}{\upsilon}\,\frac{\dot{U}_{p}}{U_{p}}\ ,

(191)

where we have also used Equation (24) to relate the mass density $\rho$ to the velocity $\upsilon$ using

\rho=\frac{\dot{M}}{\Omega R^{2}|\upsilon|}\ .

(192)

By making use of Equation (45) for the flow velocity $\upsilon$ , we can simplify Equation (191) to obtain the equivalent form

\frac{d\ln T_{p}}{dy}=-\left[\frac{3(\gamma_{p}-1)R_{g}\,c^{2}\,k_{\infty}^{2}}{R_{*}}+\frac{\dot{M}}{\Omega R_{*}m_{p}}\,\frac{\dot{U}_{p}}{n_{e}^{2}\,kT_{p}}\right]\frac{1}{\upsilon^{2}\,y^{2}}\ ,

(193)

where we have also transformed the independent variable from $R$ to $y\equiv R/R_{*}$ , and $\dot{U}_{p}$ is evaluated using Equation (187).

In order to numerically integrate Equation (193) to obtain the global solution for the proton temperature, we also need to specify a boundary value for $T_{p}$ . We do this in two different ways for the two sources treated here. In the case of Her X-1, we simply set $T_{p}=T_{e}$ at the top of the accretion column. In the case of X Per, we set $T_{p}$ equal to the post-shock value, $T_{p+}$ , obtained on the downstream side of the standing shock located at the top of the column. The standard Rankine-Hugoniot conditions give the downstream proton temperature as

kT_{p+}=m_{p}\upsilon_{-}^{2}\frac{[2+(\gamma_{p}-1){\cal M}^{2}][1+\gamma_{p}(2{\cal M}^{2}-1)]}{\gamma_{p}(1+\gamma_{p})^{2}{\cal M}^{4}}\ ,

(194)

where $\upsilon_{-}$ is the upstream velocity, which is equal to the Newtonian free-fall velocity, and ${\cal M}$ is the shock Mach number. Setting $\gamma_{p}=5/3$ and assuming that ${\cal M}\gg 1$ for a strong shock yields

kT_{p+}=\frac{3}{16}\,m_{p}\upsilon_{-}^{2}\ ,

(195)

which agrees with Equation (16c) from Langer & Rappaport (1982). With the boundary value for $T_{p}$ specified, numerical integration of Equation (193) yields the global solution for $T_{p}$ as a function of the dimensionless radius $y$ .

9.1.2 Hydrodynamic Structure

Once the global solution for the proton temperature, $T_{p}$ , is determined via numerical integration of Equation (193), we can solve for the proton pressure, $P_{p}$ , using the ideal gas relation

P_{p}=\frac{\rho\,kT_{p}}{m_{p}}=\frac{\dot{M}}{\Omega R^{2}|\upsilon|}\,\frac{kT_{p}}{m_{p}}\ ,

(196)

where we have substituted for the mass density, $\rho$ , using Equation (192), and $\upsilon$ is computed using Equation (45). We can now use the proton pressure distribution to obtain the solution for the “hydrodynamical” flow velocity, denoted by $\upsilon^{*}$ , which satisfies the hydrodynamical acceleration equation

\upsilon^{*}\frac{d\upsilon^{*}}{dR}=-\frac{GM_{*}}{R^{2}}-\frac{1}{\rho}\,\frac{dP_{p}}{dR}-\frac{1}{\rho}\,\frac{dP_{e}}{dR}-\frac{1}{\rho}\,\frac{dP_{r}}{dR}\ ,

(197)

where $P_{r}$ denotes the radiation pressure, computed using the solution for the radiation spectrum, and the electron pressure, $P_{e}$ , is given by

P_{e}=n_{e}kT_{e}\ .

(198)

Since the electron temperature, $T_{e}$ , is a constant in our model, it follows that the electron pressure, $P_{e}$ , varies in proportion to the electron number density, $n_{e}$ . By combining results, we can numerically integrate Equation (197) to determine the hydrodynamical flow velocity, $\upsilon^{*}$ , associated with the various pressure gradients acting on the gas in our model. Ideally, we would find that $\upsilon^{*}=\upsilon$ , where $\upsilon$ is the input velocity for the model computed using Equation (45) based on the selected values for the constants $k_{0}$ and $k_{\infty}$ . In practice, we can’t expect perfect agreement between the two velocity profiles, but by making a comparison between the input velocity $\upsilon$ and the hydrodynamical velocity $\upsilon^{*}$ , we can evaluate the validity of our models for the accretion columns in Her X-1 and X Per.

In Figure 11 we plot the proton temperature distributions obtained in our models for the Her X-1 and X Per accretion columns. As expected, in the case of Her X-1, the proton temperature tracks the electron temperature closely, due to efficient Coulomb coupling. Since the electrons are assumed to be isothermal, the resulting proton temperature is also constant. On the other hand, in the case of X Per, we note that the proton temperature is significantly higher than the electron temperature, due to the relative inefficiency of the Coulomb coupling between the protons and the electron in the tenuous plasma in this source (Langer & Rappaport, 1982). In Figure 12 we plot the solutions for the hydrodynamical flow velocity, $\upsilon^{*}$ , obtained for each source by numerically integrating Equation (197), and compare them with the corresponding model flow velocity profiles, $\upsilon$ , computed using Equation (45). The reasonably close agreement between $\upsilon^{*}$ and $\upsilon$ indicates that the hydrodynamics of the models for the accretion columns in Her X-1 and X Per are reasonably self-consistent.

9.2 Electron Scattering Cross Sections

The results we have obtained for the electron scattering cross sections $\sigma_{\perp}$ , $\sigma_{{}_{||}}$ , and $\overline{\sigma}$ in our model are very different in the cases of Her X-1 and X Per, as indicated in Table 2. In the case of Her X-1, the results obtained are $\sigma_{\perp}/\sigma_{{}_{\rm T}}=0.537$ , $\sigma_{{}_{||}}/\sigma_{{}_{\rm T}}=6.68\times 10^{-5}$ , and $\overline{\sigma}/\sigma_{{}_{\rm T}}=5.90\times 10^{-4}$ , and in the case of X Per, the results obtained are $\sigma_{\perp}/\sigma_{{}_{\rm T}}=848$ , $\sigma_{{}_{||}}/\sigma_{{}_{\rm T}}=2.46$ , and $\overline{\sigma}/\sigma_{{}_{\rm T}}=4.31$ . We explore this issue in detail here, and we will show that the disparity between the cross section values in the two sources stems from the strong effect of vacuum polarization in the accretion column in X Per. The low plasma density near the base of the accretion column in X Per, with $n_{e}\sim 10^{18}\,{\rm cm}^{-3}$ , results in a very low value for the vacuum energy, $\epsilon_{\rm vac}\sim 0.1\,$ keV (see Equation (11)). Conversely, the electron number density near the base of the accretion column in Her X-1 is much larger, $n_{e}\sim 10^{24}\,{\rm cm}^{-3}$ , and this yields for the vacuum energy $\epsilon_{\rm vac}\sim 100\,$ keV. Since the effect of vacuum polarization is only important for photons with energy $\epsilon\gtrsim\epsilon_{\rm vac}$ , we conclude that vacuum polarization is completely unimportant for the X-ray continuum in Her X-1 due to the much larger plasma density in that source.

Vacuum polarization not only affects the magnitudes of the total scattering cross sections, but it also has a profound effect on the mode-change cross sections that govern the transition between the ordinary and extraordinary modes that can occur in a single scattering. We plot the energy and angular dependences of the mode-change cross sections in the accretion columns in Her X-1 and X Per using Equation (29) from Meszaros & Ventura (1979) in Figure 13. The left-hand panel in Figure 13 utilizes parameters appropriate for the surface of the thermal mound in the Her X-1 accretion column, with electron number density $n_{e}=1.36\times 10^{24}\,{\rm cm}^{-3}$ and magnetic field strength $B=3.26\times 10^{12}\,$ G. The right-hand panel utilizes parameters appropriate for the thermal mound in the X Per accretion column (which is located at the stellar surface), with $n_{e}=1.47\times 10^{18}\,{\rm cm}^{-3}$ and $B=3.00\times 10^{12}\,$ G. Figure 13 provides a detailed description of the nature of the mode-to-mode transitions that can take place inside the accretion columns in the two sources.

According to the mode-change scattering cross sections for Her X-1 plotted on the left-hand side of Figure 13, when vacuum polarization is unimportant, as in the case of Her X-1, both the ordinary and extraordinary mode photons will tend to scatter coherently, meaning that they are unlikely to change mode in a single scattering. We recall that in the case of Her X-1, the observed spectrum is dominated by Comptonized bremsstrahlung emission, which produces ordinary-mode photons because it is a non-resonant process. Since the bremsstrahlung seed photons are emitted in the ordinary mode, and they tend to scatter in a mode-coherent manner, this suggests that the photon population in the Her X-1 accretion column is dominated by ordinary mode photons. In this case, the dependence on the propagation angle $\theta$ for photon energy $\epsilon\ll\epsilon_{c}$ indicated by the red dashed curve on the left-hand side of Figure 3 implies that in the perpendicular direction, photons will scatter with cross section $\sigma_{\perp}\sim\sigma_{{}_{\rm T}}$ , and in the parallel direction, they will scatter with cross section $\sigma_{{}_{||}}\sim 10^{-3}\,\sigma_{{}_{\rm T}}$ . These cross section estimates agree reasonably well with the cross section values obtained in our model for Her X-1 (see Table 2).

The mode-change scattering cross sections for X Per plotted on the right-hand side of Figure 13 display a qualitatively different behavior than the results for Her X-1. In particular, we note that ordinary mode photons with propagation angle $\theta\lesssim 30^{\circ}$ and energy $\epsilon\sim 0.5\,\epsilon_{c}$ (denoted by the green curves) are very likely to mode-change into extraordinary mode photons. However, the corresponding extraordinary mode photons (denoted by the blue curves) are likely to remain as extraordinary mode photons, regardless of the value of $\theta$ . These facts combine to suggest that “mode pumping” may occur in the X Per accretion column, which is a process in which the photons tend to become increasingly concentrated in the extraordinary mode. The regions of mode pumping are indicated by the magenta ellipses in Figure 13, which are the regions where the dotted curves (mode-change cross section) lie above the dashed curves (mode-coherent cross section). This process is physically significant, because the buildup of photons in the extraordinary mode leads to super-Thomson values for the total scattering cross section for photons with energy $\epsilon\sim\epsilon_{c}$ , with a dependence on $\theta$ that yields values for the perpendicular cross section on the order of $\sigma_{\perp}\sim 10^{3}\,\sigma_{{}_{\rm T}}$ and values for the parallel cross section on the order of $\sigma_{{}_{||}}\sim\sigma_{{}_{\rm T}}$ , according to the solid blue curve on the right-hand side of Figure 3. As indicated in Table 2, these cross section values agree reasonably well with the results obtained in our model for X Per. Note that super-Thomson cross section values for X Per are apparent in Figure 3, but not in Figure 13, because the photon energies considered in Figure 3 are closer to $\epsilon_{c}$ .

We must also consider the implications of the variation of the dipole magnetic field for the cyclotron energy, $\epsilon_{c}$ , and how this may cause the resonance in the cross section to sweep through the X-ray continuum for each of the two sources considered here. In connection with this, we point out that the observed X-ray spectrum in X Per is dominated by the column-top component. Hence most of the photons experience a large transition in radius, from $R\sim R_{*}$ at the stellar surface, up to the column top located at radius $R\sim 2\,R_{*}$ . For a dipole magnetic field, the increase in radius by a factor of 2 from the bottom to the top of the accretion column implies a drop in the magnetic field strength, $B$ , (and hence the cyclotron energy $\epsilon_{c}$ ) by a factor of 8, which will sweep the cyclotron resonance into the center of the X-ray continuum for X Per. This naturally leads to the super-Thomson scattering cross sections for this source indicated in Table 2. On the other hand, in the case of Her X-1, the escaping radiation is dominated by the column-wall component, and most of the photons escape near the base of the column. It follows that there is not much variation in $B$ or $\epsilon_{c}$ throughout the emission region in this source. Hence, in the case of Her X-1, the cyclotron resonance will not sweep through the continuum, and therefore most of the X-ray photons have energy $\epsilon\ll\epsilon_{c}$ , which is consistent with the sub-Thomson values for the scattering cross sections listed in Table 2.

10 DISCUSSION AND CONCLUSION

We have developed a new model for the radiative transfer occurring in accretion-powered X-ray pulsars based on a conical geometry for the accretion column, that also incorporates a number of additional enhancements relative to the previous model presented by the authors (BW07). In particular, the new model utilizes a very flexible form for the flow velocity profile (Equation (45)) that includes two free parameters, $k_{0}$ and $k_{\infty}$ , where $k_{0}$ represents the ratio of the flow velocity divided by the Newtonian free-fall velocity at the stellar surface, and $k_{\infty}$ is the same ratio but evaluated as $R\to\infty$ (see Equations (46) and (47)). Many different scenarios can be modeled using this approach. For example, if we set $k_{0}=0$ and $k_{\infty}=1$ , then the velocity equals zero at the stellar surface, and merges smoothly with the Newtonian free-fall profile far from the star. This velocity profile, depicted in Figure 8, is probably an appropriate choice for luminous sources such as Her X-1, in which radiation pressure decelerates the matter to rest at the base of the accretion column. On the other hand, in low-luminosity sources such as X Per, radiation pressure is probably insufficient to decelerate the material to rest at the stellar surface, and the gas may impact onto the star with a substantial residual velocity, after passing through a standing, gas-mediated shock at the top of the column. We use this scenario to model the X Per accretion column, depicted in Figure 10, using the parameter values $k_{0}\sim 0.3$ and $k_{\infty}\sim 0.25$ . In this case, the surface impact velocity is $\upsilon_{*}=-0.19\,c$ and the velocity profile at the top of the column is equal to $1/4$ the free-fall value, in accordance with the strong-shock jump condition.

The central result of the paper is the closed-form analytical solution for the Green’s function describing the spectrum of the radiation inside the accretion column as a function of the scattering optical depth, $\tau$ , and the photon energy, $\epsilon$ , given by Equation (114). This solution represents the response to the continual injection of monoenergetic seed photons with energy $\epsilon_{0}$ from a source located at optical depth $\tau_{0}$ inside a conical accretion column with velocity profile given by Equation (45). Based on this result, we proceed to derive the corresponding Green’s functions for the column-integrated spectrum escaping through the column walls (Equation (126)), and also the Green’s function for the spectrum escaping through the column top (Equation (136)). Since the fundamental transport equation (Equation (25)) includes an escape term to describe the diffusion of photons through the column walls, and the upper surface (the column top) is treated as a free-streaming surface, it follows that our results for the two escaping radiation components are computed self-consistently.

The model includes the implementation of source photon distributions corresponding to blackbody, cyclotron, and bremsstrahlung emission. We confirm that the number of photons escaping from the column per unit time is exactly equal to the number injected, for each radiation mechanism, as required in the steady-state scenario considered here. The contributions to the observed X-ray spectrum due to emission from the walls and top of the accretion column are computed using Equations (177) and (180), respectively, and the approximate phase-averaged spectrum is computed using Equation (182). In Figures 7 and 9, we display the qualitative fits to the phase-averaged X-ray spectra of Her X-1 and X Per, respectively, obtained using the new model, and we note that the agreement between the theory and the data is reasonably close for both sources. We find that in the case of Her X-1, the observed phase-averaged X-ray spectrum is dominated by Comptonized bremsstrahlung emission escaping through the side walls of the accretion column, and in the case of X Per, it is dominated by Comptonized blackbody emission escaping primarily through the top of the column. The associated velocity profiles for the two sources are presented in Figures 8 and 10, respectively, and these profiles are validated via comparison with detailed hydrodynamical calculations in Figure 12. The set of results presented here represents the first time that a single overarching model for spectral formation in X-ray pulsars has been able to successfully reproduce the observed X-ray spectra for two sources spanning about three orders of magnitude in X-ray luminosity, from $L_{X}\sim 10^{34}\,{\rm ergs\,s}^{-1}$ for X Per to $L_{X}\sim 10^{37}\,{\rm ergs\,s}^{-1}$ for Her X-1.

In the case of Her X-1, the velocity profile plotted in Figure 8 begins with a Newtonian free-fall shape far from the star, and smoothly transitions into a radiative, radiation-dominated shock as the gas approaches the base of the accretion column. The radiation-dominated shock decelerates the gas to rest at the stellar surface, and the kinetic energy of the flow is radiated away primarily through the walls of the accretion column. On the other hand, the velocity profile for X Per (Figure 10) exhibits a very different shape. In this source, the flow passes through a standing, gas-mediated discontinuous shock at the top of the accretion column, where the velocity drops by a factor of four, after which it proceeds to reaccelerate as the gas approaches the stellar surface. The gas collides with the neutron star with a residual velocity $\sim 0.19\,c$ (see Table 2).

10.1 Model Parameters

The process of obtaining the qualitative fits to the X-ray spectra displayed in Figures 7 and 9 begins with the selection of values for the stellar mass, $M_{*}$ , the stellar radius, $R_{*}$ , the accretion rate, $\dot{M}$ , and the source distance, $D$ . We use canonical values for the stellar mass and radius, with $M_{*}=1.4\,M_{\odot}$ and $R_{*}=10\,$ km, and the values for the distance $D$ and the accretion rate $\dot{M}$ are taken from published estimates. After these parameters are set, the spectral fits are obtained by varying the fundamental free parameters $\alpha$ , $\xi$ , $\psi$ , $B$ , $k_{0}$ , $k_{\infty}$ , $\Theta_{1}$ , $\Theta_{2}$ , $T_{e}$ , and $y_{\rm top}$ . The process for determining the values for the free parameters includes a comparison with the observed X-ray spectrum, as well as a consideration of the thermodynamic and hydrodynamic structure of the accretion column, as discussed in Section 9. Once satisfactory qualitative spectral fits are obtained, and the structure of the column is deemed sufficiently self-consistent, then we can compute the associated values for the scattering cross sections $\sigma_{{}_{||}}$ , $\sigma_{\perp}$ , and $\overline{\sigma}$ using Equation (139), (140), and (141), respectively. The results obtained for the fundamental free parameters for Her X-1 and X Per are listed in Table 1, and the corresponding values obtained for the scattering cross sections are reported in Table 2.

The values for the scattering cross sections $\sigma_{{}_{||}}$ , $\sigma_{\perp}$ , and $\overline{\sigma}$ obtained by applying our model to Her X-1 and X Per are very different, and this has motivated further investigation. We argue that the difference in the cross sections stems from the effect of vacuum polarization, which is negligible in the case of Her X-1, but very important in the application to X Per, due to the different values for the electron number density, $n_{e}$ , in the two sources. In X Per, the density at the thermal mound surface is $n_{e}\sim 10^{18}\,{\rm cm}^{-3}$ , whereas in the case of Her X-1, we find that $n_{e}\sim 10^{24}\,{\rm cm}^{-3}$ . According to Equation (11), the corresponding results obtained for the vacuum energy, $\epsilon_{\rm vac}$ , in the two sources are $\epsilon_{\rm vac}\sim 0.1\,$ keV for X Per and $\epsilon_{\rm vac}\sim 100\,$ keV for Her X-1. Since vacuum polarization profoundly influences the electron scattering cross sections for photons with energy $\epsilon\gtrsim\epsilon_{\rm vac}$ , it is clear the this process plays a crucial role in determining the shape of the X-ray continuum in X Per. Conversely, in the case of Her X-1, vacuum polarization is unimportant.

The detailed implications of vacuum polarization are explored in Section 9, where we plot the mode-change cross sections for extraordinary and ordinary mode photons in Figure 13, using parameters appropriate for X Per and Her X-1. The plots clearly indicate that in the case of Her X-1, most of the photons will remain in the ordinary mode, because bremsstrahlung emission (which dominates the seed photon production in Her X-1) primarily produces ordinary mode photons, and there is little propensity for these photons to switch to the extraordinary mode, as indicated on the left-hand side of Figure 13. On the other hand, in the case of X Per, the situation is quite different. In particular, in this source, there is a strong likelihood for ordinary mode photons to switch into the extraordinary mode, and this will lead to “mode pumping,” in which photons tend to accumulate in the extraordinary mode, as indicated on the right-hand side of Figure 13.

The concentration of photons in the extraordinary polarization mode in X Per, combined with the fact that the cyclotron resonance will sweep through the continuum in this source due to the escape of most of the photons through the top of the column (at radius $R_{\rm top}\sim 2\,R_{*}$ ), will tend to increase the cross sections, yielding the super-Thomson values listed in Table 2 for this source. Hence in the case of X Per, we find that $\sigma_{\perp}\gg\overline{\sigma}\gtrsim\sigma_{{}_{||}}\gtrsim\sigma_{{}_{\rm T}}$ . On the other hand, in the case of Her X-1, the cyclotron resonance will not sweep through the continuum as strongly, because most of the radiation leaks out through the column walls near the base of the flow. The resulting cross sections in this source are generally sub-Thomson, and we find that $\sigma_{||}\ll\overline{\sigma}\ll\sigma_{\perp}\sim\sigma_{{}_{\rm T}}$ for Her X-1. The similarity between the values for $\sigma_{{}_{||}}$ and $\overline{\sigma}$ obtained in the application to X Per, along with the much larger value obtained for $\sigma_{\perp}$ , suggests a strong anisotropy in the radiation field, with most of the X-rays beamed along the column axis. We note that this behavior is consistent with the dominance of the column-top emission versus the wall emission that we observe in our model for X Per (see Figure 9).

10.2 Dependence of Spectrum on Luminosity

The cyclotron resonance absorption feature in the X-ray spectrum of Her X-1 exhibits an apparent dependence on the X-ray luminosity which has been studied by several authors (e.g., Staubert et al., 2019, 2020). In general, the observations indicate that the power-law index of the continuum emission reduces, and the spectrum becomes harder, as the luminosity increases. The same behavior is observed in both the long-term variability of the source, as well as on pulse-to-pulse timescales (Klochkov et al., 2011). In the context of the model developed here, this type of spectral variability with changing luminosity is a direct consequence of changes in the radiation pressure, which ultimately controls the impact velocity of the gas onto the surface of the neutron star. In particular, an increase in the luminosity will amplify the radiation pressure, leading to a decrease in the impact velocity at the stellar surface. The decreased surface velocity leads to an increase in the gas density, and the resulting enhanced compression of the gas acts as a “piston,” which powers an increase in the mechanical $P$ d $V$ work done on the radiation field by the gas via bulk Comptonization. The increased work performed on the radiation field naturally leads to a flatter, harder spectrum in the high-luminosity sources.

The scenario described above is appropriate for supercritical sources, such as Her X-1, but it must be modified in the case of subcritical sources, such as X Per, in which radiation pressure is probably insufficient to accomplish complete deceleration of the gas to rest at the stellar surface. Instead, the gas may ram into the star with a significant residual velocity, perhaps as large as $\sim 0.25\,c$ . In this situation, the gas will do significantly less $P$ d $V$ work on the radiation field, due to the smaller amount of compression relative to what would occur if the gas were decelerated completely to rest at the stellar surface. The decreased amount of compression implies a decreased amount of bulk Comptonization, which leads to a steeper spectral slope. In this scenario, the final deceleration occurs in the last few cm above the stellar surface as a result of Coulomb collisions, which heat the material, powering the blackbody hot-spot that dominates the production of seed photons (Sokolova-Lapa et al., 2021). We believe that this simple mechanism may provide a simple physical explanation for the steeper slopes observed in the X-ray continuum for low-luminosity sources such as X Per.

10.3 Conclusion

We have demonstrated that bulk and thermal Comptonization naturally leads to emergent spectra in accretion-powered X-ray pulsars that are in good agreement with the observational data for both high-luminosity and low-luminosity sources. Furthermore, we have shown that the spectrum of the high-luminosity source treated here, Her X-1, is dominated by reprocessed bremsstrahlung emission, escaping primarily through the walls of the accretion column. This result agrees with the findings of West et al. (2017a, b). Conversely, for the low-luminosity source treated here, X Per, we have demonstrated that the observed X-ray spectrum is dominated by Comptonized blackbody emission, escaping through the top of the accretion column, rather than through the walls. In principle, the difference between the dominance of the column-wall emission versus the column-top emission in these two sources could suggest a difference in the pulse profiles and the phase-resolved spectra. We plan to investigate this question in future work, using a general relativistic framework such as the one developed by Riffert & Meszaros (1988).

In recent years, there has been a significant increase in the development of space-based X-ray polarimetry observatories designed to study the polarization properties of the high-energy emission observed from compact sources, including accretion-powered X-ray pulsars. These missions include NASA’s IXPE which launched in 2021, and the Indian mission XPoSat, with a planned launch in 2022. The availability of these new observatories complements recent advances in theoretical work, which suggest that there may be a complex, phase-dependent signature of polarization contained in the radiation detected from X-ray pulsars (e.g., Caiazzo & Heyl, 2021a, b). In fact, a possible detection of polarization in the X-ray spectrum of Her X-1 has already been reported by Doroshenko et al. (2022). We note that it may be possible to develop new predictions for X-ray polarization as a function of pulse phase and photon energy using the new model developed here, although such a calculation is beyond the scope of this paper. However, we can make a straightforward prediction regarding the differences in the overall polarization properties of the spectra emitted by Her X-1 and X Per, since our model implies that the emission from Her X-1 is dominated by ordinary mode photons, and the emission from X Per is dominated by extraordinary mode photons, as viewed in the frame of the star. However, in the frame of the Earth, the situation is complicated by the spin of the star, which will tend to homogenize the polarization properties between the extraordinary and ordinary modes (Doroshenko et al., 2022).

We also note the possible existence of an additional high-energy component in the spectrum of X Per, which manifests as a broad hump centered at photon energy $\sim 70\,$ keV (Doroshenko et al., 2012). This has been interpreted as a possible consequence of the emission of cyclotron photons produced via the excitation of electrons into the $n=1$ Landau state due to collisions with protons in the Coulomb-stopping layer just above the stellar surface. However, no-self consistent model for this process has been developed yet (Mushtukov & Tsygankov, 2022). In future work, we intend to investigate this scenario, as well as the alternative possibility that the excess high-energy emission is directly connected with the presence of the standing shock at the top of the accretion column (Langer & Rappaport, 1982).

The authors are grateful for useful discussions with Ekaterina Sokolova-Lapa, Joern Wilms, Katja Pottschmidt, Kent Wood, Jason Wong, Brent West, and Ken Wolfram. This work was supported in part by NASA through the NuSTAR Guest Observer Program and the Astrophysics Explorers Program. The authors would also like to acknowledge helpful comments from the anonymous referee. This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. We acknowledge use of NASA’s Astrophysics Data System (ADS) bibliographic services and the ArXiv.

Model	$\sigma_{\perp}/\sigma_{{}_{\rm T}}$	$\sigma_{{}_{\|\|}}/\sigma_{{}_{\rm T}}$	$\overline{\sigma}/\sigma_{{}_{\rm T}}$	$r_{0}\,$ (m)	$z_{\rm th}\,$ (cm)	$T_{\rm th}\,$ (K)	$\upsilon_{*}/c$	$\upsilon_{\rm th}/c$	$\tau_{\rm th}$	$\tau_{\rm top}$
1	$0.537$	$6.68\times 10^{-5}$	$5.90\times 10^{-4}$	$55.0$	692	$6.07\times 10^{7}$	0.0	-0.03	0.08	2.91
2	$848.0$	$2.46$	$4.31$	$966.0$	0.0	$8.04\times 10^{6}$	-0.19	-0.19	0.0	1.75

Model	Object	$\alpha$	$\xi$	$\psi$	$k_{0}$	$k_{\infty}$	$\Theta_{1}$	$\Theta_{2}$	$y_{\rm top}$	$\!\!\!\dot{M}$ (g s^-1)	$T_{e}$ (K)	$B$ (G)	$D$ (kpc)
1	Her X-1	0.35	1.14	1.42	0.0	1.0	0.315	0	2.15	$1.90\times 10^{17}$	$5.50\times 10^{7}$	$3.26\times 10^{12}$	6.60
2	X Per	0.08	1.20	0.75	0.3	0.244	5.536	0	$2.20$	$4.16\times 10^{14}$	$1.20\times 10^{8}$	$3.00\times 10^{12}$	0.725

A GENERALIZED ANALYTICAL MODEL FOR THERMAL AND BULK COMPTONIZATION IN ACCRETION-POWERED X-RAY PULSARS

Abstract

1 INTRODUCTION

2 RADIATIVE PROCESSES

2.1 Radiative Transfer in Magnetized Media

2.1.1 Electron Scattering Cross Sections in Magnetized Plasma

2.1.2 Effects of Vacuum Polarization

2.1.3 Simplified Scattering Cross Sections

3 RADIATIVE TRANSFER

3.1 Transport Equation

3.2 Separability

3.3 Optical Depth Variation

4 SOLUTION FOR THE GREEN’S FUNCTION

4.1 Spatial Boundary Conditions

4.1.1 Free-Streaming Upper Boundary Condition

4.1.2 Zero Diffusion Flux Lower Boundary Condition

4.1.3 Zero Number Flux Lower Boundary Condition

4.1.4 Zero Energy Flux Lower Boundary Condition

4.2 Eigenvalues and Spatial Eigenfunctions

4.2.1 Upper Boundary Condition

4.2.2 Lower Boundary Condition

4.2.3 Spatial Eigenfunctions

4.3 Eigenfunction Orthogonality

4.4 Energy Eigenfunctions

4.5 Green’s Function Expansion

4.6 Radiation Distribution for Arbitrary Source Term

5 EMITTED RADIATION SPECTRUM

5.1 Green’s Function for Column Wall Spectrum

5.2 Wall Spectrum for a General Source

5.3 Green’s Function for Column Top Spectrum

5.4 Column Top Spectrum for a General Source

6 MODEL PARAMETERS AND CONSTRAINTS

6.1 Scattering Cross Sections

6.2 Thermal Mound Properties

7 PHOTON SOURCES AND ASSOCIATED SPECTRA

7.1 Seed Photon Source Functions

7.2 Cyclotron Radiation

7.3 Blackbody Radiation

7.4 Bremsstrahlung Radiation

8 ASTROPHYSICAL APPLICATIONS

8.1 Her X-1

8.2 X Per

9 MODEL SELF-CONSISTENCY

9.1 Column Structure

9.1.1 Thermodynamic Structure

9.1.2 Hydrodynamic Structure

9.2 Electron Scattering Cross Sections

10 DISCUSSION AND CONCLUSION

10.1 Model Parameters

10.2 Dependence of Spectrum on Luminosity

10.3 Conclusion

References