Gas and Radiation Pressure Effects on Accretion Column Properties

# Dynamical and Radiative Properties of X-ray Pulsar Accretion Columns: Phase-Averaged Spectra

## Abstract

The availability of the unprecedented spectral resolution provided by modern X-ray observatories is opening up new areas for study involving the coupled formation of the continuum emission and the cyclotron absorption features in accretion-powered X-ray pulsar spectra. Previous research focusing on the dynamics and the associated formation of the observed spectra has largely been confined to the single-fluid model, in which the super-Eddington luminosity inside the column decelerates the flow to rest at the stellar surface, while the dynamical effect of gas pressure is ignored. In a companion paper, we have presented a detailed analysis of the hydrodynamic and thermodynamic structure of the accretion column obtained using a new self-consistent model that includes the effects of both gas and radiation pressure. In this paper, we explore the formation of the associated X-ray spectra using a rigorous photon transport equation that is consistent with the hydrodynamic and thermodynamic structure of the column. We use the new model to obtain phase-averaged spectra and partially-occulted spectra for Her X-1, Cen X-3, and LMC X-4. We also use the new model to constrain the emission geometry, and compare the resulting parameters with those obtained using previously published models. Our model sheds new light on the structure of the column, the relationship between the ionized gas and the photons, the competition between diffusive and advective transport, and the magnitude of the energy-averaged cyclotron scattering cross section.

X-ray pulsar — accretion — accretion columns
1

## 1 Introduction

Nearly a half-century has passed since the first observations of X-ray emission from accretion-powered X-ray pulsars. In these sources, the gravitational potential energy of the accreting gas is efficiently converted into kinetic energy in an accretion column that is collimated by the strong magnetic field (G). Since the surface of the neutron star represents a nearly impenetrable barrier, the kinetic energy is converted into thermal energy at the bottom of the column, and escapes through the column walls in the form of X-rays. Hence the emergent X-ray luminosity, , is essentially equal to the accretion luminosity, , where and denote the stellar mass and radius, respectively, and is the accretion rate. The X-ray emission is generated along the column and in the vicinity of one or both of the magnetic polar caps, thereby forming “hot spots” on the stellar surface. The X-rays were initially thought to emerge as fan-shaped beams near the base of the accretion column, but it soon became clear that a pencil beam component was sometimes necessary in order to obtain adequate agreement with the observed pulse profiles (Tsuruta & Rees 1974; Bisnovatyi-Kogan & Komberg 1976; Tsuruta 1975).

Analysis of X-ray pulsar spectra is usually performed by averaging over many rotation periods to obtain phase-averaged profiles. The resulting spectra are characterized by a power-law continuum extending up to a high-energy exponential (thermal) cutoff at keV. The spectra also frequently display putative cyclotron absorption features (used to estimate the magnetic field strength), and a broadened iron emission line at keV (e.g., White et al. 1983; Coburn et al. 2002).

Until recently, attempts to simulate the spectra of accretion-powered X-ray pulsars using first-principles physical models did not yield good agreement with the observed spectra. The earliest spectral models were based on the emission of blackbody radiation from the hot spots, and these were unable to reproduce the nonthermal power-law continuum. The presence of the cyclotron absorption feature led to the development of more sophisticated models based on a static slab geometry, in which the emitted spectrum is strongly influenced by cyclotron scattering (e.g., Mészáros & Nagel 1985a,b; Nagel 1981; Yahel 1980b). While the magnetized slab models are able to roughly fit the shape of the observed cyclotron absorption features, they remained unable to reproduce the nonthermal power-law X-ray continuum.

The lack of a physics-based model for the formation of X-ray pulsar spectra led to the characterization of the observed spectra using parameters derived from a variety of ad hoc functional forms. Coburn et al. (2002) describe in detail three analytical functions commonly used to empirically model the X-ray continuum. The first function is the power law with high-energy cutoff (PLCUT), given by

 PLCUT(E)=AE−Γ{1 ,(E≤Ecut) ,e−(E−Ecut)/Efold ,(E>Ecut) , (1)

where is the X-ray energy, is the photon spectral index, and and denote the cutoff and folding energies, respectively. The second function (Tanaka 1986) uses the same power law index , combined with a Fermi-Dirac form for the high-energy cutoff (FDCO), defined by

 FDCO(E)=AE−Γ11+e(E−Ecut)/Efold . (2)

The third function (Mihara 1995) uses two power laws, and , in combination with an exponential cutoff (NPEX), written as

 NPEX(E)=(AE−Γ1+BE+Γ2)e−E/Efold . (3)

Although these functional forms have no physical basis, they are often used to characterize the observed spectral shapes.

In later work, it eventually became clear that the observed power-law continuum was the result of a combination of bulk and thermal Comptonization occurring inside the accretion column. The first physically-motivated model based on these principles, that successfully described the shape of the X-ray continuum in accretion-powered pulsars, was developed by Becker & Wolff (2007, hereafter BW07). This new model allowed for the first time the computation of the X-ray spectrum emitted through the walls of the accretion column based on the solution of a fundamental radiation-transport equation. The BW07 model was recently ported to XSPEC using a model described by Wolff et al. (2016). Another variation of the original BW07 model was developed and utilized by Farinelli et al. (2012; 2016, hereafter F16), which added an additional second-order scattering term to the transport equation related to the bulk flow. It is important to note that the Wolff et al. (2016) and Farinelli et al. (2012, 2016) implementations were premised on the utilization of an approximate form for the accretion velocity profile, which was assumed to stagnate at the stellar surface while varying in proportion to the scattering optical depth measured from the stellar surface. This approximate velocity profile stagnates at the stellar surface as required, but it does not follow the exact expected variation with altitude (Becker 1998).

While the BW07 and F16 models have demonstrated notable success in fitting the observed X-ray spectra for higher luminosity sources, these analytical models do not include all of the fundamental hydrodynamic and thermodynamic processes occurring within the accretion column. It is therefore natural to ask how the incorporation of additional physics would impact the ability of the BW07-type models to fit the observed X-ray spectra. This has motivated us to investigate the importance of additional radiative and hydrodynamical physics within the context of a detailed numerical simulation. The new simulation includes the implementation of a realistic dipole geometry, rigorous physical boundary conditions, and self-consistent energy transfer between electrons, ions, and radiation. In particular, we explore the dynamical effect of the gas pressure, which was neglected by BW07 and F16, in order to examine how the radiation and gas pressures combine to determine the dynamical structure of the accretion column and the associated radius-dependent and vertically-integrated, partially-occulted X-ray spectra.

This is the second in a series of two papers in which we describe the new coupled radiative-hydrodynamical model. The integrated approach involves an iteration between an ODE-based hydrodynamical code that determines the dynamical structure, and a PDE-based radiative code that computes the radiation spectrum. The iterative process converges rapidly to yield a self-consistent and unique description of the dynamical structure of the accretion column and the energy distribution in the emergent radiation field. In West et al. (2016, hereafter Paper I), we focus on solving the coupled hydrodynamical conservation equations to describe the underlying gas and radiation hydrodynamics within the accretion column. We employed Mathematica to solve the set of nonlinear hydrodynamical equations in dipole coordinates to determine the accreting bulk fluid velocity, the electron and ion temperature profiles, and the variation of the energy flux.

In this paper (Paper II) we focus on solving the fundamental PDE photon transport equation to describe the production of the observed radius-dependent and phase-averaged X-ray spectra. The bulk velocity and electron temperature profiles developed using the Mathematica calculation described in Paper I are used as input for the analysis conducted here, in which we solve the PDE photon transport equation using the COMSOL multiphysics module. The iterative dual-platform calculation is based on the transfer of information between the COMSOL radiation spectrum calculation and the Mathematica hydrodynamical structure calculation. Using this iterative approach, we are able to model realistic X-ray pulsar accretion columns in dipole geometry by solving the photon transport equation using self-consistent spatial distributions for the accretion velocity profile and the electron temperature profile .

## 2 Hydrodynamical Model Review

Our self-consistent model for the hydrodynamics and the radiative transfer occurring in X-ray pulsar accretion flows is based on a fundamental set of conservation equations governing the flow velocity, , the bulk fluid mass density, , the radiation energy density, , the ion energy density, , the electron energy density, , and the total radial energy transport rate, , where is the radius measured from the center of the neutron star. In our hydrodynamical model, the field-aligned flow is decelerated by the combined pressure of the ions, the electrons, and the radiation. The total pressure is therefore given by the sum

 Ptot=Pe+Pi+Pr , (4)

where , and denote the electron , ion, and radiation pressures, respectively. The equations-of-state for the material components are

 Pi=nikTi ,Pe=nekTe . (5)

We note that in the X-ray pulsar application, the value of cannot be computed using a thermal relation because in general the radiation field is not in thermal equilibrium with either the matter or itself. Hence the radiation pressure must be computed using an appropriate conservation relation. The corresponding energy densities are related to the pressure components via

 Pi=(γi−1)Ui ,Pe=(γe−1)Ue ,Pr=(γr−1)Ur , (6)

where (see Paper I) we set and . The ratio of specific heats for the radiation field is given by .

As shown in Paper I, the mathematical model can be reduced to a set of five coupled, first-order, nonlinear ordinary differential equations satisfied by , , and the ion, electron, and radiation sound speeds, , , and , respectively, defined by

 a2i=γiPiρ ,a2e=γePeρ ,a2r=γrPrρ . (7)

The ion and electron temperatures, and , respectively, are related to the corresponding sound speeds via

 a2i=γikTimtot ,    a2e=γekTemtot , (8)

where and denote the ion and electron masses, respectively, and . We assume here that the accreting gas is composed of pure, fully-ionized hydrogen.

### 2.1 Conservation Equations

Here, we briefly review the hydrodynamical conservation equations and the underlying model assumptions. The fundamental independent variables in our model are the photon energy, , and the radius, , measured from the center of the star. Since our model includes only one explicit spatial dimension, the structure across the column at a given value of is not treated in detail. Therefore, all of the physical quantities (velocity, temperature, pressure, flux, and density) developed here are interpreted as mean values across the column (see Appendix A).

The mass continuity equation can be written in dipole geometry as (e.g., Langer & Rappaport 1982)

 ∂ρ(r)∂t=−1A(r)∂∂r[A(r)ρ(r)v(r)] , (9)

where denotes the radial inflow velocity, and the cross-sectional area of the column, , is related to the solid angle subtended by the dipole column at radius , by

 A(r)=r2Ω(r)=r3Ω∗R∗ . (10)

The solid angle subtended by the column at the stellar surface is given by

 Ω∗=2π(cosθ1−cosθ2) , (11)

where and denote the inner and outer polar angles of the accretion column at the stellar surface, respectively. When , the column is completely filled; otherwise, the central portion of the column is empty. We therefore find that the solid angle varies with radius according to

 Ω(r)=(rR∗)Ω∗ . (12)

In a steady-state situation, the mass accretion rate is conserved, and is related to the density and velocity via

 ˙M=Ωr2ρ|v| . (13)

Combining Equations (12) and (13), we find that the variation of the mass density is given by

 ρ(r)=˙MR∗Ω∗r3|v| . (14)

Under the assumption that the accreting gas is composed of pure, fully-ionized hydrogen, the electron and ion number densities are equal and can be computed using

 ne=ni=ρmtot=˙MR∗mtotΩ∗r3|v| , (15)

where .

The total energy flux in the radial direction, averaged over the column cross-section at radius , can be written as

 F(r)=12ρv3+v(Pi+Pe+Pr+Ui+Ue+Ur)−cneσ∥∂Pr∂r−GM∗ρvr , (16)

where denotes the electron-scattering cross section for photons propagating parallel to the magnetic field. The flux is defined to be negative for energy flow in the downward direction, towards the stellar surface, and the corresponding accretion velocity is negative (). The terms on the right-hand side of Equation (16) represent the kinetic energy flux, the enthalpy flux, the radiation diffusion flux, and the gravitational energy flux, respectively. The diffusion of radiation energy through the walls is given by and is modeled using an escape-probability formalism. The value of the energy transport rate varies as a function of the radius in response to the escape of radiation perpendicular to the magnetic field direction. We therefore utilize a total energy conservation equation of the form

 ∂∂t(12ρv2+Ui+Ue+Ur−GM∗ρr)=−1A(r)∂˙E∂r−Urtesc , (17)

where is the total energy transport rate in the radial direction, defined by

 ˙E≡A(r)F(r) , (18)

and the mean escape timescale, , is discussed in Section 3 and in Appendix A.

It is convenient to work in terms of non-dimensional radius, flow velocity, sound speed, and total energy transport rate variables by introducing the quantities

 ~r=rRg ,  ~v=vc ,  ~ai=aic ,  ~ae=aec ,  ~ar=arc ,  ~E=˙E˙Mc2 , (19)

where is the gravitational radius, defined by

 Rg≡GM∗c2 . (20)

The computational domain extends from the top of the accretion column, at radius , down to the stellar surface at dimensionless radius , which is based on canonical values for the stellar mass, , and the stellar radius, km.

The set of five fundamental hydrodynamical differential equations that must be solved simultaneously using Mathematica comprises Equations (21), (22), (23), (24), and (25), which were described in detail in Paper I and are shown here as

 d~ard~r =~ar2(3~r+1~vd~vd~r)−σ∥Rg2mtotc˙MAγr~ar(~E+~v22+~a2eγe−1+~a2iγi−1+~a2rγr−1−1~r), (21) d~aid~r =(1−γi)~ai2(3~r+1~vd~vd~r+Rgc2A˙Mγi˙Ui~a2i), (22) d~aed~r =(1−γe)~ae2(3~r+1~vd~vd~r+Rgc2A˙Mγe˙Ue~a2e), (23) d~Ed~r =~a2rc(ℓ2−ℓ1)γr(γr−1)~v(R3∗Rg~r3)1/2min(c,cτ⊥), (24) d~vd~r =~v~v2−~a2i−~a2e{3(~a2i+~a2e)~r−1~r2+σ∥Rgmtotc˙MA(~E+~v22+~a2iγi−1+~a2eγe−1 +~a2rγr−1−1~r)+Rgc2A˙M[(γi−1)˙Ui+(γe−1)˙Ue]} , (25)

where the cross-sectional area of the dipole accretion column is given by Equation (10). The ions do not radiate appreciably, and therefore the ion energy density is only affected by adiabatic compression and Coulomb energy exchange with the electrons (see Langer & Rappaport 1982). Conversely, the electrons experience free-free, cyclotron, and Compton heating and cooling, in addition to adiabatic compression and the Coulomb exchange of energy with the ions. The radiation energy density responds to both the creation and absorption of photons via free-free and cyclotron processes, in addition to Compton scattering with the electrons.

The thermal coupling terms in Equations (22), (23), and (25) describe a comprehensive set of electron and ion heating and cooling processes which are broken down as follows,

 ˙Ue = ˙Uemitbrem+˙Uabsbrem+˙Uemitcyc+˙Uabscyc+˙UComp+˙Uei , ˙Ui = −˙Uei . (26)

The terms in the expression for denote, respectively, thermal bremsstrahlung (free-free) emission and absorption, cyclotron emission and absorption, the Compton exchange of energy between the electrons and photons, and electron-ion Coulomb energy exchange. In our sign convention, a heating term is positive and a cooling term is negative. The mathematical descriptions of each term in Equation (26) are provided in Paper I. We note that the term plays a fundamental role in the energy exchange between photons and electrons, as discussed in detail below and in Section 3.

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. The electrons in the accretion column experience heating due to free-free absorption of low-frequency radiation, and they also experience heating and cooling due to the absorption and emission of thermal cyclotron radiation. These emission and absorption processes can play an important role in regulating the temperature of the gas. On average, however, cyclotron absorption does not result in net heating of the gas, due to the subsequently rapid radiative de-excitation, and we therefore set in our dynamical calculations. Near the surface of the accretion column, however, photons scattered out of the outwardly directed beam are not replaced, and this leads to the formation of the observed cyclotron absorption feature in a process that is very analogous to the formation of absorption lines in the solar spectrum (Ventura et al. 1979). While cyclotron absorption does not result in net heating of the gas, due to the rapid radiative de-excitation, cyclotron emission will cool the gas. The electrons can also be heated or cooled via Coulomb collisions with the protons, depending on whether the electron temperature exceeds the ion temperature .

Compton scattering plays a fundamental role in the formation of the emergent X-ray spectrum. It is critically important in establishing the radial variation of the electron temperature profile through the exchange of energy between photons and electrons. The Compton cooling rate, prior to dimensionless units conversion, is expressed in terms of the mass density, , the electron sound speed, , and the radiation sound speed, . The complete derivation of the electron cooling rate due to Compton scattering is given in Section 3.4, and the final result is given by (see Equation (76))

 ˙UComp(r)=ne¯σcmec2 4kTe[g(r)−1]Ur(r) , (27)

where the temperature ratio function, , is defined by

 g(r)≡TIC(r)Te(r) . (28)

The derivation of is provided later in Section 3.4. Equation (27) can be stated in terms of the sound speeds as

 ˙UComp=4¯σ(g−1)mecγeγr(γr−1) ρ2a2ra2e . (29)

The sign of depends on the value of . The electrons experience Compton cooling if (i.e. ); otherwise, the electrons are heated via inverse-Compton scattering.

Our task is to solve the five coupled hydrodynamic conservation equations, shown in Equations (21)-(25) above, to determine the radial profiles of the dynamic variables , and , subject to the boundary conditions discussed in Section 2.2 and Section 3.1. Once these profiles are available, the electron temperature can be computed from the electron sound speed using the relation (see Equation (8))

 Te=mtotc2γek~a2e . (30)

The solutions for and are the fundamental output of the solved system of conservation equations. These are used as input to the COMSOL finite-element environment in order to compute the photon distribution function, , inside the column, which is the focus of this paper.

Solving the five coupled conservation equations to determine the radial profiles of the quantities , , , , and requires an iterative approach, because the rate of Compton energy exchange between the photons and the electrons depends on the relationship between the electron temperature, , and the inverse-Compton temperature, , which is determined by the shape of the radiation distribution (see Equation (77)). In order to achieve a self-consistent solution for all of the flow variables, while taking into account the feedback loop between the dynamical calculation and the radiative transfer calculation, the simulation must iterate through a specific sequence of steps. The steps required in a single iteration are: (1) the computation of the initial estimate of the accretion column dynamical structure, which is obtained by solving the five conservation equations described in Paper I under the assumption that ; (2) calculation of the associated radiation distribution function by solving the radiative transfer equation described here in Paper II; (3) computation of the inverse-Compton temperature profile based on integration of the radiation energy distribution using Equation (77); and (4) re-computation of the dynamical structure using the new estimate for . The iteration continues until adequate convergence is achieved between successive temperature profiles for all values of . The iterative process is discussed in detail in Paper I, and also in Section 4.1 of this paper.

### 2.2 Hydrodynamic Boundary Conditions

Here we summarize the set of boundary conditions utilized in the Mathematica integration of the set of hydrodynamical conservation equations, required in order to determine the structure of the accretion column. At the top of the accretion column (radius ), the inflow velocity equals the local free-fall velocity, so that

 vtop≡vff(rtop)=−(2GM∗rtop)1/2 . (31)

We also assume that the local acceleration of the gas is equal to the gravitational value, and therefore

 dvdr∣∣∣r=rtop=(GM∗2r3top)1/2 . (32)

At the upper surface of the dipole-shaped accretion funnel, the radiation escapes freely, forming a “pencil-beam” component in the observed radiation field. This physical principle can be used to derive a useful boundary condition at the upper surface of the accretion column. The radial component of the radiation energy flux, averaged over the cross-section of the column at radius , is given in the diffusion approximation by (see Equation (71))

 Fr(r)=−c3neσ∥dUrdr+43vUr , (33)

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 towards the stellar surface (with ). The fact that the top of the accretion column is the last scattering surface implies that the photon transport makes a transition from diffusion to free streaming at , and therefore we make the following replacement in Equation (33) at the upper boundary,

 −c3neσ∥dUrdr→cUr ,r→rtop . (34)

We incorporate Equation (34) into Equation (33) to conclude that the radiation energy flux at the upper surface is given by

 Fr(r)∣∣∣r=rtop=(c+43vtop)Ur(rtop) , (35)

where is the free-fall velocity given by Equation (31).

The ionized gas flows downward after entering the top of the accretion funnel at radius , 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 sum of the diffusion and the advection components 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 we conclude the radiation energy flux as . We refer to this as the “mirror” surface boundary condition, which can be written as

 Fr(r)∣∣∣r=R∗=0 . (36)

We approximate stagnation at the stellar surface in our simulations using the condition

 limr→R∗|v(r)|<∼ 0.01c . (37)

## 3 Photon Transport Equation

The equation describing photon transport is a second order, non-linear, partial differential equation (PDE) of elliptic type, which can be solved numerically with appropriate boundary conditions using the finite-element method (FEM). The time-dependent transport equation for an isotropic photon distribution function (phase-space density of photons) experiencing Compton scattering, spatial diffusion, and bulk advection in an X-ray pulsar accretion column is given by the generalized Kompaneets (1957) equation (e.g., Becker 1992; Becker & Wolff 2007)

 ∂f∂t=ne¯σcmec21ϵ2∂∂ϵ[ϵ4(f+kTe∂f∂ϵ)]−→∇⋅→Fph+˙fprod+˙fabs−13ϵ2∂∂ϵ(ϵ3→v⋅→∇f) , (38)

where the specific flux, or “streaming,” of the photons in the radial direction is given by (Gleeson & Axford 1967; Skilling 1975)

 →Fph=−κ→∇f−ϵ→v3∂f∂ϵ , (39)

and represents the spatial diffusion coefficient. The terms on the right-hand side of Equation (38) describe, respectively, the transport of the photons through the energy space due to stochastic electron scattering (the Kompaneets operator); the spatial diffusion of radiation; the injection of seed photons; the absorption of radiation by the accreting gas; and the differential work performed on the photons due to the convergence of the accreting gas (Gleeson & Webb 1978; Cowsik & Lee 1982). It is important to note that Equation (38) must be supplemented by suitable radiation boundary conditions for the photon distribution function , imposed at the column walls, the stellar surface, and at the top of the column. This is further discussed in Section 3.1.

Equations (38) and (39) can be combined to obtain the single equivalent equation (e.g., Becker 2003)

 ∂f∂t+→v⋅→∇f=ne¯σcmec21ϵ2∂∂ϵ[ϵ4(f+kTe∂f∂ϵ)]+→∇⋅(κ→∇f) +13(→∇⋅→v)ϵ∂f∂ϵ+˙fprod+˙fabs . (40)

In order to solve this equation, we must adopt a specific geometrical model. The accretion flow in an X-ray pulsar is channeled by the strong magnetic field, and therefore we will employ a dipolar geometry for the accretion column, in which the cross-sectional area of the column, , is given by Equation (10). The angle-dependent operators in Equation (40) can be removed by averaging across the column cross-section, and by assuming azimuthal symmetry around the magnetic field axis. The resulting transport equation, satisfied by the isotropic photon distribution function can be written in the dipole geometry as

 ∂f∂t+v∂f∂r=ne¯σcmec21ϵ2∂∂ϵ[ϵ4(f+kTe∂f∂ϵ)]+1A(r)∂∂r[A(r)κ∂f∂r] +13A(r)∂[A(r)v]∂rϵ∂f∂ϵ+˙fprod+˙fabs+˙fesc , (41)

where the spatial diffusion coefficient in the radial direction is given by

 κ(r)=c3ne(r)σ∥ . (42)

We demonstrate in Appendix A that the photon escape term can be implemented using an escape probability formalism by writing

 ˙fesc(r,ϵ)=−f(r,ϵ)tesc(r) , (43)

where the mean escape time at radius , denoted by , is given by

 tesc(r)=ℓesc(r)w⊥(r) . (44)

Here, denotes the perpendicular escape distance across the column at radius , computed using

 ℓesc(r)=(ℓ2−ℓ1)(rR∗)3/2 , (45)

and is the perpendicular diffusion velocity, given by

 w⊥(r)=min[c,cτ⊥(r)] ,τ⊥(r)=ne(r)σ⊥ℓesc(r) , (46)

where denotes the electron-scattering cross section for photons propagating perpendicular to the magnetic field, and is the corresponding perpendicular optical thickness of the accretion column at radius . The parameters and in Equation (45) represent the inner and outer arc-length surface radii at the base of the accretion column, respectively (see Equation (11)). When , the column is completely filled; otherwise, the center of the column is empty, which may reflect the manner in which the gas is entrained onto the field lines at large radii, or alternatively, it may reflect the possible presence of a quadrupole field component (see Section 6.4). Note that if , then transitions to the free-streaming value, , as required.

Significant variability of the X-ray emission on pulse-to-pulse timescales has been observed from some X-ray pulsars, including Her X-1 (Becker et al. 2012; Staubert et al. 2007; Staubert et al. 2014), in which the pulsation period is about one second. The pulsation period is several orders of magnitude larger than the dynamical timescale for accretion onto the star, and therefore, to first approximation, the accretion flow can be considered to be steady-state in the frame of the star. It is therefore sufficient for our purposes here to focus on solving the time-independent version of Equation (41), although the steady-state assumption may be relaxed in future work. In order to reduce the complexity of the simulations while retaining the essential physical details, we will assume in this paper that all of the physical quantities represent averages across the column cross-section at a given value of the spherical radius , measured from the center of the star. In future work, we plan to relax this assumption and treat the structure of the column using a more realistic two-dimensional model, but that is beyond the scope of the present paper. In the steady-state situation of interest here, Equation (41) can be written in the flux-conservation form

 ∂∂r[ϵ2r3(−κ∂f∂r−vϵ3∂f∂ϵ)]+∂∂ϵ{−ϵ2r3[ne¯σcmec2ϵ2(f+kTe∂f∂ϵ)−ϵv3∂f∂r]} =ϵ2r3(˙fprod+˙fabs−ftesc) , (47)

where we have used Equations (10) and (43) to substitute for the cross-sectional area and the escape term , respectively.

The solution space of the radiation transport equation contains spatial and energy dimensions, forming a 2D grid within the FEM-based COMSOL multiphysics environment. In the radial dimension, the computational domain extends from  km (the stellar surface) up to  km (the top of the accretion column). The precise value for is determined self-consistently for each source, as discussed in Paper I. In the energy dimension, the computational domain extends from a minimum photon energy  keV to a maximum photon energy  keV. The final form of the transport equation is found by converting Equation (47) to the dimensionless radius and velocity variables and (see Equation (19)). The result obtained is

 ∂∂~r[ϵ2R2g~r3(−c3neσ∥Rg∂f∂~r−13c~vϵ∂f∂ϵ)] +∂∂ϵ{R3g~r3ϵ2[−ne¯σcmec2ϵ2(f+kTe∂f∂ϵ)+ϵc~v3Rg∂f∂~r]} =ϵ2R3g~r3(˙fprod+˙fabs−ftesc) , (48)

where we have also utilized Equation (42). Since COMSOL employs the finite-element method, it is well-suited to the iterative solution method required to solve the photon transport Equation (48), which is a second order, elliptic PDE. The associated boundary conditions are discussed in detail below, along with the forms utilized for the terms describing radiation absorption, injection, and escape.

### 3.1 Photon Transport Boundary Conditions

In order to solve the partial differential transport equation (Equation (48)) for the photon distribution function in the COMSOL finite-element environment, we are obliged to develop and apply a suitable set of physical boundary conditions in radius and energy. The upper surface of the accretion column, located at radius , represents the last-scattering surface for photons diffusing up from deeper layers in the column. Hence at this radius, the photon transport makes a transition from classical diffusion to free-streaming, as occurs in the scattering layer above the photosphere in a stellar atmosphere. In the dipole geometry considered here, the specific flux or “streaming” in the radial direction is given by (cf. Equation (39))

 Fph,r≡−κ∂f∂r−ϵv3∂f∂ϵ . (49)

Using Equation (42) to substitute for the spatial diffusion coefficient , we find that the transition from diffusion to free-streaming implies that (cf. Equation (34))

 −c3neσ∥dfdr→cf ,r→rtop , (50)

and therefore

 Fph,r→cf−ϵv3∂f∂ϵ ,r→rtop . (51)

This relation expresses the spatial boundary condition applied to the photon distribution function at the column top.

We must also develop and apply an appropriate boundary condition at the stellar surface, where . In our idealized, one-dimensional model, the density of the accreting gas formally diverges as it settles onto the star, and therefore no radiation flux can penetrate the stellar surface. This concept is physically manifest in the “mirror” boundary condition, which states that the radiation streaming in the radial direction must vanish as . The mathematical implementation of this condition is given by

 Fph,r→0 ,r→R∗ . (52)

The mirror condition described above provides the spatial boundary condition that must be satisfied by at the stellar surface.

The lowest and highest energies treated in our simulations are denoted by keV and keV, respectively. The boundary conditions applied at these two energies are determined by the requirement that the radiation transport rate along the energy axis must vanish in the limit , and also in the limit . The numerical values we choose for and are such that we can assume that the radiation transport rate vanishes at these two energies, which implies that no photons cross the boundaries of the computational domain at and . The high-energy boundary condition can be written as (cf. Equation (47))

 ne¯σcmec2ϵ2(f+kTe∂f∂ϵ)−ϵv3∂f∂r=0 ,ϵ=ϵmax , (53)

and likewise, the low-energy boundary condition is given by

 ne¯σcmec2ϵ2(f+kTe∂f∂ϵ)−ϵv3∂f∂r=0 ,ϵ=ϵmin . (54)

Taken together, Equations (51), (52), (53), and (54) comprise the set of four boundary conditions that we apply on the four edges of the computational domain employed in the COMSOL environment.

### 3.2 Photon Sources

The injection of seed photons due to the various emission mechanisms operating inside the accretion column is represented by the term in Equation (48). Following BW07, we can express the photon production rate by writing

 ˙fprod=Qprodr2Ω(r)=Qbrem+Qcyc+Qbbr2Ω(r) , (55)

where the source functions , , and correspond to bremsstrahlung, cyclotron, and blackbody emission, respectively. The source functions are normalized so that gives the number of photons injected per unit time between radii and , with energy between and . The source function is therefore related to the volume emissivity, , via

 ϵ2Q(r,ϵ)dϵdr=r2Ω(r)˙nϵ(r,ϵ)dϵdr , (56)

where gives the number of photons injected per unit time per unit volume in the energy range between and .

Computation of the thermal bremsstrahlung (free-free) source function, , is based upon Equation (5.14b) from Rybicki & Lightman (1979), which gives for the free-free volume emissivity (in cgs units)

 ˙nffϵ=3.7×1036ρ2T−1/2eϵ−1e−ϵ/kTe . (57)

We can combine Equations (56) and (57) to obtain for the bremsstrahlung source function

 Qbrem=3.7×1036Ωr2ρ2T−1/2eϵ−3e−ϵ/kTe . (58)

The cyclotron volume emissivity, , gives the production rate of cyclotron photons per unit volume per unit energy. For the case of pure, fully-ionized hydrogen, we can employ Equations (7) and (11) from Arons et al. (1987) and substitute into Equation (56) to obtain

 ˙ncycϵ=2.1×1036ρ2B−3/212H(ϵcyckTe)e−ϵcyc/kTeδ(ϵ−ϵcyc) , (59)

where keV denotes the cyclotron energy, , and is a piecewise function defined by

 H(x)≡{0.15√7.5 ,x≥7.5 ,0.15√x ,x<7.5 . (60)

The radial variation of the magnetic field is evaluated using Equation (4) from Paper I, with , which gives

 B(r)=B∗(rR∗)−3 , (61)

where denotes the field strength at the center of the magnetic pole on the stellar surface. The cyclotron source function, , is found by substituting Equation (59) into Equation (56), which gives

 Qcyc=2.1×1036Ωr2ρ2B−3/212H(ϵcyckTe)ϵ−2e−ϵcyc/kTeδ(ϵ−ϵcyc) . (62)

In our numerical simulations, the -function in Equation (62) is approximated by a narrow Gaussian, centered on the cyclotron energy, with a standard deviation of 1 keV.

The blackbody source function, , is defined by writing

 Qbb=S(ϵ)δ(r−rth) , (63)

where is the radius at the upper surface of the thermal mound, and the function is related to the Planck distribution according to

 ϵ3S(ϵ)dϵ=r2thΩ(rth)πBϵ(ϵ)dϵ , (64)

with denoting the blackbody intensity, given by

 Bϵ(ϵ)=2ϵ3c2h31eϵ/kTth−1 . (65)

We can combine Equations (63), (64), and (65) to obtain for the blackbody photon source function the result

 Qbb(r,ϵ)=2πr2thΩ(rth)c2h3δ(r−rth)eϵ/kTth−1 . (66)

The radial -function in Equation (66) is approximated using a narrow Gaussian feature, centered on the thermal mound, with a standard deviation of km.

### 3.3 Photon Absorption

The photon absorption term, , in Equation (48) can be written in the form

 ˙fabs=−ftabs , (67)

where the mean absorption time, , describes the average time for a photon to be reabsorbed through the bremsstrahlung absorption process (free-free absorption), given by

 1tabs=cαR . (68)

Here, denotes the Rosseland mean of the absorption coefficient for fully ionized hydrogen, expressed in cgs units by (Rybicki & Lightman 1979)

 αR=1.7×10−25T−7/2en2e  ∝  cm−1 . (69)

The requirement of self-consistency in the treatment of the hydrodynamics and the radiative transfer is a major goal of the model developed here and in Paper I, and therefore it is essential that we utilize the (frequency-independent) Rosseland mean absorption coefficient in Equation (68).

### 3.4 Transport Phenomena and Photon-Electron Energetics

The angle-averaged (isotropic) photon distribution function, is normalized so that gives the number of photons per unit volume with energy between and inside the accretion column at radius , measured from the center of the star. From knowledge of , we can compute the radiation number density and energy density using the relations

 nr(r)=∫∞0ϵ2f(r,ϵ)dϵ ,     Ur(r)=∫∞0ϵ3f(r,ϵ)dϵ . (70)

We derive the ordinary differential equations satisfied by and in Appendix B by integrating the full transport equation using Equations (70). The radiation energy flux, (erg cm s), is defined by

 Fr(r)=−c3neσ∥dUrdr+43vUr , (71)

which is obtained by integrating the radial component of the specific flux, , via (see Equation (39))

 Fr(r)=∫∞0ϵ3Fph(r,ϵ)dϵ . (72)

Compton scattering plays a dominant role in determining the nature of the thermal balance between the electrons and the radiation field. In a single scattering, the mean value of the energy transferred from the photon to the electron is given by (Rybicki & Lightman 1979)

 ⟨Δϵ⟩=ϵmec2(4kTe−ϵ) , (73)

where is the photon energy. The term in Equation (26) is the rate of change of the electron energy density due to Compton scattering. Hence this quantity can be computed by writing

 ˙UComp(r)=−ne¯σc∫∞0ϵ2f(r,ϵ)⟨Δϵ⟩dϵ , (74)

where the negative sign appears because the electrons gain energy whenever for the photons. Combining Equations (73) and (74) yields

 ˙UComp(r)=ne¯σcmec2[∫∞0ϵ4f(r,ϵ)dϵ−4kTe∫∞0ϵ3f(r,ϵ)dϵ] , (75)

which can be written in the equivalent form

 ˙UComp(r)=ne¯σcmec2 4k</