A New Two-fluid Radiation-hydrodynamical Model for X-Ray Pulsar Accretion Columns

The ultimate power source for the observed X-ray emission from accretion-powered pulsars is gravity, and therefore the total power available is equal to the accretion luminosity, defined by

$$\begin{eqnarray}&&{L}_{\mathrm{acc}}\equiv \displaystyle \frac{{{GM}}_{* }\dot{M}}{{R}_{* }},\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant, $\dot{M}$ denotes the accretion rate, and M_* and R_* are the stellar mass and radius, respectively. If no kinetic or thermal energy enters the star (Lenzen & Trümper 1978), then the X-ray luminosity L_X is given by the relation

$$\begin{eqnarray}&&{L}_{{\rm{X}}}={L}_{\mathrm{acc}},\end{eqnarray} \tag{ 2 }$$

although we note that Basko & Sunyaev (1976) have argued that some energy may diffuse down into the star.

Becker et al. (2012) have shown there is a critical luminosity, L_crit, below which the pressure of the radiation alone is insufficient to bring the matter to rest, and therefore Coulomb interactions must cause the final deceleration to stagnate at the stellar surface. Additionally, very low-luminosity sources (L_X ≲ 10³⁴ erg s⁻¹) can potentially exhibit the presence of a gas-mediated (discontinuous) shock downstream from the (smooth) radiation shock (Langer & Rappaport 1982). The precise locations of the radiation and gas shocks largely depend upon the source luminosity and the upstream and downstream boundary conditions. Hence, in order to fully understand the hydrodynamic and thermodynamic processes that determine the structure of the accretion column over the full range of observed luminosities (L_X ∼ 10^34–38 erg s⁻¹), it is essential to include the effect of gas pressure in the model. In this paper, we focus on treating three well-known luminous X-ray pulsars, and we defer discussion of low-luminosity sources, such as X Persei, to a later paper.

The magnetic field surrounding a neutron star is well approximated by a dipole configuration, with spherical vector components given by (e.g., Jackson 1962)

$$\begin{eqnarray}&&{B}_{r}=-\displaystyle \frac{{B}_{* }{R}_{* }^{3}\cos \theta }{{r}^{3}},\quad {B}_{\theta }=-\displaystyle \frac{{B}_{* }{R}_{* }^{3}\sin \theta }{2{r}^{3}},\quad {B}_{\phi }=0,\end{eqnarray} \tag{ 3 }$$

where the polar angle θ is measured from the magnetic field axis, B_* denotes the field strength measured at the magnetic pole on the surface of the star, and R_* is the stellar radius. The magnitude of the field, $| {\boldsymbol{B}}|$, varies with the spherical radius r according to

$$\begin{eqnarray}&&| {\boldsymbol{B}}| =\displaystyle \frac{{B}_{* }{R}_{* }^{3}}{2{r}^{3}}\sqrt{1+3{\cos }^{2}\theta },\end{eqnarray} \tag{ 4 }$$

and therefore the field strength decreases by a factor of two between the magnetic pole (θ = 0) and the magnetic equator (θ = π/2), so that

$$\begin{eqnarray}&&{B}_{\mathrm{eq}}=\displaystyle \frac{1}{2}{B}_{* },\end{eqnarray} \tag{ 5 }$$

where B_eq denotes the magnitude of the field at the stellar surface along the magnetic equator.

In the scenario considered here, the accreting gas is entrained onto magnetic field lines from the surrounding disk, and fed onto the magnetic poles of the star. The detailed density distribution inside the accretion column is influenced by a variety of unknown geometrical factors, such as the angle between the star's rotation and magnetic axes (Lamb et al. 1973; Elsner & Lamb 1977; Ghosh et al. 1977). In some sources, the entrainment of matter from the disk results in a partially filled column, but in other sources, such as Her X-1, an alternate accretion mode seems to be at play, in which the gas is introduced into the polar cap region from a dense atmosphere concentrated above the cap, and the accretion column is completely filled (Boroson et al. 2001).

We define the physical extent of the accretion column at the stellar surface using the polar angles θ₁ and θ₂, which are measured from the magnetic field axis and delineate the inner and outer boundaries of the dipole accretion column at the stellar surface, respectively. The corresponding inner and outer arc-length surface radii, denoted by ℓ₁ and ℓ₂, respectively, are given by (see Figure 2)

$$\begin{eqnarray}&&{{\ell }}_{1}={\theta }_{1}\,{R}_{* },\ \ \ {{\ell }}_{2}={\theta }_{2}\,{R}_{* }.\end{eqnarray} \tag{ 6 }$$

Note that the column is partially hollow if 0 < ℓ₁ < ℓ₂, and it is completely filled if ℓ₁ = 0. The solid angle subtended by the accretion column at the stellar surface, Ω_*, is related to θ₁ and θ₂ via

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{* }=2\pi (\cos {\theta }_{1}-\cos {\theta }_{2}).\end{eqnarray} \tag{ 7 }$$

The variable solid angle, Ω(r), subtended by the accretion column at radius r increases in proportion to r in the dipole field geometry, so that

$$\begin{eqnarray}&&{\rm{\Omega }}(r)=\displaystyle \frac{r}{{R}_{* }}\ {{\rm{\Omega }}}_{* }.\end{eqnarray} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

Lamb et al. (1973) provide some insight into the upper limit of the outer polar cap arc-radius ℓ₂. The stellar surface "hot spot" has an area that must be less than or equal to $\pi ({R}_{* }/{R}_{{\rm{A}}}){R}_{* }^{2}$, and therefore

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{2}{R}_{* }^{2}\lesssim \pi \left(\displaystyle \frac{{R}_{* }}{{R}_{{\rm{A}}}}\right){R}_{* }^{2},\end{eqnarray} \tag{ 9 }$$

where R_A is the Alfvén radius and Ω₂ is the solid angle subtended by a filled polar cap of radius ℓ₂ on the surface of the star, given by

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{2}\equiv 2\pi (1-\cos {\theta }_{2}).\end{eqnarray} \tag{ 10 }$$

It follows that the solid angle of the polar cap is restricted by the condition

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{2}\lesssim \pi \,\displaystyle \frac{{R}_{* }}{{R}_{{\rm{A}}}}.\end{eqnarray} \tag{ 11 }$$

Since the stellar radius is much larger than the polar cap radius (R_* ≫ ℓ₂), and the Alfvén radius is much larger that the stellar radius (R_A ≫ R_*), we can use the small angle approximation, sin θ ≈ θ, along with Equations (6), (10), and (11), to conclude that the outer polar cap arc-radius, ℓ₂, is constrained by the condition

$$\begin{eqnarray}&&{{\ell }}_{2}\leqslant {R}_{* }{\left(\displaystyle \frac{{R}_{* }}{{R}_{{\rm{A}}}}\right)}^{1/2}.\end{eqnarray} \tag{ 12 }$$

In the dipole field geometry, the radius r along a field line (which is also a flow streamline in the pulsar application) is a function of the angle θ measured from the magnetic pole via

$$\begin{eqnarray}&&r(\theta )={R}_{\mathrm{eq}}{\sin }^{2}\theta ,\end{eqnarray} \tag{ 13 }$$

where R_eq is the radius of the field line in the magnetic-equatorial plane (θ = π/2). The field lines connected with the inner and outer surfaces of the accretion column have magnetic-equatorial radii R_eq equal to R₁ and R₂, respectively, where R₁ > R₂ (see Figure 3). We can relate R₁ and R₂ to the corresponding polar angles, θ₁ and θ₂, respectively, by setting r = R_* in Equation (13), which yields

$$\begin{eqnarray}&&{R}_{1}=\displaystyle \frac{{R}_{* }}{{\sin }^{2}{\theta }_{1}},\ \ \ {R}_{2}=\displaystyle \frac{{R}_{* }}{{\sin }^{2}{\theta }_{2}}.\end{eqnarray} \tag{ 14 }$$

Zoom In Zoom Out Reset image size

We define the coordinate z as the altitude above the magnetic-equatorial plane, measured along the field line that connects with the outer wall of the accretion funnel and with the inner accretion radius in the Keplerian disk. We therefore have

$$\begin{eqnarray}&&z(\theta )=r(\theta )\cos \theta ={R}_{2}{\sin }^{2}\theta \cos \theta ,\end{eqnarray} \tag{ 15 }$$

where the final result follows from Equation (13). The dipole field reaches its maximum altitude, z_c, at the critical angle θ = θ_c, and then turns over to extend downward toward the disk. By setting the derivative of Equation (15) with respect to θ equal to zero, we find that the critical angle is given by

$$\begin{eqnarray}&&{\theta }_{c}={\cos }^{-1}\left(\displaystyle \frac{1}{\sqrt{3}}\right)=57.{74}^{^\circ }.\end{eqnarray} \tag{ 16 }$$

The corresponding maximum altitude is therefore

$$\begin{eqnarray}&&{z}_{c}={r}_{c}\cos {\theta }_{c}=\displaystyle \frac{2}{3\sqrt{3}}\,{R}_{2},\end{eqnarray} \tag{ 17 }$$

where the corresponding spherical radius, r_c, is related to the magnetic-equatorial plane radius, R₂, via (see Equation (13))

$$\begin{eqnarray}&&{r}_{c}=\displaystyle \frac{2}{3}\ {R}_{2}.\end{eqnarray} \tag{ 18 }$$

A fundamental geometrical restriction of our model is that the spherical radius at the top of the accretion funnel, denoted by r_top, must be below the dipole turnover radius, r_c, associated with the outer-wall field line. Hence we must satisfy the condition

$$\begin{eqnarray}&&{r}_{\mathrm{top}}\leqslant {r}_{c}.\end{eqnarray} \tag{ 19 }$$

The pulsations observed from an X-ray pulsar result from a misalignment between the magnetic and rotation axes of the star. The angle between these two axes is denoted by φ in our model. The misalignment causes the magnetic field at the surface of the star in the plane of the accretion disk, B_disk, to sweep between minimum and maximum values during the star's rotation, as observed from a standard reference direction, which we take to be the direction to the companion star. Based on Equation (4), we find that

$$\begin{eqnarray}&&{B}_{\mathrm{disk}}=\displaystyle \frac{{B}_{* }}{2}\sqrt{1+3{\sin }^{2}\alpha },\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{\pi }{2}-\theta \end{eqnarray} \tag{ 21 }$$

represents the magnetic latitude in the accretion disk (in the direction toward the companion star), which varies between ±φ as the star rotates, such that −φ ≤ α ≤ φ (see Figure 3).

Matter is picked up from the disk and entrained onto the magnetic field lines at the Alfvén radius, R_A, located where the pressure of the magnetic field balances the ram pressure of the accreting gas (Lamb et al. 1973). Outside this radius, the magnetic field of the neutron star is effectively shielded, and therefore it does not significantly influence the flow structure. Inside the Alfvén radius, the strong magnetic field channels the plasma onto the magnetic poles of the star. Due to the complex structure of the pulsar magnetosphere and the uncertainties regarding its interaction with the matter in the disk, it is difficult to precisely compute the value of R_A (e.g., Romanova et al. 2003). However, a useful estimate is provided by Lamb et al. (1973), who find that

$$\begin{eqnarray}{R}_{{\rm{A}}} & \sim & 2.6\times {10}^{8}\,\mathrm{cm}\,{\left(\displaystyle \frac{{B}_{\mathrm{disk}}}{{10}^{12}{\rm{G}}}\right)}^{4/7}{\left(\displaystyle \frac{{R}_{* }}{10\mathrm{km}}\right)}^{10/7}\\ & & \times \,{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{1/7}{\left(\displaystyle \frac{{L}_{{\rm{X}}}}{{10}^{37}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{-2/7}\xi ,\end{eqnarray} \tag{ 22 }$$

where ξ is a constant of order unity. Based on Equations (20) and (22), we observe that the oscillation of the disk-plane surface magnetic field, B_disk, in the direction toward the companion star, will generate a corresponding oscillation in the Alfvén radius, R_A, in the same direction. Since the matter is picked up from the accretion disk and entrained onto the magnetic field lines at radius R_A, it follows that the pick-up radius in the disk oscillates between minimum and maximum values as the star rotates.

We denote the radii where the magnetic field lines connected with the inner and outer walls of the accretion column cross the accretion disk as ${R}_{1,\mathrm{disk}}$, and ${R}_{2,\mathrm{disk}}$, respectively, where ${R}_{1,\mathrm{disk}}\gt {R}_{2,\mathrm{disk}}$. The corresponding radii at which these field lines cross the equatorial plane of the magnetic dipole are R₁ and R₂, respectively. By setting the magnetic-equatorial crossing radius, R_eq, equal to either R₁ or R₂ in Equation (13), we find that the corresponding disk-crossing radii for the two field lines in question are given by

$$\begin{eqnarray}{R}_{1,\mathrm{disk}} & = & {R}_{1}{\sin }^{2}\theta ={R}_{1}{\cos }^{2}\alpha ,\\ {R}_{2,\mathrm{disk}} & = & {R}_{2}{\sin }^{2}\theta ={R}_{2}{\cos }^{2}\alpha ,\end{eqnarray} \tag{ 23 }$$

where α is the magnetic latitude in the accretion disk, and the final results follow from Equation (21). Equation (23) indicates that the disk-crossing radii ${R}_{1,\mathrm{disk}}$ and ${R}_{2,\mathrm{disk}}$ oscillate as the star rotates and α varies between ±φ. By combining Equations (14) and (23), we can eliminate R₁ and R₂ to express the disk-crossing radii in terms of the angles θ₁ and θ₂, which yield

$$\begin{eqnarray}&&{R}_{1,\mathrm{disk}}={R}_{* }{\left(\displaystyle \frac{\cos \alpha }{\sin {\theta }_{1}}\right)}^{2},\qquad {R}_{2,\mathrm{disk}}={R}_{* }{\left(\displaystyle \frac{\cos \alpha }{\sin {\theta }_{2}}\right)}^{2}.\end{eqnarray} \tag{ 24 }$$

Equation (24) allows us to study the variation of the two disk-crossing radii as the star spins and the disk-plane latitude α oscillates between ±φ. This is important because the matter is picked up from the disk at the Alfvén radius, R_A, and therefore material is fed onto the inner and outer walls of the accretion column when ${R}_{1,\mathrm{disk}}={R}_{{\rm{A}}}$ and ${R}_{2,\mathrm{disk}}={R}_{{\rm{A}}}$, respectively. For intermediate values, the matter is fed into the central part of the column, between the inner and outer walls. Hence, as the star rotates, matter is cyclically fed into the entire volume of the accretion column.

In order to close the system and ensure that we are generating self-consistent models for the pulsar accretion column and its connection with the surrounding accretion disk, we must therefore set ${R}_{1,\mathrm{disk}}$ and ${R}_{2,\mathrm{disk}}$ equal to the maximum and minimum values for the oscillating Alfvén radius, which is obtained by combining Equations (22) and (20). Essentially, we must find that during the spin of the star, R_A varies in the range of

$$\begin{eqnarray}&&{R}_{2,\mathrm{disk}}\lesssim {R}_{{\rm{A}}}\lesssim {R}_{1,\mathrm{disk}}.\end{eqnarray} \tag{ 25 }$$

We should emphasize that our model does not include a complete description of the entire pulsar magnetosphere and the associated accretion disk, and therefore we must interpret expressions such as Equation (25) as approximations, rather than strict quantitative relations. However, these expressions are nonetheless valuable in assessing the overall validity of our model and the related parameters, which we will discuss in more detail in Section 7.3.

Quantum mechanical effects play an important role in the strong magnetic fields (B ∼ 10¹² G) inherent to X-ray pulsars because the cyclotron energy, _cyc, separating the ground state from the first excited Landau level,

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{cyc}}=\displaystyle \frac{{eBh}}{2\pi {m}_{e}c}\approx 11.57\left(\displaystyle \frac{B}{{10}^{12}\,{\rm{G}}}\right)\,\mathrm{keV},\end{eqnarray} \tag{ 26 }$$

is in the range of _cyc ∼ 10–50 keV, where m_e, e, h, and c denote the electron mass and charge, Planck's constant, and the speed of light, respectively. The resulting cyclotron absorption feature can be clearly identified in many X-ray pulsar spectra (e.g., White et al. 1983).

The strong magnetic field inside the accretion column differentiates the photons into ordinary and extraordinary polarization modes. In the case of the ordinary mode, the electric field vector is oriented in the plane formed by the magnetic field and the photon propagation direction. In the case of the extraordinary mode, the electric field vector is aligned perpendicular to this plane. The details of the photon–electron scattering process depend on the relationship between the photon energy and the cyclotron energy _cyc, and also on the propagation direction and polarization state of the photon (Chanan et al. 1979; Ventura 1979; Nagel 1980).

In the ordinary polarization mode (m = 1), the scattering cross-section is given by

$$\begin{eqnarray}&&{\sigma }_{s}^{m=1}={\sigma }_{{\rm{T}}}[{\sin }^{2}{\theta }_{s}+{f}_{s}(\epsilon ){\cos }^{2}{\theta }_{s}],\end{eqnarray} \tag{ 27 }$$

and the extraordinary mode (m = 2) scattering cross-section can be written as

$$\begin{eqnarray}&&{\sigma }_{s}^{m=2}={\sigma }_{{\rm{T}}}{f}_{s}(\epsilon )+{\rm{\Psi }},\end{eqnarray} \tag{ 28 }$$

where σ_T is the Thomson cross-section, θ_s is the angle between the photon propagation direction and the magnetic field, Ψ is the resonant contribution, and the function f_s() is defined in terms of the cyclotron energy, _cyc, by

$$\begin{eqnarray}{f}_{s}(\epsilon )\equiv \left\{\begin{array}{ll}1, & \epsilon \geqslant {\epsilon }_{\mathrm{cyc}},\\ {(\epsilon /{\epsilon }_{\mathrm{cyc}})}^{2}, & \epsilon \lt {\epsilon }_{\mathrm{cyc}}.\end{array}\right.\end{eqnarray} \tag{ 29 }$$

A complete treatment of the energy and angular dependence of the scattering of the ordinary and extraordinary mode photons is beyond the scope of this paper, and therefore we follow Wang & Frank (1981) and BW07 by splitting the photons into two populations: those propagating either parallel or perpendicular to the magnetic field direction.

Photons propagating perpendicular to the magnetic field (θ_s = π/2) are dominated by the ordinary polarization mode (m = 1) if is below the cyclotron energy, _cyc, because in this case the resonant portion of Equation (28) makes no contribution, and we find that ${\sigma }_{s}^{m=2}\lt {\sigma }_{s}^{m=1}={\sigma }_{{\rm{T}}}$. In this situation, we can therefore set the perpendicular scattering cross-section equal to the Thomson value (Ventura 1979; Becker 1998),

$$\begin{eqnarray}&&{\sigma }_{\perp }={\sigma }_{{\rm{T}}}.\end{eqnarray} \tag{ 30 }$$

For photons propagating parallel to the magnetic field (θ_s = 0), with energy  < _cyc, both modes see the Thomson cross-section reduced by the ratio ${(\epsilon /{\epsilon }_{\mathrm{cyc}})}^{2}$. In this case, we follow Arons et al. (1987) and remove the energy dependence of the parallel scattering cross-section by replacing with the radius-dependent mean photon energy, $\bar{\epsilon }(r)$, so that σ_∥(r) varies as

$$\begin{eqnarray}{\sigma }_{\parallel }\approx \left\{\begin{array}{ll}{\sigma }_{{\rm{T}}}, & \bar{\epsilon }\geqslant {\epsilon }_{\mathrm{cyc}},\\ {\sigma }_{{\rm{T}}}{(\bar{\epsilon }/{\epsilon }_{\mathrm{cyc}})}^{2}, & \bar{\epsilon }\lt {\epsilon }_{\mathrm{cyc}}.\end{array}\right.\end{eqnarray} \tag{ 31 }$$

In our computational approach, the value of σ_∥ is obtained as part of an iterative parameter variation procedure in which we self-consistently compute the radiation spectrum and the hydrodynamic structure of the accretion column, and attempt to fit the observational spectral data with adherence to the appropriate boundary conditions (see Section 4). However, as a check on the validity of the model parameters, we will refer to Equation (31) in our discussion in Section 7 in order to verify that the resulting values for σ_∥ are physically reasonable. We also require that the angle-averaged cross-section, $\overline{\sigma }$, used in the solution of the photon transport equation, must satisfy the constraint ${\sigma }_{\parallel }\lt \overline{\sigma }\lt {\sigma }_{\perp }$ (Canuto et al. 1971; BW07).

The magnetic field near the surface of an accreting neutron star is so strong that the cyclotron energy given by Equation (26) becomes comparable to the thermal energy of the electrons. Consequently, the electron energy distribution is quantized in the direction perpendicular to the magnetic field, and therefore the electrons possess a one-dimensional Maxwellian distribution along the magnetic field direction, with a mean thermal energy equal to (1/2) kT_e, where k is Boltzmann's constant. On the other hand, the proton energy is not quantized, and therefore the protons are described by a three-dimensional Maxwellian distribution, with a mean thermal energy equal to (3/2) kT_i. The ion and electron internal energy densities are therefore given by

$$\begin{eqnarray}&&{U}_{i}=\displaystyle \frac{3}{2}\,{n}_{i}{{kT}}_{i},\ \ \ {U}_{e}=\displaystyle \frac{1}{2}\,{n}_{e}{{kT}}_{e},\end{eqnarray} \tag{ 32 }$$

where n_i and n_e denote the ion and electron number densities, respectively. In principle, the ion and electron temperatures T_i and T_e are not necessarily equal, and therefore in our two-temperature model we implement separate energy equations for each species, including a term describing their Coulomb coupling.

The magnetic field pressure is orders of magnitude stronger than either the gas pressure or the radiation pressure in an X-ray pulsar accretion column, and therefore the charged particles are constrained to follow the curved dipole magnetic field as the plasma flows downward toward the stellar surface. Charge neutrality ensures that n_i = n_e at all locations. From the point of view of the accretion hydrodynamics, the relevant pressure is the total pressure parallel to the local magnetic field direction, given by the sum of the electron, ion, and radiation components,

$$\begin{eqnarray}&&{P}_{\mathrm{tot}}={P}_{e}+{P}_{i}+{P}_{r},\end{eqnarray} \tag{ 33 }$$

where

$$\begin{eqnarray}&&{P}_{i}={n}_{i}{{kT}}_{i},\qquad {P}_{e}={n}_{e}{{kT}}_{e},\end{eqnarray} \tag{ 34 }$$

denote the ion and electron pressures, respectively. The radiation pressure, P_r, is not given by a thermal formula since the X-ray pulsar radiation field is nonthermal. Hence the radiation pressure must be computed using a conservation relation. The pressure components are related to their corresponding energy densities via

$$\begin{eqnarray}&&{P}_{i}=({\gamma }_{i}-1)\,{U}_{i},\quad {P}_{e}=({\gamma }_{e}-1)\,{U}_{e},\quad {P}_{r}=({\gamma }_{r}-1)\,{U}_{r},\end{eqnarray} \tag{ 35 }$$

where it follows from Equations (32), (34), and (35) that γ_e = 3 and γ_i = 5/3. The ratio of specific heats for the radiation is γ_r = 4/3.

3.0.1. Mass Flux

In the one-dimensional case considered here, the cross-sectional structure of the accretion column is not considered in detail, and all of the densities and temperatures represent averages across the column at a given radius r. Hence the mass continuity equation can be written in dipole geometry as (e.g., Langer & Rappaport 1982)

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho (r)}{\partial t}=-\displaystyle \frac{1}{A(r)}\displaystyle \frac{\partial }{\partial r}[A(r)\rho (r)v(r)],\end{eqnarray} \tag{ 38 }$$

where v < 0 denotes the radial inflow velocity, and the cross-sectional area of the column, A(r), is related to the solid angle, Ω(r) = (r/R_*)Ω_*, by

$$\begin{eqnarray}&&A(r)={r}^{2}{\rm{\Omega }}(r)=\displaystyle \frac{{r}^{3}{{\rm{\Omega }}}_{* }}{{R}_{* }}.\end{eqnarray} \tag{ 39 }$$

In a steady state, we see from Equation (38) that A(r)ρ(r)v(r) is a conserved quantity, i.e., the mass accretion rate $\dot{M}$ is conserved and is related to the density ρ and velocity v via

$$\begin{eqnarray}&&\dot{M}={\rm{\Omega }}\,{r}^{2}\rho | v| ,\end{eqnarray} \tag{ 40 }$$

which can be combined with Equation (8) to obtain for the mass density

$$\begin{eqnarray}&&\rho =\displaystyle \frac{\dot{M}{R}_{* }}{{{\rm{\Omega }}}_{* }{r}^{3}| v| }.\end{eqnarray} \tag{ 41 }$$

This algebraic relation for the density is used to supplement the set of differential conservation equations in our hydrodynamical model for the column structure.

We assume that the accreting gas is composed of pure, fully ionized hydrogen, and therefore the electron and ion number densities are given by

$$\begin{eqnarray}&&{n}_{e}={n}_{i}=\displaystyle \frac{\rho }{{m}_{\mathrm{tot}}}=\displaystyle \frac{\dot{M}{R}_{* }}{{m}_{\mathrm{tot}}{{\rm{\Omega }}}_{* }{r}^{3}| v| },\end{eqnarray} \tag{ 42 }$$

where ${m}_{\mathrm{tot}}={m}_{e}+{m}_{i}$.

3.0.2. Total Energy Flux

The total energy flux in the radial direction, averaged over the column cross-section at radius r, is given by

$$\begin{eqnarray}F(r) & = & \displaystyle \frac{1}{2}\rho {v}^{3}+v({P}_{i}+{P}_{e}+{P}_{r}+{U}_{i}+{U}_{e}+{U}_{r})\\ & & -\displaystyle \frac{c}{{n}_{e}{\sigma }_{\parallel }}\displaystyle \frac{\partial {P}_{r}}{\partial r}-\displaystyle \frac{{{GM}}_{* }\rho v}{r},\end{eqnarray} \tag{ 43 }$$

where the energy flux is defined to be negative for energy flow in the downward direction, and the accretion velocity v is negative (v < 0). The terms on the right-hand side of Equation (43) represent the kinetic energy flux, the enthalpy flux, the radiation diffusion flux, and the gravitational energy flux, respectively. The total energy flux F is related to the total energy transport rate in the radial direction, denoted by $\dot{E}$, via

$$\begin{eqnarray}&&\dot{E}(r)=A(r)\,F(r)\ \ \propto \ \ \mathrm{erg}\ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 44 }$$

where the column cross-sectional area A(r) is given by Equation (39).

We can derive a first-order differential equation for the radiation sound speed, a_r, by substituting for the energy densities and pressures in Equation (43) using Equations (35) and (36), substituting for the electron number density n_e using Equation (42), and substituting for F using Equation (44). After some algebra, we obtain in the steady-state case

$$\begin{eqnarray}\displaystyle \frac{{{da}}_{r}}{{dr}} & = & \displaystyle \frac{3{a}_{r}}{2r}+\displaystyle \frac{{a}_{r}}{2v}\displaystyle \frac{{dv}}{{dr}}-\displaystyle \frac{1}{2}\displaystyle \frac{{\sigma }_{\parallel }{\gamma }_{r}\dot{M}}{{m}_{\mathrm{tot}}{{ca}}_{r}{\rm{\Omega }}{r}^{2}}\left(\displaystyle \frac{\dot{E}}{\dot{M}}+\displaystyle \frac{{v}^{2}}{2}\right.\\ & & +\left.\displaystyle \frac{{a}_{i}^{2}}{{\gamma }_{i}-1}+\displaystyle \frac{{a}_{e}^{2}}{{\gamma }_{e}-1}+\displaystyle \frac{{a}_{r}^{2}}{{\gamma }_{r}-1}-\displaystyle \frac{{{GM}}_{* }}{r}\right).\end{eqnarray} \tag{ 45 }$$

3.0.3. Ion and Electron Energy Equations

The variation of the internal energy density of the ionized gas is influenced by adiabatic heating, energy exchange between the ions and electrons, and the emission and absorption of radiation. Averaging over the cross-section of the column at radius r, the energy equations for the ions and electrons can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{{DU}}_{i}}{{Dt}}={\gamma }_{i}\displaystyle \frac{{U}_{i}}{\rho }\displaystyle \frac{D\rho }{{Dt}}+{\dot{U}}_{i},\ \ \ \ \ \ \displaystyle \frac{{{DU}}_{e}}{{Dt}}={\gamma }_{e}\displaystyle \frac{{U}_{e}}{\rho }\displaystyle \frac{D\rho }{{Dt}}+{\dot{U}}_{e},\end{eqnarray} \tag{ 46 }$$

respectively, where the first terms on the right-hand side represent adiabatic compression, the final terms represent thermal coupling with the other species, and the comoving (Lagrangian) time derivative D/Dt is defined by

$$\begin{eqnarray}&&\displaystyle \frac{D}{{Dt}}\equiv \displaystyle \frac{\partial }{\partial t}+v\displaystyle \frac{\partial }{\partial r}.\end{eqnarray} \tag{ 47 }$$

The thermal coupling terms appearing in Equation (46) represent the net heating due to a variety of combined processes, which are broken down as follows,

$$\begin{eqnarray}{\dot{U}}_{e} & = & {\dot{U}}_{\mathrm{brem}}^{\mathrm{emit}}+{\dot{U}}_{\mathrm{brem}}^{\mathrm{abs}}+{\dot{U}}_{\mathrm{cyc}}^{\mathrm{emit}}+{\dot{U}}_{\mathrm{cyc}}^{\mathrm{abs}}+{\dot{U}}_{\mathrm{Comp}}+{\dot{U}}_{\mathrm{ei}},\\ {\dot{U}}_{i} & = & -\,{\dot{U}}_{\mathrm{ei}}.\end{eqnarray} \tag{ 48 }$$

The terms in the expression for ${\dot{U}}_{e}$ denote, respectively, bremsstrahlung (free–free) emission and absorption, cyclotron emission and absorption, photon–electron Comptonization, and electron–ion Coulomb energy exchange. The ions do not radiate appreciably, and therefore they only experience adiabatic compression and Coulomb energy exchange (see Langer & Rappaport 1982). In our sign convention, a heating term is positive and a cooling term is negative. These energy transfer rates are discussed in more detail in Section 3.1.

In a steady state, Equation (46) can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{{dU}}_{i}}{{dr}}={\gamma }_{i}\displaystyle \frac{{U}_{i}}{\rho }\displaystyle \frac{d\rho }{{dr}}+\displaystyle \frac{{\dot{U}}_{i}}{v},\ \ \ \ \ \ \displaystyle \frac{{{dU}}_{e}}{{dr}}={\gamma }_{e}\displaystyle \frac{{U}_{e}}{\rho }\displaystyle \frac{d\rho }{{dr}}+\displaystyle \frac{{\dot{U}}_{e}}{v}.\end{eqnarray} \tag{ 49 }$$

We can derive equivalent differential equations satisfied by the electron and ion sound speeds by using Equations (35), (36), and (41) to substitute for the energy and mass densities in Equation (49), obtaining

$$\begin{eqnarray}&&\displaystyle \frac{{{da}}_{i}}{{dr}}=-({\gamma }_{i}-1)\left(\displaystyle \frac{{a}_{i}}{2v}\displaystyle \frac{{dv}}{{dr}}+\displaystyle \frac{3{a}_{i}}{2r}\right)-\displaystyle \frac{1}{2}{\gamma }_{i}({\gamma }_{i}-1)\displaystyle \frac{{\rm{\Omega }}{r}^{2}}{\dot{M}}\displaystyle \frac{{\dot{U}}_{i}}{{a}_{i}},\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{da}}_{e}}{{dr}}=-({\gamma }_{e}-1)\left(\displaystyle \frac{{a}_{e}}{2v}\displaystyle \frac{{dv}}{{dr}}+\displaystyle \frac{3{a}_{e}}{2r}\right)-\displaystyle \frac{1}{2}{\gamma }_{e}({\gamma }_{e}-1)\displaystyle \frac{{\rm{\Omega }}{r}^{2}}{\dot{M}}\displaystyle \frac{{\dot{U}}_{e}}{{a}_{e}}.\end{eqnarray} \tag{ 51 }$$

3.0.4. Momentum Equation

The ionized, accreting gas is constrained to spiral around the magnetic field lines by the Lorentz force. Since there is no component of the Lorentz force parallel to the local B-field, the remaining acceleration in the parallel direction is due to the total pressure gradient and the gravitational field of the neutron star. If we average over the cross-section of the accretion column at radius r, then the comoving acceleration in the radial direction can be written as (e.g., Langer & Rappaport 1982),

$$\begin{eqnarray}&&\displaystyle \frac{{Dv}}{{Dt}}=-\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial {P}_{\mathrm{tot}}}{\partial r}-\displaystyle \frac{{{GM}}_{* }}{{r}^{2}},\end{eqnarray} \tag{ 52 }$$

where ${P}_{\mathrm{tot}}={P}_{r}+{P}_{i}+{P}_{e}$ is the total pressure, and the Lagrangian time derivative D/Dt is defined by Equation (47). Substituting for the mass density ρ and the pressure components P_i, P_e, and P_r using Equations (41) and (36), respectively, we can derive a first-order differential equation satisfied by the fluid velocity v involving the sound speeds a_i, a_e, and a_r, and the energy transport rate $\dot{E}$. After some algebra, the result obtained in a steady state is

$$\begin{eqnarray}\displaystyle \frac{{dv}}{{dr}} & = & \displaystyle \frac{v}{{v}^{2}-({a}_{i}^{2}+{a}_{e}^{2})}\left\{\displaystyle \frac{3({a}_{i}^{2}+{a}_{e}^{2})}{r}\right.\\ & & -\displaystyle \frac{{{GM}}_{* }}{{r}^{2}}+\displaystyle \frac{{\sigma }_{\parallel }\dot{M}}{{m}_{\mathrm{tot}}c\,{\rm{\Omega }}{r}^{2}}\left(\displaystyle \frac{\dot{E}}{\dot{M}}+\displaystyle \frac{{v}^{2}}{2}\right.\\ & & +\left.\displaystyle \frac{{a}_{i}^{2}}{{\gamma }_{i}-1}+\displaystyle \frac{{a}_{e}^{2}}{{\gamma }_{e}-1}+\displaystyle \frac{{a}_{r}^{2}}{{\gamma }_{r}-1}-\displaystyle \frac{{{GM}}_{* }}{r}\right)\\ & & +\left.\displaystyle \frac{{\rm{\Omega }}{r}^{2}}{\dot{M}}[({\gamma }_{i}-1){\dot{U}}_{i}+({\gamma }_{e}-1){\dot{U}}_{e}]\right\},\end{eqnarray} \tag{ 53 }$$

where we have also made use of Equations (45), (50), and (51).

3.0.5. Radiative Losses

The value of the energy transport rate $\dot{E}$ (Equation (44)) varies as a function of the radius r in response to the escape of radiation energy through the walls of the accretion column, perpendicular to the magnetic field direction. In our one-dimensional model, all quantities are averaged over the cross-section of the column, and therefore we use an escape-probability formalism to account for the diffusion of radiation through the walls of the column. We therefore utilize a total energy conservation equation of the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{1}{2}\,\rho {v}^{2}+{U}_{i}+{U}_{e}+{U}_{r}-\displaystyle \frac{{{GM}}_{* }\rho }{r}\right)\\ &&\quad =\,-\,\displaystyle \frac{1}{A(r)}\displaystyle \frac{\partial }{\partial r}[A(r)F(r)]+{\dot{U}}_{\mathrm{esc}},\end{eqnarray} \tag{ 54 }$$

where the total energy flux is given by $F(r)=\dot{E}(r)/A(r)$ (see Equation (44)), and the energy escape rate per unit volume is given by

$$\begin{eqnarray}&&{\dot{U}}_{\mathrm{esc}}=-\displaystyle \frac{{U}_{r}}{{t}_{\mathrm{esc}}},\ \ \ \ {t}_{\mathrm{esc}}=\displaystyle \frac{{{\ell }}_{\mathrm{esc}}}{{w}_{\perp }}.\end{eqnarray} \tag{ 55 }$$

Here, t_esc(r) represents the mean escape time for photons to diffuse across the column and escape through the walls, w_⊥(r) is the perpendicular diffusion velocity, and ℓ_esc(r) denotes the perpendicular escape distance across the column at radius r, computed using (see Figure 2)

$$\begin{eqnarray}&&{{\ell }}_{\mathrm{esc}}(r)=({{\ell }}_{2}-{{\ell }}_{1}){\left(\displaystyle \frac{r}{{R}_{* }}\right)}^{3/2},\end{eqnarray} \tag{ 56 }$$

so that at the stellar surface, we obtain ${{\ell }}_{\mathrm{esc}}={{\ell }}_{2}-{{\ell }}_{1}$, as required. The perpendicular diffusion velocity w_⊥ cannot exceed the speed of light, and therefore we compute it using the constrained formula

$$\begin{eqnarray}&&{w}_{\perp }=\min \left(c,\displaystyle \frac{c}{{\tau }_{\perp }}\right),\ \ \ \ {\tau }_{\perp }={n}_{e}{\sigma }_{\perp }{{\ell }}_{\mathrm{esc}},\end{eqnarray} \tag{ 57 }$$

where τ_⊥ denotes the perpendicular optical thickness of the column at radius r.

In a steady state, Equation (54) reduces to

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}{\rm{\Omega }}(r)}\displaystyle \frac{d\dot{E}}{{dr}}=-\displaystyle \frac{{U}_{r}}{{{\ell }}_{\mathrm{esc}}}\,\min \left(c,\displaystyle \frac{c}{{\tau }_{\perp }}\right).\end{eqnarray} \tag{ 58 }$$

By combining Equations (35), (36), (40), and (58), we can obtain the final form for the energy transport differential equation,

$$\begin{eqnarray}&&\displaystyle \frac{d\dot{E}}{{dr}}=-\displaystyle \frac{{a}_{r}^{2}\dot{M}}{{\gamma }_{r}({\gamma }_{r}-1)\,{{\ell }}_{\mathrm{esc}}\,| v| }\,\min \left(c,\displaystyle \frac{c}{{\tau }_{\perp }}\right).\end{eqnarray} \tag{ 59 }$$

The energy exchange rates per unit volume introduced in Equations (48), denoted by ${\dot{U}}_{\mathrm{brem}}^{\mathrm{emit}}$, ${\dot{U}}_{\mathrm{brem}}^{\mathrm{abs}}$, ${\dot{U}}_{\mathrm{cyc}}^{\mathrm{emit}}$, ${\dot{U}}_{\mathrm{cyc}}^{\mathrm{abs}}$, ${\dot{U}}_{\mathrm{Comp}}$, and ${\dot{U}}_{\mathrm{ei}}$, describe a comprehensive set of heating and cooling processes experienced by the gas and radiation, including Coulomb coupling between the ions and electrons, the Compton exchange of energy between the electrons and photons, and the emission and absorption of radiation energy via thermal bremsstrahlung and cyclotron. In this section, we provide additional details regarding the computation of these various rates.

3.1.1. Bremsstrahlung Emission and Absorption

Thermal bremsstrahlung emission plays a significant role in cooling the ionized gas, and in the case of luminous X-ray pulsars, it also provides the majority of the seed photons that are subsequently Compton scattered to form the emergent X-ray spectrum (BW07). Assuming a fully ionized hydrogen composition for the accreting gas, with n_e = n_i, the total power per unit volume emitted by the electrons is given by (see Rybicki & Lightman 1979, Equation (5.14)),

$$\begin{eqnarray}&&{\dot{U}}_{\mathrm{brem}}^{\mathrm{emit}}=-{\left(\displaystyle \frac{2\pi {{kT}}_{e}}{3{m}_{e}}\right)}^{1/2}\displaystyle \frac{{2}^{5}\pi {e}^{6}}{3{{hm}}_{e}{c}^{3}}\,{n}_{e}^{2},\end{eqnarray} \tag{ 60 }$$

where we have set the Gaunt factor equal to unity. The negative sign appears in Equation (60) because this term represents a cooling process in which heat is removed from the electrons. We can write an equivalent expression for the bremsstrahlung cooling rate in terms of the electron sound speed, a_e, by using Equation (37) to eliminate the electron temperature T_e in Equation (60), thereby obtaining, in cgs units,

$$\begin{eqnarray}&&{\dot{U}}_{\mathrm{brem}}^{\mathrm{emit}}=-3.2\times {10}^{16}\,{a}_{e}\,{\rho }^{2},\end{eqnarray} \tag{ 61 }$$

where we have also used Equation (42).

The electrons in the accretion column also experience heating due to free–free absorption of low-frequency radiation, which plays an important role in regulating the temperature of the gas. The heating rate per unit volume due to thermal bremsstrahlung absorption, integrated over photon frequency, is given by

$$\begin{eqnarray}&&{\dot{U}}_{\mathrm{brem}}^{\mathrm{abs}}={U}_{r}\,{\alpha }_{{\rm{R}}}\,c,\end{eqnarray} \tag{ 62 }$$

where α_R is the Rosseland mean absorption coefficient for fully ionized hydrogen, expressed in cgs units by (Rybicki & Lightman 1979)

$$\begin{eqnarray}&&{\alpha }_{{\rm{R}}}=1.7\times {10}^{-25}\,{T}_{e}^{-7/2}{n}_{e}^{2}\propto \,{\mathrm{cm}}^{-1}.\end{eqnarray} \tag{ 63 }$$

Note that we have set the Gaunt factor equal to unity and assumed that the gas is composed of fully ionized hydrogen. By combining Equations (62) and (63) and substituting for ${U}_{r}={P}_{r}/({\gamma }_{r}-1)$, n_e, and T_e, using Equations (36), (37), and (42), respectively, we obtain, in cgs units

$$\begin{eqnarray}&&{\dot{U}}_{\mathrm{brem}}^{\mathrm{abs}}=9.8\times {10}^{62}\,{a}_{e}^{-7}{a}_{r}^{2}\,{\rho }^{3}.\end{eqnarray} \tag{ 64 }$$

The sign of this quantity is positive since it represents a heating process for the electrons.

3.1.2. Cyclotron Emission and Absorption

The electrons in the accretion column also experience heating and cooling due to the emission and absorption of thermal cyclotron radiation. At any given time, most of the electrons are found in the ground state, but they can be excited to the first Landau level via collisions, or via the absorption of radiation at the cyclotron energy, _cyc. At the densities and temperatures prevalent in pulsar accretion columns, radiative excitation is followed immediately by radiative de-excitation back to the ground state, so that in net terms, cyclotron absorption can be interpreted as a resonant scattering process, which results in no net change in the angle-averaged photon distribution (Nagel 1980; Arons et al. 1987). Hence, on average, cyclotron absorption does not result in the net heating of the gas, due to the rapid radiative de-excitation, and we therefore set ${\dot{U}}_{\mathrm{cyc}}^{\mathrm{abs}}=0$ in our dynamical calculations. However, near the surface of the accretion column, photons scattered out of the outwardly directed beam are not replaced, and this leads to the formation of the observed cyclotron absorption features, in a process that is very analogous to the formation of absorption lines in the solar spectrum (Ventura et al. 1979). The formation of the cyclotron absorption features is further considered in Paper II.

While cyclotron absorption does not result in the net heating of the gas, due to the rapid radiative de-excitation, cyclotron emission will cool the gas. In this process, kinetic energy is converted into excitation energy via collisions, and the subsequent emission of cyclotron radiation removes heat from the electrons. To compute the cyclotron cooling rate, ${\dot{U}}_{\mathrm{cyc}}^{\mathrm{emit}}$, we begin with the cyclotron emissivity, ${\dot{n}}_{\epsilon }^{\mathrm{cyc}}$, which gives the production rate of cyclotron photons per unit volume per unit energy. Using Equations (7) and (11) from Arons et al. (1987), we have

$$\begin{eqnarray}&&{\dot{n}}_{\epsilon }^{\mathrm{cyc}}=2.1\times {10}^{36}{\rho }^{2}{B}_{12}^{-3/2}H\left(\displaystyle \frac{{\epsilon }_{{cyc}}}{{{kT}}_{e}}\right){e}^{-{\epsilon }_{\mathrm{cyc}}/{{kT}}_{e}}\delta (\epsilon -{\epsilon }_{\mathrm{cyc}}),\end{eqnarray} \tag{ 65 }$$

where B₁₂ = $| {\boldsymbol{B}}|$/(10¹² G) is evaluated using Equation (4) with θ = 0, and H(x) is a piecewise function defined by

$$\begin{eqnarray}H(x)\equiv \left\{\begin{array}{ll}0.15\,\sqrt{7.5}, & x\geqslant 7.5,\\ 0.15\,\sqrt{x}, & x\lt 7.5.\end{array}\right.\end{eqnarray} \tag{ 66 }$$

The total cyclotron cooling rate is obtained by multiplying Equation (65) by the photon energy and integrating over all energies, which yields, in cgs units,

$$\begin{eqnarray}&&{\dot{U}}_{\mathrm{cyc}}^{\mathrm{emit}}=-2.1\times {10}^{36}{\rho }^{2}{B}_{12}^{-3/2}{\epsilon }_{\mathrm{cyc}}\ H\left(\displaystyle \frac{{\epsilon }_{\mathrm{cyc}}}{{{kT}}_{e}}\right){e}^{-{\epsilon }_{\mathrm{cyc}}/{{kT}}_{e}},\end{eqnarray} \tag{ 67 }$$

where the negative sign indicates this is a cooling process for the electrons.

3.1.3. Compton Heating and Cooling

Compton scattering plays a fundamental role in the formation of the emergent X-ray spectrum. It is also critically important in establishing the radial variation of the electron temperature profile through the exchange of energy between the photons and electrons. Equation (7.36) from Rybicki & Lightman (1979) gives the mean change in the photon energy during a single scattering as

$$\begin{eqnarray}&&\langle {\rm{\Delta }}\epsilon \rangle =\displaystyle \frac{\epsilon }{{m}_{e}{c}^{2}}(4{{kT}}_{e}-\epsilon ),\end{eqnarray} \tag{ 68 }$$

and the associated mean rate of change of the photon energy is therefore

$$\begin{eqnarray}&&{\left.\left\langle \displaystyle \frac{d\epsilon }{{dt}}\right\rangle \right|}_{\mathrm{Comp}}={n}_{e}\bar{\sigma }c\langle {\rm{\Delta }}\epsilon \rangle ,\end{eqnarray} \tag{ 69 }$$

where ${({n}_{e}\bar{\sigma }c)}^{-1}$ denotes the mean-free time between scatterings for the photons, and $\bar{\sigma }$ is the angle-averaged electron scattering cross-section (BW07). The corresponding rate of change of the electron energy density due to Compton scattering can therefore be written as

$$\begin{eqnarray}&&{\dot{U}}_{\mathrm{Comp}}=-{n}_{e}\bar{\sigma }c{\int }_{0}^{\infty }{\epsilon }^{2}f(r,\epsilon )\langle {\rm{\Delta }}\epsilon \rangle \,d\epsilon ,\end{eqnarray} \tag{ 70 }$$

where the distribution function, f(r, ), is the solution to the photon transport equation introduced in Paper II, which is related to the total radiation number density, n_r, and energy density, U_r, via

$$\begin{eqnarray}&&{n}_{r}(r)={\int }_{0}^{\infty }{\epsilon }^{2}f(r,\epsilon )\,d\epsilon ,\ \ \ \ {U}_{r}(r)={\int }_{0}^{\infty }{\epsilon }^{3}f(r,\epsilon )\,d\epsilon .\end{eqnarray} \tag{ 71 }$$

Combining Equations (68) and (70), we find that the net Compton cooling rate for the electrons is given by

$$\begin{eqnarray}&&{\dot{U}}_{\mathrm{Comp}}=\displaystyle \frac{{n}_{e}\bar{\sigma }c}{{m}_{e}{c}^{2}}\left[{\int }_{0}^{\infty }{\epsilon }^{4}\,f(r,\epsilon )d\epsilon -4{{kT}}_{e}{\int }_{0}^{\infty }{\epsilon }^{3}f(r,\epsilon )\,d\epsilon \right],\end{eqnarray} \tag{ 72 }$$

which vanishes if the electron temperature, T_e, is equal to the inverse-Compton temperature, T_IC, defined by

$$\begin{eqnarray}&&{T}_{\mathrm{IC}}(r)\equiv \displaystyle \frac{1}{4k}\displaystyle \frac{{\int }_{0}^{\infty }{\epsilon }^{4}\,f(r,\epsilon )\,d\epsilon }{{\int }_{0}^{\infty }{\epsilon }^{3}\,f(r,\epsilon )d\epsilon }.\end{eqnarray} \tag{ 73 }$$

In the present paper, we are primarily interested in the implications of Compton scattering for the heating and cooling of the gas, and its effect on the dynamical structure of the accretion column. The electron cooling rate can be rewritten as

$$\begin{eqnarray}&&{\dot{U}}_{\mathrm{Comp}}={n}_{e}\bar{\sigma }c\ \displaystyle \frac{4{{kT}}_{e}}{{m}_{e}{c}^{2}}[g(r)-1]{U}_{r},\end{eqnarray} \tag{ 74 }$$

where we introduce g(r) as the temperature ratio function,

$$\begin{eqnarray}&&g(r)\equiv \displaystyle \frac{{T}_{\mathrm{IC}}}{{T}_{e}}.\end{eqnarray} \tag{ 75 }$$

The sign of ${\dot{U}}_{\mathrm{Comp}}$ depends on the value of g. If g < 1 (i.e., T_IC < T_e), then the electrons experience Compton cooling; otherwise, the electrons are heated via inverse-Compton scattering. We can obtain the final form for the Compton cooling rate in terms of the mass density, ρ, the electron sound speed, a_e, and the radiation sound speed, a_r by combining Equations (34)–(36), and (74), which yields

$$\begin{eqnarray}&&{\dot{U}}_{\mathrm{Comp}}=\displaystyle \frac{4\,\bar{\sigma }\,(g-1)}{{m}_{e}c{\gamma }_{e}{\gamma }_{r}({\gamma }_{r}-1)}\ {\rho }^{2}\,{a}_{r}^{2}\,{a}_{e}^{2}.\end{eqnarray} \tag{ 76 }$$

3.1.4. Electron–Ion Energy Exchange

The electrons can also be heated or cooled via Coulomb collisions with the protons, depending on whether the electron temperature T_e exceeds the ion temperature T_i. The net heating rate per unit volume for the electrons is given by (Langer & Rappaport 1982)

$$\begin{eqnarray}{\dot{U}}_{\mathrm{ei}} & = & \displaystyle \frac{3}{2}{\left(\displaystyle \frac{2}{\pi }\right)}^{1/2}{\sigma }_{{\rm{T}}}{c}^{3}{m}_{e}{n}_{e}^{2}\left(\displaystyle \frac{{m}_{e}}{{m}_{i}}\right)\left(\displaystyle \frac{{T}_{i}-{T}_{e}}{{T}_{\mathrm{eff}}}\right)\\ & & \times \,{\left(\displaystyle \frac{{m}_{e}{c}^{2}}{{{kT}}_{\mathrm{eff}}}\right)}^{1/2}\mathrm{ln}\ {{\rm{\Lambda }}}_{\mathrm{Coul}},\end{eqnarray} \tag{ 77 }$$

where

$$\begin{eqnarray}&&{T}_{\mathrm{eff}}\equiv {T}_{e}+\left(\displaystyle \frac{{m}_{e}}{{m}_{i}}\right){T}_{i},\end{eqnarray} \tag{ 78 }$$

and the Coulomb logarithm is given by

$$\begin{eqnarray}&&\mathrm{ln}\ {{\rm{\Lambda }}}_{\mathrm{Coul}}=5.41+\displaystyle \frac{1}{4}\mathrm{ln}\ \left(\displaystyle \frac{{{kT}}_{\mathrm{eff}}}{20\,\mathrm{keV}}\displaystyle \frac{B}{{10}^{12}\,{\rm{G}}}\displaystyle \frac{{10}^{20}\,{\mathrm{cm}}^{-3}}{{n}_{e}}\right).\end{eqnarray} \tag{ 79 }$$

We can further simplify Equation (77) by substituting for n_e using Equation (42) and substituting for T_e and T_i using Equation (37), obtaining

$$\begin{eqnarray}{\dot{U}}_{\mathrm{ei}} & = & \displaystyle \frac{3}{2}{\left(\displaystyle \frac{2}{\pi }\right)}^{1/2}{\sigma }_{{\rm{T}}}{c}^{4}\left(\displaystyle \frac{{m}_{e}}{{m}_{i}}\right){\left(1+\displaystyle \frac{{m}_{i}}{{m}_{e}}\right)}^{-5/2}\\ & & \times \,\left(\displaystyle \frac{{a}_{i}^{2}}{{\gamma }_{i}}-\displaystyle \frac{{a}_{e}^{2}}{{\gamma }_{e}}\right){\left(\displaystyle \frac{{a}_{e}^{2}}{{\gamma }_{e}{m}_{e}}+\displaystyle \frac{{a}_{i}^{2}}{{\gamma }_{i}{m}_{i}}\right)}^{-3/2}{\rho }^{2}\ \mathrm{ln}\ {{\rm{\Lambda }}}_{\mathrm{Coul}},\end{eqnarray} \tag{ 80 }$$

which in cgs units becomes

$$\begin{eqnarray}{\dot{U}}_{\mathrm{ei}} & = & 2.4\times {10}^{6}\left(\displaystyle \frac{{a}_{i}^{2}}{{\gamma }_{i}}-\displaystyle \frac{{a}_{e}^{2}}{{\gamma }_{e}}\right){\left(\displaystyle \frac{{a}_{e}^{2}}{{\gamma }_{e}{m}_{e}}+\displaystyle \frac{{a}_{i}^{2}}{{\gamma }_{i}{m}_{i}}\right)}^{-3/2}\\ & & \times \,{\rho }^{2}\ \mathrm{ln}\ {{\rm{\Lambda }}}_{\mathrm{Coul}}.\end{eqnarray} \tag{ 81 }$$

Note that when T_e = T_i, the second factor in Equation (81) is zero, and thus ${\dot{U}}_{{ei}}=0$, as expected. Based on the symmetry of the energy exchange between the particle species, we immediately conclude that the energy transfer rate per unit volume for the protons is given by ${\dot{U}}_{i}=-{\dot{U}}_{\mathrm{ei}}$ (see Equation (48)).

The upper surface of the dipole-shaped accretion funnel is located at radius r = r_top, which must be below the radius corresponding to the turnover height of the dipole field, r_c, as discussed in Section 2.2 (see Equation (19)). In analogy with the theory of stellar atmospheres, the top of the accretion column represents the last scattering surface for photon–electron interaction as photons travel out the top of the column, implying that the scattering optical depth from r_top to r_c should equal unity. Defining the parallel scattering optical depth, τ_∥, so that it increases in the downward direction for bulk fluid entering at the top of the column and flowing downward, from τ_∥ = 0 at r = r_top, we have

$$\begin{eqnarray}&&{\tau }_{\parallel }(r)={\int }_{r}^{{r}_{\mathrm{top}}}{n}_{e}(r^{\prime} ){\sigma }_{\parallel }{dr}^{\prime} .\end{eqnarray} \tag{ 82 }$$

Since the top of the accretion column is the last scattering surface, we can also write

$$\begin{eqnarray}&&{\int }_{{r}_{\mathrm{top}}}^{{r}_{c}}{n}_{e}(r^{\prime} ){\sigma }_{\parallel }{dr}^{\prime} =1,\end{eqnarray} \tag{ 83 }$$

where r_top < r_c.

We can use Equation (83) to constrain the radius at the top of the accretion column, r_top, as follows. We assume the gas is in free-fall above r_top, with velocity

$$\begin{eqnarray}&&v(r)={v}_{\mathrm{ff}}(r)\equiv -{\left(\displaystyle \frac{2{{GM}}_{* }}{r}\right)}^{1/2},\qquad r\gt {r}_{\mathrm{top}}.\end{eqnarray} \tag{ 84 }$$

Using Equation (84) to substitute for v in Equation (42) yields for the variation of the electron number density n_e the result

$$\begin{eqnarray}&&{n}_{e}(r)=\displaystyle \frac{\dot{M}{R}_{* }}{{m}_{\mathrm{tot}}{{\rm{\Omega }}}_{* }{r}^{3}}{\left(\displaystyle \frac{2{{GM}}_{* }}{r}\right)}^{-1/2},\end{eqnarray} \tag{ 85 }$$

where ${m}_{\mathrm{tot}}={m}_{e}+{m}_{i}$.

By utilizing Equation (85) to substitute for the electron number density n_e in Equation (83) and carrying out the radial integration, we obtain the condition

$$\begin{eqnarray}&&\displaystyle \frac{2{\sigma }_{\parallel }\dot{M}}{3{m}_{\mathrm{tot}}{(2{{GM}}_{* })}^{1/2}}\displaystyle \frac{{R}_{* }}{{{\rm{\Omega }}}_{* }}({r}_{\mathrm{top}}^{-3/2}-{r}_{c}^{-3/2})=1,\end{eqnarray} \tag{ 86 }$$

where the left-hand side is positive definite, since r_top < r_c, and the dipole turnover radius r_c is given by Equation (18). By rearranging Equation (86), we can obtain an explicit expression for r_top, given by

$$\begin{eqnarray}&&{r}_{\mathrm{top}}={\left[{r}_{c}^{-3/2}+\displaystyle \frac{3{m}_{\mathrm{tot}}{(2{{GM}}_{* })}^{1/2}}{2{\sigma }_{\parallel }\dot{M}}\displaystyle \frac{{{\rm{\Omega }}}_{* }}{{R}_{* }}\right]}^{-2/3}.\end{eqnarray} \tag{ 87 }$$

This relation allows us to self-consistently compute the value of r_top in terms of the parameters Ω_*, r_c, and σ_∥ in our model.

At the top of the accretion column, the inflow velocity v equals the local free-fall velocity, so that

$$\begin{eqnarray}&&{v}_{\mathrm{top}}\equiv {v}_{\mathrm{ff}}({r}_{\mathrm{top}})=-{\left(\displaystyle \frac{2{{GM}}_{* }}{{r}_{\mathrm{top}}}\right)}^{1/2}.\end{eqnarray} \tag{ 88 }$$

We also assume that at the top of the accretion column, the local acceleration of the gas is equal to the gravitational value, so that

$$\begin{eqnarray}&&{\left.v\displaystyle \frac{{dv}}{{dr}}\right|}_{r={r}_{\mathrm{top}}}=-\displaystyle \frac{{GM}}{{r}_{\mathrm{top}}^{2}},\end{eqnarray} \tag{ 89 }$$

which implies that

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{dv}}{{dr}}\right|}_{r={r}_{\mathrm{top}}}={\left(\displaystyle \frac{{{GM}}_{* }}{2{r}_{\mathrm{top}}^{3}}\right)}^{1/2}.\end{eqnarray} \tag{ 90 }$$

By assuming pure gravitational acceleration at the top of the accretion column, we are implicitly neglecting the effects of the radiation pressure gradient, which will partially counteract the downward gravitational force. We revisit this issue in Section 7, where we conclude that this assumption is warranted, since most of the radiation escapes out the sides of the accretion column as a fan beam in the high-luminosity sources of interest here. However, in lower-luminosity sources, a larger fraction of the radiation may escape out the top of the column via a pencil-beam component, but even in this case, the effect of radiation deceleration at the top of the column is still likely to be negligible.

Although our calculation allows for the possibility of two-temperature flow, with unequal values of T_i and T_e, in luminous X-ray pulsar accretion columns, significant deviation between the two temperatures is not expected, because the thermal equilibration timescale is much smaller than the dynamical timescale (BW07). We therefore assume that T_i = T_e for the inflowing gas at the top of the column (Elsner & Lamb 1977), so that

$$\begin{eqnarray}&&{T}_{i,\mathrm{top}}={T}_{e,\mathrm{top}}.\end{eqnarray} \tag{ 91 }$$

The electron and ion sound speeds at the top of the column are given by (see Equation (37))

$$\begin{eqnarray}&&{a}_{e,\mathrm{top}}^{2}=\displaystyle \frac{{\gamma }_{e}{{kT}}_{e,\mathrm{top}}}{{m}_{\mathrm{tot}}},\qquad {a}_{i,\mathrm{top}}^{2}=\displaystyle \frac{{\gamma }_{i}{{kT}}_{i,\mathrm{top}}}{{m}_{\mathrm{tot}}},\end{eqnarray} \tag{ 92 }$$

and therefore our assumption that ${T}_{i,\mathrm{top}}={T}_{e,\mathrm{top}}$ leads to the relation

$$\begin{eqnarray}&&{a}_{e,\mathrm{top}}={\left(\displaystyle \frac{{\gamma }_{e}}{{\gamma }_{i}}\right)}^{1/2}{a}_{i,\mathrm{top}}.\end{eqnarray} \tag{ 93 }$$

The radial component of the radiation energy flux, averaged over the cross-section of the column at radius r, is given by

$$\begin{eqnarray}&&{F}_{r}(r)=-\displaystyle \frac{c}{3{n}_{e}{\sigma }_{\parallel }}\,\displaystyle \frac{{{dU}}_{r}}{{dr}}+\displaystyle \frac{4}{3}{{vU}}_{r},\end{eqnarray} \tag{ 94 }$$

where the first term on the right-hand side represents the upward diffusion of radiation energy parallel to the magnetic field, and the second term represents the downward advection of radiation energy toward the stellar surface (with v < 0). The fact that the top of the accretion column is the last scattering surface implies photon transport makes a transition from diffusion to free streaming at r = r_top, so that we make the following replacement in Equation (94),

$$\begin{eqnarray}&&-\displaystyle \frac{c}{3{n}_{e}{\sigma }_{\parallel }}\,\displaystyle \frac{{{dU}}_{r}}{{dr}}\to c\,{U}_{r},\ \ \ \ r\to {r}_{\mathrm{top}}.\end{eqnarray} \tag{ 95 }$$

By incorporating this transition into Equation (94), we see that the radiation energy flux at the upper surface is given by

$$\begin{eqnarray}&&{F}_{r}({r}_{\mathrm{top}})=\left(c+\displaystyle \frac{4}{3}{v}_{\mathrm{top}}\right){U}_{r}({r}_{\mathrm{top}}).\end{eqnarray} \tag{ 96 }$$

The form of the total energy transport rate is derived from Equation (43), using Equations (35), (36), (39), and (40), which yields

$$\begin{eqnarray}\dot{E}(r) & = & A(r)F(r)\\ & = & \dot{M}\left(\displaystyle \frac{{F}_{r}}{\rho | v| }-\displaystyle \frac{{v}^{2}}{2}-\displaystyle \frac{{a}_{i}^{2}}{{\gamma }_{i}-1}-\displaystyle \frac{{a}_{e}^{2}}{{\gamma }_{e}-1}+\displaystyle \frac{{{GM}}_{* }}{r}\ \right).\end{eqnarray} \tag{ 97 }$$

The expression for the total energy transport rate at r = r_top is simplified once we implement the free-streaming boundary condition in Equation (96), and use Equations (35) and (36) to substitute for the radiation energy density U_r in terms of the radiation sound speed a_r. The result obtained is

$$\begin{eqnarray}{{\dot{E}}_{\mathrm{top}}\equiv \dot{E}| }_{r={r}_{\mathrm{top}}} & = & -\,\dot{M}\left[\left(\displaystyle \frac{{\gamma }_{i}}{{\gamma }_{i}-1}+\displaystyle \frac{{\gamma }_{e}}{{\gamma }_{e}-1}\right)\displaystyle \frac{{a}_{i,\mathrm{top}}^{2}}{{\gamma }_{i}}\right.\\ & & +\left.\left(\displaystyle \frac{c}{{v}_{\mathrm{top}}}+\displaystyle \frac{4}{3}\right)\displaystyle \frac{{a}_{r,\mathrm{top}}^{2}}{{\gamma }_{r}({\gamma }_{r}-1)}\right],\end{eqnarray} \tag{ 98 }$$

where we have also utilized Equations (88) and (93).

The ionized gas flows downward after entering the top of the accretion funnel at radius r = r_top, and eventually passes through a standing, radiation-dominated shock, where most of the kinetic energy is radiated away through the walls of the accretion column (Becker 1998). Below the shock, the gas passes through a sinking regime, where the remaining kinetic energy is radiated away (Basko & Sunyaev 1976). Ultimately, the flow stagnates at the stellar surface, and the accreting matter merges with the stellar crust.

The surface of the neutron star is too dense for radiation to penetrate significantly (Lenzen & Trümper 1978), and therefore the diffusion component of the radiation energy flux must vanish there. Furthermore, due to the stagnation of the flow at the stellar surface, the advection component should also vanish, and consequently we conclude that the radiation energy flux ${F}_{r}\to 0$ as $r\to {R}_{* }$. We refer to this as the "mirror" surface boundary condition, which can be written as

$$\begin{eqnarray}&&{{F}_{r}(r)| }_{r={R}_{* }}=0.\end{eqnarray} \tag{ 99 }$$

The stagnation of the flow at the stellar surface also implies there is no flux of kinetic energy into the star. Hence, at the stellar surface, the total energy transport rate, $\dot{E}$, reduces to the addition of (negative) gravitational potential energy to the star. The surface boundary condition for the total energy transport rate is therefore given by (see Equations (43) and (44))

$$\begin{eqnarray}&&{\dot{E}(r)| }_{r={R}_{* }}=\displaystyle \frac{{{GM}}_{* }\dot{M}}{{R}_{* }}.\end{eqnarray} \tag{ 100 }$$

The stagnation boundary condition formally requires that v = 0 at the stellar surface, where r = R_*. However, in practice, it is not possible to perfectly satisfy this condition due to the divergence of the mass density ρ implied by stagnation. Therefore, we approximate stagnation at the stellar surface in our simulations using the condition

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{r\to {R}_{* }}| v(r)| \lesssim \ 0.01\,c.\end{eqnarray} \tag{ 101 }$$

As the flow decelerates near the base of the accretion column, the density increases and the opacity becomes dominated by free–free absorption, leading to the formation of a dense "thermal mound" (e.g., Davidson 1973). The thermal mound, with a temperature between 10⁷ K and 10⁸ K, is the source of the blackbody seed photons that scatter throughout the column and contribute to the emergent Comptonized spectrum. The upper surface of the thermal mound is located at radius r = r_th, which is defined as the radius at which the Rosseland mean of the free–free optical depth, ${\tau }_{\parallel }^{\mathrm{ff}}$, measured from the top of the column, is equal to unity.

In general, the vertical variation of ${\tau }_{\parallel }^{\mathrm{ff}}$ is computed using the integral

$$\begin{eqnarray}&&{\tau }_{\parallel }^{\mathrm{ff}}(r)={\int }_{r}^{{r}_{\mathrm{top}}}{\alpha }_{{\rm{R}}}(r^{\prime} )\,{dr}^{\prime} ,\end{eqnarray} \tag{ 102 }$$

where α_R is the Rosseland mean free–free absorption coefficient for fully ionized hydrogen. Equation (102) implies that the Rosseland mean free–free optical depth at the top of the column is zero, so that

$$\begin{eqnarray}&&{\tau }_{\parallel }^{\mathrm{ff}}({r}_{\mathrm{top}})=0.\end{eqnarray} \tag{ 103 }$$

At the upper surface of the thermal mound, we have

$$\begin{eqnarray}&&{\tau }_{\parallel }^{\mathrm{ff}}({r}_{\mathrm{th}})={\int }_{{r}_{\mathrm{th}}}^{{r}_{\mathrm{top}}}{\alpha }_{{\rm{R}}}(r^{\prime} ){dr}^{\prime} =1.\end{eqnarray} \tag{ 104 }$$

Between the thermal mound surface and the stellar surface, ${\tau }_{\parallel }^{\mathrm{ff}}\gt 1$, leading to an approximate balance between thermal emission and absorption, though the balance is not perfect due to the escape of photons through the sides of the accretion column. The various thermal transfer rates and corresponding timescales are further discussed in Section 7.6.

The computational domain for the calculation extends from the stellar surface, at dimensionless radius $\tilde{r}=4.836$, up to the top of the accretion column, at radius $\tilde{r}={\tilde{r}}_{\mathrm{top}}$, where we have assumed a canonical stellar mass of M_* = 1.4 M_⊙ and a radius of R_* = 10 km. The attainment of a completely self-consistent description of the hydrodynamic structure of the accretion column, along with the radiation spectrum, is achieved using an iterative procedure. The coupling between the hydrodynamical simulation performed in Mathematica and the spectrum calculation performed in COMSOL is made via three vectors of information, which are passed between the two computational environments. In order to compute the dynamical structure in Mathematica, we require knowledge of the inverse-Compton temperature function, $g(\tilde{r})$ (see Equation (75)). Conversely, in order to carry out the spectrum calculation in COMSOL, we require knowledge of the velocity and electron temperature profiles, $\tilde{v}(\tilde{r})$ and ${T}_{e}(\tilde{r})$, respectively.

The iteration procedure begins with a calculation of the "0th" hydrodynamical structure in Mathematica, which is generated by arbitrarily setting $g(\tilde{r})=1$, meaning that we are initially assuming the inverse-Compton temperature ${T}_{\mathrm{IC}}(\tilde{r})$ is exactly equal to the electron temperature ${T}_{e}(\tilde{r})$ for all $\tilde{r}$ along the column. Once the six free model parameters listed in Table 1 are assigned provisional values, the system of five coupled ODEs is integrated in Mathematica to determine the first approximation of the dynamical structure of the column. The resulting accretion velocity profile, $\tilde{v}(\tilde{r})$, and electron temperature profile, ${T}_{e}(\tilde{r})$, are then exported from Mathematica and passed into the COMSOL multiphysics module in preparation for the computation of the phase-averaged radiation distribution inside the column.

The COMSOL multiphysics module is a computer environment that employs the finite element method (FEM) and is well-suited for solving the radiation transport equation, which is a second order, elliptical, nonlinear partial differential equation. COMSOL inputs the electron temperature and accretion velocity profiles from Mathematica and then solves the photon transport equation on a meshed grid using the boundary conditions discussed in Section 4. The resulting photon distribution function $f(\tilde{r},\epsilon )$ (photons cm⁻³ erg⁻³) and phase-averaged photon count rate spectrum ${F}_{\epsilon }(\epsilon )$ (photons s⁻¹ cm⁻² keV⁻¹) are obtained and discussed in Paper II, where is the photon energy. All transport phenomena are calculated using $f(\tilde{r},\epsilon )$, including the radiation flux F_r, the radiation energy density U_r, and the photon number density n_ph. By exploiting the combined strengths of Mathematica and COMSOL, we are able to solve, for the first time to our knowledge, the complete self-consistent problem of spectral formation and radiation hydrodynamics in an X-ray pulsar accretion column. We briefly discuss some aspects of the dual-platform iteration and the related convergence criteria below, but we defer complete details on the COMSOL calculation to Paper II.

Although we do not present detailed spectral results in this paper, it is important to highlight our method for treating cyclotron absorption here, since this process plays a significant role in determining the shape of the simulated spectrum, which is compared to the observational data in order to finalize the model parameters. In lieu of a detailed model for the formation of cyclotron absorption features in the envelopes of pulsar accretion columns, which has not yet been developed, we will treat the formation of the observed absorption features by supposing that the features are imprinted at a particular altitude, denoted by r_cyc. Hence the centroid energy of the absorption feature is interpreted as the cyclotron energy corresponding to the dipole magnetic field strength at radius r_cyc in the column. We argue that this approach is reasonable, provided the cyclotron imprint radius r_cyc is close to the radius at which the X-ray luminosity per unit length along the column, ${{\mathscr{L}}}_{r}$, is maximized, where ${{\mathscr{L}}}_{r}{dr}$ is the energy emitted per unit time through the walls of the dipole-shaped volume of the accretion column between positions r and r + dr.

We can derive an expression for ${{\mathscr{L}}}_{r}$ by noting that in our escape-probability formalism, the energy escaping through the walls of the accretion column between radii r and r + dr per unit time is given by

$$\begin{eqnarray}&&{{\mathscr{L}}}_{r}{dr}={\dot{U}}_{\mathrm{esc}}A(r){dr},\end{eqnarray} \tag{ 120 }$$

where ${\dot{U}}_{\mathrm{esc}}$ is given by Equation (55) and the cross-sectional area of the column is A(r) = Ω(r)r² (see Equation (39)). Solving for ${{\mathscr{L}}}_{r}$ yields

$$\begin{eqnarray}&&{{\mathscr{L}}}_{r}=\displaystyle \frac{{U}_{r}(r){\rm{\Omega }}(r){r}^{2}}{{t}_{\mathrm{esc}}(r)}.\end{eqnarray} \tag{ 121 }$$

We denote the radius of maximum X-ray emission using r_X. In our approach, we attempt to minimize the distance between r_X and the cyclotron imprint radius, r_cyc. Out of the three sources treated here, Her X-1 is the only one in which the cyclotron absorption radius r_cyc is exactly equal to r_X. In the other two sources, Cen X-3 and LMC X-4, the two radii deviate by about 10%.

Our method for determining the convergence of the solutions for the flow velocity $\tilde{v}(\tilde{r})$, the sound speeds ${\tilde{a}}_{i}(\tilde{r})$, ${\tilde{a}}_{e}(\tilde{r})$, and ${\tilde{a}}_{r}(\tilde{r})$, and the energy transport rate $\tilde{{\mathscr{E}}}(\tilde{r})$, is based on the comparison of successive iterates of the electron temperature, T_e, and the inverse-Compton temperature, T_IC. We define the convergence ratios, ${{\mathscr{R}}}_{e}$ and ${{\mathscr{R}}}_{\mathrm{IC}}$, respectively, for the electron and inverse-Compton temperatures using

$$\begin{eqnarray}&&{{\mathscr{R}}}_{e}^{n+1}\equiv \displaystyle \frac{{T}_{e}^{n+1}}{{T}_{e}^{n}},\qquad {{\mathscr{R}}}_{\mathrm{IC}}^{n+1}\equiv \displaystyle \frac{{T}_{\mathrm{IC}}^{n+1}}{{T}_{\mathrm{IC}}^{n}},\end{eqnarray} \tag{ 122 }$$

where the superscripts represent the iteration number for the corresponding solution vectors. The solutions are deemed to have converged when the vector of convergence ratios for both the electron and the inverse-Compton temperature profiles are within 1% of unity across the entire computational grid.

As explained in Section 5.1, we obtain the solution for the "0th" iteration for the dynamical structure by setting $g(\tilde{r})=1$ across the grid in the Mathematica calculation, and we then pass the resulting velocity profile $\tilde{v}(\tilde{r})$ and electron temperature profile ${T}_{e}(\tilde{r})$ into the COMSOL platform in order to obtain the corresponding "0th" iteration of the photon distribution function, $f(\tilde{r},\epsilon )$. Once the solution for $f(\tilde{r},\epsilon )$ has been obtained using COMSOL, the associated profile of the inverse-Compton temperature, ${T}_{\mathrm{IC}}(\tilde{r})$, is computed using Equation (73), which is then combined with the electron temperature profile ${T}_{e}(\tilde{r})$ to obtain the new iteration of the temperature ratio function, $g(\tilde{r})$, using Equation (75). Subsequently, the new iterate for $g(\tilde{r})$ is used as input into the Mathematica implementation to compute new results for the dynamical structure variables, and so on.

This iterative cycle is continued, and the convergence ratios between successive iterates are computed using Equation (122), until convergence is achieved, which operationally means that the convergence ratios for the electron and inverse-Compton temperature profiles differ from unity by less than 1% at all radii in the column. In the end, once convergence is achieved, we have obtained a self-consistent set of results for the radiation distribution $f(\tilde{r},\epsilon )$ and the five dynamical variables $\tilde{v}(\tilde{r})$, ${\tilde{a}}_{r}(\tilde{r})$, ${\tilde{a}}_{i}(\tilde{r})$, ${\tilde{a}}_{e}(\tilde{r})$, and $\tilde{{\mathscr{E}}}(\tilde{r})$. In the following section, we discuss the application of the method to compute the structure of the accretion column for three specific accretion-powered X-ray pulsars.

Figure 4 depicts the results obtained for the accretion column structure upon applying our model to Her X-1. The dynamical variables plotted include the bulk fluid velocity and the radiation sound speed, the gas and radiation pressures, the pressure gradients, the Mach numbers, the temperatures, the energy transport per unit mass, the bulk fluid density, and the parallel scattering and parallel absorption optical depths. We adopt for the source luminosity L_X = 2 × 10³⁷ erg s⁻¹ (Reynolds et al. 1997; Dal Fiume et al. 1998). The values of the six model free parameters in the Her X-1 simulation are $\bar{\sigma }/{\sigma }_{{\rm{T}}}=2.600\times {10}^{-3}$, σ_∥/σ_T = 1.024 × 10⁻³, ℓ₁ = 0 m, ℓ₂ = 125 m, ${{\mathscr{M}}}_{r0}=4.07$, and B_* = 6.25 × 10¹² G. Note that the top of the accretion column in the graphs of Figure 4 is located on the right side, and the stellar surface is located on the left side.

Zoom In Zoom Out Reset image size

The accretion column for Her X-1 is completely filled with inflowing plasma, which makes this the only completely filled column among the three sources investigated here. This may be reasonable, since Her X-1 is a "fast rotator," as discussed in Section 7.4. The upper limit for the outer radius is ℓ₂ ≲ 223 m, given by Equation (12), which is almost double the 125 m outer polar cap radius used in our model. In the case of Her X-1, the accretion column spans a length of 11.20 km, and the bulk free-fall velocity at the top of the column (Equation (84)) is equal to 0.442 c.

The radiation sound speed at the top of the column is derived using Equation (114), and the radiation sonic surface (where ${{\mathscr{M}}}_{r}=1$) is located at a radius of r_sonic = 11.95 km, which is where the bulk fluid slows to less than the radiation sound speed. The onset of stagnation is most noticeable when the bulk fluid enters the extended sinking regime (Basko & Sunyaev 1976), which begins approximately 700 m above the surface, and is characterized by a gradually decelerating flow, accompanied by a corresponding increase in temperature, pressure, and density. Approximate stagnation occurs at the stellar surface, with a residual bulk velocity of 0.0084 c.

It is apparent from the Mach number profiles plotted in Figure 4 that the flow remains supersonic with respect to the gas at the lower boundary of our computational domain, which is located just above the stellar surface. Hence we would expect a final discontinuous shock transition to occur as the flow merges into the stellar crust. However, the amount of residual kinetic energy converted into radiation at the discontinuous shock is negligible compared to the energy loss associated with the radiation emitted farther up in the column. Similar behavior is also observed in the cases of LMC X-4 and Cen X-3. Hence our neglect of the discontinuous shock is reasonable for the luminous sources treated here. However, the effect of the discontinuous shock is likely to be more important in lower-luminosity sources such as X Persei (e.g., Langer & Rappaport 1982).

The model for Her X-1 required 16 iterations before T_e and T_IC stabilize to less than a 1% change from the previous iteration. Electron and ion temperatures are in near thermal equilibrium throughout the column, and at the top of the column, we have T_e = T_i = 1.41 × 10⁷ K. The inverse-Compton temperature at the top of the column, T_IC = 6.82 × 10⁷ K, is almost five times larger than the electron temperature, which is the largest temperature gap between the photons and gas at any radius in the column. The electron temperature at the stellar surface is found to be 6.91 × 10⁷ K. Further discussion of the temperature distributions and the related thermal and dynamical timescales is presented in Section 7.

Energy transport per unit mass transport is plotted in Figure 4 in terms of the dimensionless quantity $\dot{E}/(\dot{M}{c}^{2})$. According to our sign convention, a negative value corresponds to energy flow downward, toward the stellar surface (see Equation (84)), and therefore the profile of the gravitational potential energy component, ${\dot{E}}_{g}$, is depicted as a positive value, given by

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{E}}_{g}}{\dot{M}{c}^{2}}=\displaystyle \frac{{{GM}}_{* }\dot{M}}{r}\displaystyle \frac{1}{\dot{M}{c}^{2}}=\displaystyle \frac{1}{\tilde{r}}.\end{eqnarray} \tag{ 123 }$$

At the stellar surface, the value of the dimensionless radius is $\tilde{r}={\tilde{r}}_{* }=4.836$, assuming canonical values for the stellar mass and radius, with ${M}_{* }=1.4\,{M}_{\odot }$ and R_* = 10 km. The kinetic energy transport component dominates over the radiation component at the top of the column, while the two are nearly equal at the radiation sonic surface. However, in the sinking regime, the kinetic energy is negligible, and we see that the energy transport is dominated by radiation advection and diffusion.

According to Equation (100), at the surface of the star, we expect the total energy transport rate to reduce to the gravitational component only, so that $\dot{E}/(\dot{M}{c}^{2})=1/{\tilde{r}}_{* }=0.2068$. Hence the radiation energy flux should vanish at the stellar surface (the "mirror" condition), and that is the boundary condition we attempt to enforce. Her X-1 is the only source in which the mirror condition at the stellar surface is slightly relaxed. The Her X-1 results depicted in Figure 4 indicate that the advective and diffusive components in Equation (94) do not exactly cancel at the surface, and we obtain for the residual total energy transport rate $\dot{E}/(\dot{M}{c}^{2})=0.228$. This represents an error of ∼10% from the purely gravitational component. Her X-1 is also the only source in which the radius at which the cyclotron absorption feature is imprinted, r_cyc = 11.74 km, is exactly equal to the radius of maximum emission, r_X. See Section 5.2 for further details.

The scattering and absorption optical depths, given by Equations (82) and (104), respectively, are plotted in Figure 4 for photons propagating parallel to the magnetic field. The top of the column is the last scattering surface before photons freely escape in the vertical direction, and therefore the scattering optical depth is low in that region, and gradually increases toward the bottom of the column. The scattering optical depth diverges exponentially in the sinking regime, reflecting the pileup of material at the stellar surface where the plasma density becomes extremely large. On the other hand, the parallel absorption optical depth never reaches unity in the case of Her X-1. Therefore, the thermal mound is located just above the stellar surface and blackbody photons are produced at the electron surface temperature, which may help explain the positive radiation diffusion flux at the stellar surface for this source.

The results we obtain for the dimensions of the hot spot at the magnetic pole in Her X-1 are significantly different than those obtained by BW07. In particular, we find that the outer polar cap radius is 125 m, whereas BW07 found that their cylindrical polar cap radius was r₀ = 44 m, which is about three times smaller than our result. Furthermore, BW07 assumed a constant electron temperature of 6.25 × 10⁷ K, which is about 9.5% lower than our stellar surface temperature of 6.91 × 10⁷ K. Another significant difference is that our stellar surface B-field strength of B_* = 6.25 × 10¹² G is nearly double the BW07 value of B_* = 3.80 × 10¹² G. We believe the differences between our results and those of BW07 probably reflect the fact that BW07 assumed constant values for T_e, r₀, and B, which means that they should be interpreted as average values in an actual accretion column. On the other hand, our model implements a realistic dipole magnetic field geometry, with varying electron and ion temperatures, and therefore our surface results would be expected to exceed the mean values taken from the BW07 model.

The profiles of the dynamical variables obtained in our application to Cen X-3 are depicted in Figure 5. The six free parameter values used in the Cen X-3 simulation are $\bar{\sigma }/{\sigma }_{{\rm{T}}}=3.000\times {10}^{-3}$, σ_∥/σ_T = 7.510 × 10⁻³, ℓ₁ = 657 m, ℓ₂ = 750 m, ${{\mathscr{M}}}_{r0}=6.151$, and B_* = 3.6 × 10¹² G. The source luminosity L_X = 2.8 × 10³⁸ erg s⁻¹ is the same value used by BW07, which provides an opportunity to directly compare our model results with theirs using the same accretion rate. We find that the accretion column in Cen X-3 is a hollow cavity, with a wall thickness of 93 m at the stellar surface. The upper limit for the outer radius is ℓ₂ ≲ 761 m, according to Equation (12). The accretion column spans a length of 14.40 km, and the bulk velocity at the top of the column is equal to the free-fall velocity of 0.412 c. The radiation sonic surface is located at a radius of r_sonic = 12.21 km, and the bulk fluid enters the sinking regime at an altitude of 1.1 km above the surface. A thermal mound exists for Cen X-3 at an altitude of 590 m above the stellar surface, where the parallel optical depth exceeds unity, and the electron temperature is 5.73 × 10⁷ K. Approximate bulk stagnation occurs at the stellar surface, with a residual velocity of 0.0081 c and a surface electron temperature equal to 6.97 × 10⁷ K, in contrast to the constant electron temperature T_e = 3.40 × 10⁷ K used in the corresponding BW07 model. Eight iterations were required before T_e and T_IC converged in our model. The stellar surface mirror condition is satisfied, so that the radiation energy flux essentially vanishes, and the total energy flux reduces to the gravitational component only. In the case of Cen X-3, the radius at which the cyclotron absorption feature is imprinted on the spectrum is r_cyc = 10.94 km, whereas the radius of maximum emission is r_X = 12.02 km.

Zoom In Zoom Out Reset image size

Figure 6 depicts the results we obtain for the dynamical profiles upon applying our model to LMC X-4, using for the source luminosity L_X = 3.9 × 10³⁸ erg s⁻¹ (La Barbera et al. 2001). The six model free parameter values are $\bar{\sigma }/{\sigma }_{{\rm{T}}}=2.500\times {10}^{-3}$, σ_∥/σ_T = 4.176 × 10⁻³, ℓ₁ = 547 m, ℓ₂ = 650 m, ${{\mathscr{M}}}_{r0}=2.761$, and B_* = 8.00 × 10¹² G. The accretion column for LMC X-4 is hollow, exhibiting a geometry similar to that found in Cen X-3. The upper limit for the outer radius is ℓ₂ ≲ 633 m, according to Equation (12). The accretion column spans a length of 11.30 km, and the bulk velocity at the top of the column has a local free-fall velocity equal to 0.44 c. The radiation sonic point is located at radius r_sonic = 13.21 km, and the bulk fluid enters the sinking regime at an altitude of 1.4 km above the stellar surface. The thermal mound is located 530 m above the surface, where the electron temperature is equal to 7.59 × 10⁷ K. Approximate stagnation occurs at the surface with a residual velocity of 0.0098 c and a surface temperature equal to 8.75 × 10⁷ K, and the stellar surface mirror condition is satisfied. Our value for the electron surface temperature is somewhat higher than the BW07 model, which used the constant value T_e = 5.90 × 10⁷ K. The model required eight iterations to converge T_e and T_IC, and the radius of maximum emission occurs at r_X = 12.94 km, whereas the cyclotron absorption feature is imprinted at radius r_cyc = 14.26 km.

Zoom In Zoom Out Reset image size

Our model uses the canonically accepted values of R_* = 10⁶ cm for stellar radius, M_* = 1.4 M_⊙ for stellar mass, and ${\sigma }_{\perp }={\sigma }_{{\rm{T}}}=6.652\times {10}^{-25}\,{\mathrm{cm}}^{2}$ for the perpendicular electron scattering cross-section. Six free parameters uniquely determine the dynamical properties of the column, namely (1) the angle-averaged scattering cross-section $\bar{\sigma }$, (2) the inner polar cap arc-radius ℓ₁, (3) the outer polar cap arc-radius ℓ₂, (4) the radius at the top of the accretion column ${\tilde{r}}_{\mathrm{top}}$, (5) the incident radiation Mach number ${{\mathscr{M}}}_{r0}$, and (6) the stellar surface magnetic field strength, B_*, which is established at the magnetic pole. We carry out a self-consistent and iterative calculation of a set of five coupled conservation equations in order to self-consistently compute the hydrodynamic structure of the accretion column and the emergent radiation spectrum.

The dynamical structure is determined by using Mathematica to solve the dynamical equations, combined with appropriate physical boundary conditions. In this paper, we have presented results describing the detailed structure of the accretion columns in the well-known and high-luminosity X-ray pulsars Her X-1, Cen X-3, and LMC X-4. The solutions shown in Figures 4–6 provide the radial profiles of the bulk flow velocity $\tilde{v}$, the radiation sound speed ${\tilde{a}}_{r}$, the ion sound speed ${\tilde{a}}_{i}$, the electron sound speed ${\tilde{a}}_{e}$, and the total energy transport rate per unit mass $\tilde{{\mathscr{E}}}$.

We find that the dynamical effects of gas pressure are negligible and that radiation pressure decelerates the gas to rest at the stellar surface in all three sources considered here, which all have relatively high accretion rates ($\dot{M}\gtrsim {10}^{16}\,{\rm{g}}\,{{\rm{s}}}^{-1}$). However, the inclusion of the gas energy equation is essential in our calculation of the ion and electron temperature profiles, which show significant deviations from the profile of the inverse-Compton temperature, especially near the top of the accretion column, where conditions are somewhat farther from equilibrium than in the deeper portions of the column.

A noticeable effect of increasing the surface magnetic pole field strength, B_*, is an increase of the electron temperature at the stellar surface. We shall see in Paper II that this tends to harden the phase-averaged spectrum, and raises the energy of the exponential cutoff due to thermal Comptonization. The value of B_* was determined by minimizing the distance between the location of the maximum emission from the walls of the column and the location at which the cyclotron absorption feature is imprinted on the spectrum. This is, admittedly, a rather crude criterion, but it is adequate for our purposes here, pending the availability of a generalized model that includes a rigorous treatment of the formation of the cyclotron absorption feature in the outer sheath of the column (Schönherr et al. 2008, 2014). Furthermore, our magnetic field follows the correct dipole variation with radius, whereas BW07 assumed a constant magnetic field. Since BW07 essentially utilized an average value of the magnetic field within the column, it is reasonable that our surface field value would exceed their (constant) value.

It should also be emphasized that the distance between the radius of maximum emission, r_X, and the radius at which the cyclotron absorption feature is imprinted, r_cyc, only agree closely in the case of Her X-1. In the other two sources, the disagreement between these two radii can be as large as ∼10%. A complete treatment of this issue will require the development of a more generalized code that simultaneously treats the hydrodynamics and the effect of cyclotron absorption occurring along the entire length of the column, which is beyond the scope of the present paper.

We find that the vertical extent of the accretion column is between 10 and 12 km for all three sources. This is obviously comparable to the stellar radius, and therefore it is essential to implement the dipole variation of the magnetic field, both in order to compute the correct magnetic field variation, and also to properly compute the cross-sectional area of the column. Equation (87) establishes the boundary condition at the upper surface of the column relating r_top and σ_∥, which is given by

$$\begin{eqnarray}&&{r}_{\mathrm{top}}={\left[{r}_{c}^{-3/2}+\displaystyle \frac{3{m}_{\mathrm{tot}}{(2{{GM}}_{* })}^{1/2}}{2{\sigma }_{\parallel }\dot{M}}\displaystyle \frac{{{\rm{\Omega }}}_{* }}{{R}_{* }}\right]}^{-2/3}.\end{eqnarray} \tag{ 124 }$$

The top of the cylindrical column in the BW07 model is given by combining their Equations (26) and (80) to yield

$$\begin{eqnarray}&&{r}_{\max }^{\mathrm{BW}}=\displaystyle \frac{{R}_{* }}{2}\left[{\left(1+\displaystyle \frac{4{m}_{p}{{GM}}_{* }\pi {r}_{0}^{2}}{\alpha c{\sigma }_{\parallel }\dot{M}{R}_{* }^{2}}\right)}^{1/2}-1\right],\end{eqnarray} \tag{ 125 }$$

where r₀ is the radius of the cylinder and α ∼ 0.3–0.4 is a constant. We compare Equation (124) to Equation (125) and note the similarities with the purely cylindrical model from BW07 with respect to dependence on r₀ = ℓ₂, $\dot{M}$, and σ_∥. Our model differs because we allow for a hollow-column geometry that alters the value of Ω_* according to Equation (7). Additionally, Equation (125) relies on the value of α, which was required to be introduced in the BW07 velocity profile in order to solve the photon transport equation using the separation of variables method (Lyubarskii & Syunyaev 1982).

The column-top radii for the three BW07 models using Equation (125) are ${r}_{\max }^{\mathrm{BW}}=13.72$ km (Her X-1), 29.97 km (Cen X-3), and 32.44 km (LMC X-4), respectively. The larger differences in accretion column length for Cen X-3 and LMC X-4 can be attributed to B-field geometry. Whereas the BW07 model is based strictly on a cylindrical geometry, our new model is based on a realistic dipole geometry that matches not only the free-fall velocity, but also the derivative of the free-fall velocity, thereby satisfying the momentum conservation equation. The dynamical profiles for all three of our sources show that the radiation sonic surface and the length of the sinking regime increase proportionally to the source luminosity (Basko & Sunyaev 1976).

The free-streaming boundary condition at the top of the column can be verified by comparing the forces of gravity and radiation acting on the inflowing gas. The assumption at the top of the column is that the last photon–electron scattering events occur prior to photons freely escaping in the vertical direction. Therefore, we expect the upward radiation force on the incoming electrons to be much smaller than the gravitational pull on the bulk fluid. The Newtonian gravitational force per electron–ion couple at the top of the column is

$$\begin{eqnarray}&&{F}_{g}=-\displaystyle \frac{{{GM}}_{* }{m}_{\mathrm{tot}}}{{r}_{\mathrm{top}}^{2}},\end{eqnarray} \tag{ 126 }$$

where ${m}_{\mathrm{tot}}={m}_{i}+{m}_{e}$ is the combined mass of the electron and ion. The radiation force on incoming bulk fluid is equal to the product of the electron parallel scattering cross-section and the radiation pressure,

$$\begin{eqnarray}&&{F}_{r}={\sigma }_{\parallel }{P}_{r}.\end{eqnarray} \tag{ 127 }$$

We calculate the two forces from Equations (126) and (127) using the parallel scattering cross-section from Table 4, the column-top radii r_top from Table 6, and the radiation pressure for each source, respectively, at the top of the columns shown in Figures 4–6. The force relationships calculated are shown in Table 8. As expected, the upward radiation force is dominated by the downward gravitational force. This result strengthens the free-streaming argument in which the top of the column is the last scattering surface.

A second method to verify the free-streaming surface is to compare the bulk fluid ram pressure and the radiation pressure. The ram pressure at the top of the column is

$$\begin{eqnarray}&&{P}_{\mathrm{ram}}=\rho {v}_{\mathrm{ff}}^{2}.\end{eqnarray} \tag{ 128 }$$

The pressure relationships calculated are shown in Table 9. The ram pressure dominates at the top of the column for all three sources. Here, too, we conclude that the free-streaming condition is a valid assumption for all three sources.

As discussed in Section 2.3, the field lines connected to the inner and outer walls of the accretion column cross the accretion disk at the radii ${R}_{2,\mathrm{disk}}$ and ${R}_{1,\mathrm{disk}}$, respectively. Matter is picked up from the disk and entrained onto the magnetosphere at the Alfvén radius, R_A. As the star rotates, the magnetic latitude in the disk plane, α, oscillates between ±φ, where φ is the inclination angle between the magnetic and rotation axes of the star. As a result of the variation in α, the Alfvén radius in the disk, R_A, oscillates due to the oscillation of the magnetic field strength in the disk plane, B_disk (see Equations (22) and (20)). In addition, the geometry of the inclined pulsar magnetosphere causes the inner and outer disk-crossing radii, ${R}_{2,\mathrm{disk}}$ and ${R}_{1,\mathrm{disk}}$, respectively, to oscillate as well. The combination of these oscillations causes matter to be fed into different parts of the accretion column as the star rotates. Hence, matter is fed to the inner wall of the accretion column when ${R}_{1,\mathrm{disk}}={R}_{{\rm{A}}}$, and to the outer wall of the column when ${R}_{2,\mathrm{disk}}={R}_{{\rm{A}}}$.

By using Equations (22), (20), and (24) to evaluate R_A, ${R}_{1,\mathrm{disk}}$, and ${R}_{2,\mathrm{disk}}$ as functions of the disk-plane latitude α, we can attempt to determine the inclination angle of the system, φ, such that ${R}_{2,\mathrm{disk}}\lesssim {R}_{{\rm{A}}}\lesssim {R}_{1,\mathrm{disk}}$ during one spin of the star. We carry out this procedure for Cen X-3 and LMC X-4 in Figure 7. The process also requires selecting a value for the normalization constant ξ appearing in Equation (22), such that the minimum value of the Alfvén radius equals the maximum value of ${R}_{2,\mathrm{disk}}$, corresponding to α = 0. We find that ξ = 1.21 and ξ = 1.12 for Cen X-3 and LMC X-4, respectively, as depicted by the red-dashed curves in Figure 7. For comparison, we also plot the results obtained when ξ = 1.00 for both sources, which are indicated by the orange-dashed curves. Once the value of ξ is determined, the inclination angle φ is obtained by requiring that ${R}_{{\rm{A}}}={R}_{1,\mathrm{disk}}$ when α = ±φ. We find that φ = 22.6° and φ = 26.0° for Cen X-3 and LMC X-4, respectively. Our estimate of φ = 22.6° for Cen X-3 is close to the estimate of 18° provided by Kraus et al. (1996). Gas is transferred from the disk to the pulsar magnetosphere in the range of magnetic latitude labeled as the "capture region" in Figure 7. We denote the mean values of the oscillating radii R_A, ${R}_{1,\mathrm{disk}}$, and ${R}_{2,\mathrm{disk}}$ using $\langle {R}_{{\rm{A}}}\rangle$, $\langle {R}_{1,\mathrm{disk}}\rangle$, and $\langle {R}_{2,\mathrm{disk}}\rangle$, respectively, and we present values for these quantities in Table 12. The Her X-1 values in Table 12 were computed using ξ = 1.00.

Zoom In Zoom Out Reset image size

The size of the hot spot on the stellar surface must be understood in terms of the dynamics between accreted plasma gas and the neutron star magnetosphere. There is vast literature on the topic, which includes, but is not limited to, Lamb et al. (1973), Arons & Lea (1976, 1980), Elsner & Lamb (1976, 1977, 1984), Michel (1977a, 1977b, 1977c), Ghosh et al. (1977), Petterson (1977a, 1977b, 1977c), Ghosh & Lamb (1978, 1979a, 1979b), Lai (1999), Romanova et al. (2003), Pfeiffer & Lai (2004), Ikhsanov et al. (2012), and Kulkarni & Romanova (2013).

The general picture depends on the accretion scenario, which is often proposed in models pertaining to spherical (radial) accretion, or to inflow via a Keplerian disk. Mass transfer occurs across, or through, the magnetosphere, and is eventually aligned with the polar cap region via entrainment in the magnetic field lines, or, in the "fast" rotators, such as Her X-1, the plasma may be directly deposited far above the top of the accretion column from the plasma in a dense atmosphere. This picture is consistent with our results for Her X-1, which show that the accretion column is completely filled, whereas the columns for Cen X-3 and LMC X-4 are partially hollow. There are a number of additional factors to consider in these two differing topologies, which we discuss in further detail below.

Our model assumes the inflowing electrons and ions at the top of the accretion column are equilibrated, so that T_i0 ≈ T_e0 (Arons & Lea 1976). This situation is treated by Equation (5) from Arons & Lea (1976), which gives the infalling gas Mach number at the top of the column, for r ≪ R_A, as

$$\begin{eqnarray}&&{{\mathscr{M}}}_{\mathrm{infall}}=1.36\times {10}^{9}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{1/2}{(r{T}_{e})}^{-1/2},\end{eqnarray} \tag{ 129 }$$

where M_* is the stellar mass, M_⊙ is one solar mass, r is the radial distance (cm) from the X-ray source, and T_e is the electron temperature (K) at the top of the column. Table 10 compares the infalling Mach numbers (Equation (129)) with the actual value of the incident ion Mach number using Equation (155), where we set T_i0 = T_e0. The error is less than 10% for all three sources. We conclude that the large incident flow velocities in our models are expected, and momentum conservation and incident free-fall velocity are still rigorously implemented in our model.

Table 10.  Infalling Gas Mach Number

Source	rtop	Ti0	${{\mathscr{M}}}_{\mathrm{infall}}$	${{\mathscr{M}}}_{i0}$	Error
	km	107 K	Equation (129)	Equation (155)	%
Her X-1	21.19	1.95	250	274	8.8
Cen X-3	24.40	5.71	430	435	1.1
LMC X-4	21.30	1.95	249	238	4.6

Download table as:  ASCIITypeset image

We shall review the fundamental orbital parameters for each source. The corotation radius, R_co, is defined as the radius at which the Keplerian orbital period is equal to the star's spin period, P_rot, so that (Lamb et al. 1973; Elsner & Lamb 1977)

$$\begin{eqnarray}&&{R}_{\mathrm{co}}\equiv {\left(\displaystyle \frac{{{GM}}_{* }}{{{\rm{\Omega }}}_{\mathrm{rot}}^{2}}\right)}^{1/3}=1.5\times {10}^{8}{P}_{\mathrm{rot}}^{2/3}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{1/3}\,\mathrm{cm},\end{eqnarray} \tag{ 130 }$$

where ${{\rm{\Omega }}}_{\mathrm{rot}}=2\pi /{P}_{\mathrm{rot}}$. The Keplerian angular velocity at the (mean) Alfvén radius is given by

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{K}}}(\langle {R}_{{\rm{A}}}\rangle )={\left(\displaystyle \frac{{{GM}}_{* }}{{\langle {R}_{{\rm{A}}}\rangle }^{3}}\right)}^{1/2}.\end{eqnarray} \tag{ 131 }$$

The influence of the stellar rotation on the accreting gas in the magnetosphere depends on the relationship between the angular velocities Ω_rot and ${{\rm{\Omega }}}_{{\rm{K}}}(\langle {R}_{{\rm{A}}}\rangle )$. Slow rotators have ${{\rm{\Omega }}}_{\mathrm{rot}}\lesssim {{\rm{\Omega }}}_{{\rm{K}}}(\langle {R}_{{\rm{A}}}\rangle )$, and therefore ${R}_{\mathrm{co}}\gtrsim \langle {R}_{{\rm{A}}}\rangle$, which indicates that the force of gravity at $\langle {R}_{{\rm{A}}}\rangle$ is larger than the centripetal force. The opposite condition is satisfied by fast rotators, in which ${{\rm{\Omega }}}_{\mathrm{rot}}\gtrsim {{\rm{\Omega }}}_{{\rm{K}}}(\langle {R}_{{\rm{A}}}\rangle )$ and ${R}_{\mathrm{co}}\lesssim \langle {R}_{{\rm{A}}}\rangle$. The rotation properties of each source are shown in Table 11. The spin period, P_rot, is 1.24 s for Her X-1 (Vasco et al. 2013), 4.82 s for Cen X-3 (Raichur & Paul 2008), and 13.5 s for LMC X-4 (Levine et al. 2000). We see that Her X-1 is the only "fast" rotator, and therefore we expect that rotation will affect its accretion flow differently than in the slow rotators, such as Cen X-3 and LMC X-4 (Elsner & Lamb 1976).

Table 11.  Stellar Rotation Properties

Source	Prot	$\langle {R}_{{\rm{A}}}\rangle$	Rco	${{\rm{\Omega }}}_{{\rm{K}}}(\langle {R}_{{\rm{A}}}\rangle )$	Ωrot	Rotation
	s	km	km	rad s−1	rad s−1
Her X-1	1.24	5334	1940	1.11	5.07	fast
Cen X-3	4.82	1879	4780	5.29	1.30	slow
LMC X-4	13.50	2535	9510	3.38	0.46	slow

Download table as:  ASCIITypeset image

The simulated rotation of Cen X-3, with a magnetic inclination of 226, over half of a spin period, is depicted in Figure 8. We define the phase β = 0° as the direction to the companion star, so that the magnetic field axis and the axis of rotation are exactly aligned when viewed from the companion star's point of view. Therefore, when β = 0°, the magnetic pole is pointed toward the companion star as nearly as possible. Figure 8 displays eight instantaneous "snapshots," with the rotation phase angle increasing by 22.5° each time. The plots illustrate how the Alfvén radius, as well as the inner and outer disk-crossing radii, oscillate as the star rotates.

Zoom In Zoom Out Reset image size

Her X-1 has a completely filled column, and therefore the magnetic field associated with the inner polar angle (θ₁ = 0) extends to the accretion capture radius, r_acc, given by (Boroson et al. 2001)

$$\begin{eqnarray}&&1\times {10}^{6}\,\mathrm{km}\lesssim {r}_{\mathrm{acc}}\equiv \displaystyle \frac{2{{GM}}_{* }}{{v}_{\mathrm{rel}}^{2}}\lesssim 4\times {10}^{6}\,\mathrm{km},\end{eqnarray} \tag{ 132 }$$

where ${v}_{\mathrm{rel}}={({v}_{\mathrm{wind}}^{2}+{v}_{\mathrm{ns}}^{2})}^{1/2}$, the orbital velocity v_ns ∼ 169 km s⁻¹, and the wind velocity v_wind ∼ 300–600 km s⁻¹. The results in Table 12 suggest a high magnetic field inclination angle for Her X-1, in which the plane of the accretion disk is more likely aligned over the polar caps, rather than close to the stellar equatorial plane.

Table 12.  Magnetospheric Geometry and Disk Inclination Angle

Source	$\langle {R}_{1,\mathrm{disk}}\rangle$	$\langle {R}_{2,\mathrm{disk}}\rangle$	$\langle {R}_{{\rm{A}}}\rangle$	θ1	θ2	φ
	km	km	km	degrees	degrees	degrees
Her X-1	N/A	32000	5331	0	0.72	N/A
Cen X-3	2149	1650	1879	3.76	4.30	22.6
LMC X-4	3023	2142	2535	3.13	3.72	26.0

Download table as:  ASCIITypeset image

In the hollow-column slow rotators (Cen X-3 and LMC X-4), the supply of accreting plasma is confined to the magnetic field equatorial region (i.e., the angular momentum vectors of the B-field and accretion disk are nearly aligned), and the plasma may be forced to squeeze between the magnetic field lines in the equatorial plane via the Rayleigh–Taylor instability. In this scenario, the gas is eventually entrained onto the magnetic field as it flows toward the polar cap (Elsner & Lamb 1976). On the other hand, in the fast rotators, such as Her X-1, the accretion column is completely filled, which may be the result of descent of the polar cusps (Arons & Lea 1976; Michel 1977c; Elsner & Lamb 1984), in which plasma enters the cusp region due to density buildup and the related pressure gradient. Michel (1977c) showed that the external gas pressure at the top of the polar cusp is expected to be at least five times greater than that at the equator, and therefore the magnetosphere shape for Her X-1, if it does have a larger plasma atmosphere above the polar caps, may be very different than those for Cen X-3 and LMC X-4.

We also note the luminosity of Her X-1 is sufficiently high (>10³⁶ erg s⁻¹) that the reconnection of entrained magnetic fields may distort the large-scale structure of the magnetosphere to produce narrow open clefts (Arons & Lea 1976), thereby allowing plasma to flood the full column width above the polar cap. Furthermore, Pfeiffer & Lai (2004) showed that accretion disks can be highly warped with saturated tilt angles at high latitudes where the inner radius of the disk may naturally rotate directly above the polar region. Finally, Ikhsanov et al. (2012) conducted a study of spherical accretion scenarios in which the magnetic field of a companion star distorts the B-field flowing with the accreting material, thereby controlling the process via turbulent diffusion, which could potentially allow for accretion flow above the polar cap in fast rotators such as Her X-1.

A thorough analysis of the timescale and energy rate dynamics would not be possible without the fundamental inclusion of gas pressure in the conservation equations. Recall from Section 2.5 that the electrons are essentially confined along the magnetic field with a 1D Maxwell–Boltzmann energy distribution and an energy density per degree of freedom given by (1/2)n_ekT_e. The column dynamics are best understood by investigating the energy transfer rates and associated timescales in this formalism. The itemized list below provides definitions for eight of the key physical processes governing the thermodynamics of the energy transfer between the ions, electrons, and radiation. The timescales include the following:

• 1.  
  Escape timescale:

  $$\begin{eqnarray}&&{t}_{\mathrm{esc}}={{\ell }}_{\mathrm{esc}}/{w}_{\perp }.\end{eqnarray} \tag{ 133 }$$

  See Equation (55) for more discussion on the escape timescale. A photon will travel escape distance ℓ_esc, with diffusion velocity w_⊥, from the centerline of the accretion channel to the exterior wall.
• 2.  
  Thermal bremsstrahlng absorption timescale:

  $$\begin{eqnarray}&&{t}_{\mathrm{abs}}=1/({\alpha }_{{\rm{R}}}\,c).\end{eqnarray} \tag{ 134 }$$

  See Equation (63) for the Rosseland mean absorption coefficient α_R. The absorption timescale defines the average time before a photon is absorbed via bremsstrahlung thermal free–free absorption.
• 3.  
  Bulk Comptonization timescale:

  $$\begin{eqnarray}&&{t}_{\mathrm{bulk}\mathrm{Comp}}=1/({\boldsymbol{\nabla }}\cdot {\boldsymbol{v}}).\end{eqnarray} \tag{ 135 }$$

  Timescale over which the bulk fluid with velocity ${\boldsymbol{v}}$ is compressed and adds heat to the accreting gas. The divergence operator captures the compression dynamics.
• 4.  
  Electron scattering timescale:

  $$\begin{eqnarray}&&{t}_{\mathrm{scat}}=1/({n}_{e}{\sigma }_{{\rm{T}}}c).\end{eqnarray} \tag{ 136 }$$

  Photons scatter isotropically off electrons with Thomson cross-section ${\sigma }_{{\rm{T}}}$ inside the accretion column.
• 5.  
  Comptonization timescale:

  $$\begin{eqnarray}&&{t}_{\mathrm{Comp}}=\displaystyle \frac{(1/2)\,{n}_{e}{{kT}}_{e}}{{\dot{U}}_{\mathrm{Comp}}}.\end{eqnarray} \tag{ 137 }$$

  Energy is added or removed from electrons due to interactions with photons within a Comptonization time interval t_Comp. The rate of energy transfer, ${\dot{U}}_{\mathrm{Comp}}$, is given by Equation (74).
• 6.  
  Bremsstrahlung emission timescale:

  $$\begin{eqnarray}&&{t}_{\mathrm{brem}}=\displaystyle \frac{(1/2)\,{n}_{e}{{kT}}_{e}}{{\dot{U}}_{\mathrm{ff}}}.\end{eqnarray} \tag{ 138 }$$

  Energy is removed from the electrons due to the bremsstrahlung cooling rate ${\dot{U}}_{\mathrm{ff}}$ given by Equation (61).
• 7.  
  Cyclotron emission timescale:

  $$\begin{eqnarray}&&{t}_{\mathrm{cyc}}=\displaystyle \frac{(1/2)\,{n}_{e}{{kT}}_{e}}{{\dot{U}}_{\mathrm{cyc}}}.\end{eqnarray} \tag{ 139 }$$

  Energy is removed from the electrons due to the cyclotron cooling rate ${\dot{U}}_{\mathrm{cyc}}$ given by Equation (67).
• 8.  
  The electron–ion equilibration timescale (based on Equations (32) and (103) from Arons et al. 1987):

  $$\begin{eqnarray}{t}_{\mathrm{ei}} & = & 3.2\times {10}^{-10}\\ & & \times {\left(\displaystyle \frac{{n}_{e}}{{10}^{22}{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\displaystyle \frac{{{kT}}_{\mathrm{eff}}}{10\mathrm{keV}}\right)}^{3/2}\left(\displaystyle \frac{9}{\mathrm{ln}\,{\rm{\Lambda }}}\right)\,{\rm{s}},\end{eqnarray} \tag{ 140 }$$

  where T_eff is given by Equation (78), and

  $$\begin{eqnarray}&&{\rm{\Lambda }}=9.1-\displaystyle \frac{1}{2}\mathrm{ln}\left(\displaystyle \frac{{n}_{e}}{{10}^{22}\,{\mathrm{cm}}^{-3}}\right)+\mathrm{ln}\left(\displaystyle \frac{{{kT}}_{\mathrm{eff}}}{10\,\mathrm{keV}}\right).\end{eqnarray} \tag{ 141 }$$

An analysis of Equation (140) shows that the electron–ion equilibration timescale is at least an order of magnitude smaller than all other timescales, which explains why the electron and ion are nearly in thermal equilibrium throughout the column. The ionized gas is too hot and too dense for a significant deviation to occur as the electrons and ions thermally equilibrate more quickly than any other process.

The profiles for the timescales and energy transfer rates are shown in Figure 9. The three sources exhibit similar behavior in three basic regions, which we generalize here. The electron–ion equilibration timescale (see Equation (140)) is fast enough to ensure that they are in near thermal equilibrium at all times. The corresponding energy transfer rate, ${\dot{U}}_{\mathrm{ei}}$, is orders of magnitude smaller than the four $\dot{U}$ terms shown in Figure 9 and is therefore not included. The three regions are described as follows:

• 1.  
  Region 1 begins at the top of the column and extends over a majority of the column length. The cyclotron cooling timescale initially dominates the bremsstrahlung absorption timescale and is slightly faster than photon Comptonization. The inverse-Compton temperature is initially much larger than the gas temperature, and energy added to the gas by hotter photons is nearly offset by cyclotron emission. These two processes compete against each other as thermal bremsstrahlung absorption continually adds heat. Near the end of Region 1 the surplus heat reservoir in the photons is emptied while thermal bremsstrahlung absorption accelerates, thereby resulting in a role reversal between the photons and electrons at the entrance to Region 2.
• 2.  
  The gas begins to heat the photons via inverse-Compton scattering at the Region 2 entrance. Bulk fluid pressure and density rises rapidly as the accreting material builds-up near the surface. Electron scattering becomes the dominant timescale and bremsstrahlung absorption becomes the dominant heat transfer process. The energy transfer rates and the temperatures exhibit rising exponential behavior inside the extended sinking regime.
• 3.  
  Region 3 begins at the top of the thermal mound for Cen X-3 and LMC X-4, where the parallel absorption optical depth exceeds unity. Thermal bremsstrahlung absorption and Comptonization dominate all other processes. Most photons are either absorbed or scatter within the mound and cannot escape. Mass density and pressure increase as the fluid stagnates near the stellar surface.

Zoom In Zoom Out Reset image size

We conclude that the Comptonization timescale and energy transfer rate determine two markedly different regions in the column. The upper region (Region 1) starts at the free-streaming surface and reaches deep into the column to between 1.5 km and 2.5 km above the stellar surface. It is characterized by slowly changing energy transfer from photons to the gas, where cyclotron and Comptonization processes dominate, and the relative positions among the various timescales are mostly unchanged. Beyond the column half-way point the bulk fluid density begins to play a very important role as the accreting material decelerates, fluid kinetic energy dissipates, and density rises to the point where scattering becomes the dominant timescale. The additional collisions result in accelerated energy addition to the gas due to bremsstrahlung thermal free–free absorption.

The lowest ∼500 m of Region 1 is characterized by slowing energy transfer from the photons to the electrons as they approach thermal equilibrium. The radiation sonic surface is located within a few hundred meters of the thermal mound surface. The dynamics change dramatically at the entrance to Region 2, where the gas temperature exceeds the inverse-Compton temperature. Here, the energy transfer rates are orders of magnitude larger than those in the upper region (Region 1) and the dominant heat mechanisms are bremsstrahlung thermal absorption and inverse-Compton scattering.

The values of the surface magnetic field at the magnetic pole, B_*, obtained in each of our three source models are larger than the corresponding values obtained in the BW07 models. The magnetic field strengths in the BW07 models were set as constants. The field strength in our model, however, quickly drops due to the dipole implementation, where B ∝ r⁻³, thereby resulting in a larger ratio /_cyc at the cyclotron energy _cyc(r) for a photon at energy . The difference in approaches becomes especially relevant when we compare the electron scattering cross-sections with the theoretical values from Equation (31). We examine this more closely in Paper II where the average photon energy, $\bar{\epsilon }$, is computed along the column. Here, we compare our parallel and angle-averaged scattering cross-sections (labeled as WWB) with those used in the BW07 model, which are shown in Tables 13 and 14. We notice our values are an order of magnitude larger. Canuto et al. (1971) showed that the magnitude of the electron scattering cross-section is reduced from the Thomson cross-section when  < _cyc, over all propagation angles θ with respect to the magnetic field direction, according to

$$\begin{eqnarray}&&\displaystyle \frac{\sigma (\theta )}{{\sigma }_{{\rm{T}}}}\propto {\left[\displaystyle \frac{\epsilon }{B(r)}\right]}^{2}.\end{eqnarray} \tag{ 142 }$$

We see from Equation (142) that the ratio σ(θ)/σ_T in our model will be larger than those of BW07 because B(r) is lower with increasing r, whereas the BW07 model maintained a constant B-field throughout.

Table 13.  Parallel Electron Scattering Cross-section

Source	σ∥/σT	σ∥/σT	Ratio
	WWB	BW07	${}^{\mathrm{WWB}}/{}_{\mathrm{BW}07}$
Her X-1	2.60 × 10−3	2.93 × 10−4	8.9
Cen X-3	3.00 × 10−3	4.51 × 10−4	6.7
LMC X-4	2.50 × 10−3	3.98 × 10−4	6.3

Download table as:  ASCIITypeset image

Table 14.  Angle-averaged Electron Scattering Cross-sections

Source	$\bar{\sigma }/{\sigma }_{{\rm{T}}}$	$\bar{\sigma }/{\sigma }_{{\rm{T}}}$	Ratio
	WWB	BW07	${}^{\mathrm{WWB}}/{}_{\mathrm{BW}07}$
Her X-1	1.02 × 10−3	4.15 × 10−5	24.6
Cen X-3	7.51 × 10−4	8.30 × 10−5	9.0
LMC X-4	4.18 × 10−4	4.85 × 10−5	8.6

Download table as:  ASCIITypeset image

The capabilities of our model, in its current state, permit a detailed physical study over the full length of an X-ray pulsar accretion column, including the radial dependences of many key phenomena, such as the photon emission, the hydrodynamic and thermodynamic structure, the imprint of the cyclotron absorption feature, and the stellar surface magnetic field strength. In this paper, we have modeled the dynamical behavior of three high-luminosity sources spanning a range of luminosities from L_X ∼ 10³⁷–10³⁸ erg s⁻¹. We find that the gas pressure does not play a dynamically significant role in these sources. However, this question needs to be reexamined in the context of lower-luminosity sources, such as X Persei, in which the pressure of the gas could play a more significant role, as discussed by Langer & Rappaport (1982). It is also likely that in these sources, the effect of the gas-mediated, discontinuous shock at the base of the column will be more significant than in the high-luminosity pulsars studied here (Becker et al. 2012). We plan to pursue these questions in future work.

Our model provides a new tool for estimating the surface magnetic field strength, and it may also facilitate studies of the variation of the cyclotron absorption energy ${\epsilon }_{\mathrm{cyc}}$ as the luminosity is varied (Mihara et al. 1995; Staubert et al. 2007; Becker et al. 2012). In particular, our model can be used to investigate observed changes in the Her X-1 pulse phase-averaged cyclotron line energy (Staubert et al. 2007) by changing the mass flow rate. More recent observations by Staubert et al. (2014) show a reduction in the line energy over a span of 16 years. We can perform a parameter study using our model to yield possible explanations for this behavior. This may also include experimental modification of the boundary conditions due to stellar surface magnetic field compression or other anomalies such as energy transport into or out of the star (Payne & Melatos 2007; Mukherjee et al. 2013).

We also plan to eventually implement energy-dependent electron scattering cross-sections, relativistic corrections for gravity, a cyclotron absorption term to observe the imprint of the cyclotron resonant scattering feature in the calculated spectrum, the energy-dependent bremsstrahlung thermal free–free absorption coefficient instead of the Rosseland mean coefficient, and a second spatial dimension to observe 2D dynamical behavior. We also intend to couple the five dynamical conservation equations and the photon transport PDE into a single-platform simulation within the COMSOL finite element environment. This may also facilitate studies of the effect of an imbedded gas-mediated shock, thereby allowing us to observe how the electron and ion temperatures react to the presence of a density discontinuity near the base of the accretion column. This would permit a direct comparison with the work of Langer & Rappaport (1982) and Canalle et al. (2005), where the discontinuous shock near the stellar surface has a profound effect on the thermodynamics in the column. The authors would like to acknowledge several useful comments from the anonymous referee that helped to significantly improve the presentation of the material.