Theoretical Analysis of the RX J0209.6-7427 X-ray Spectrum During Giant Outburst
Brent F. West Peter A. Becker Georgios Vasilopoulos




0009-0007-1219-7508]Brent F. West
Department of Physics, United States Naval Academy, Annapolis, MD 21402, USA

0000-0002-3718-1293]Peter A. Becker
Department of Physics and Astronomy, George Mason University, Fairfax, VA 20030, USA

0000-0003-3902-3915]Georgios Vasilopoulos
Department of Physics, National and Kapodistrian University of Athens, University Campus Zografos, GR 15784, Athens, Greece
Institute of Accelerating Systems & Applications, University Campus Zografos, Athens, Greece



We model the spectral formation occurring in the binary X-ray pulsar RX J0209.6-7427 during the 2019 super-Eddington outburst. Using a theoretical model previously developed by the authors, we are able to produce spectra that closely resemble the phase-averaged X-ray spectra observed using NuSTAR and Insight-HXMT during low and high luminosity states of the outburst, respectively. The theoretical model simulates the accretion of fully ionized gas in a dipole magnetic field, and includes a complete description of the radiation hydrodynamics, matter distribution, and spectral formation. Type II X-ray outbursts provide an opportunity to study accretion over a large range of luminosities for the same neutron star. The analysis performed here represents the first time both the outburst low and high states of an accretion-powered X-ray pulsar are modeled using a physics-based model rather than standard phenomenological fitting with arbitrary mathematical functions. We find the outer polar cap radius remains constant and the column is more fully-filled with increasing luminosity, Comptonized bremsstrahlung dominates the formation of the phase-averaged X-ray spectrum, and a negative correlation exists between cyclotron centroid energy and luminosity, as expected. The super-Eddington nature of the outburst is rendered possible due to the low scattering cross section for photons propagating parallel to the magnetic field. We also find emission through the column top dominates in both the low and high states, implying the pulse profiles should have a roughly sinusoidal shape, which agrees with observed properties of ultra-luminous X-ray pulsars. 




























§ INTRODUCTION
 

Accretion-powered X-ray pulsars (XRPs) are intriguing systems that enable the study of complex radiative processes in the presence of strong magnetic fields <cit.>. 
The accreting material is entrained from the surrounding accretion disk by the strong magnetic field, and then channeled onto a magnetic pole of the neutron star (NS), forming an accretion column (AC). The luminosity of these systems range from ∼10^34 up to ∼38, and even above the Eddington limit in extreme cases. The classical spherical Eddington limit, L_ Edd=1.25× 10^38M/M_⊙  erg s^-1, is the maximum luminosity a star of mass M can have, in the Thomson scattering limit, and not eject the accreting material by radiation pressure, where M_⊙ is the mass of the sun.

Be X-ray binaries are especially interesting <cit.>. The accretion disk from the companion Be star in these systems periodically provides enhanced mass transfer onto the NS. The variable mass accretion rate converts gravitational potential energy into a highly luminous release of radiation, designated as a Type II, or giant, X-ray outburst <cit.>. The X-ray luminosity can exceed ∼ 10^39  erg s^-1; hence, Type II outbursts from BeXRBs provide a unique window into the behavior of the closest ultra-luminous X-ray sources <cit.> and establish an important connection between typical XRPs and the ULX pulsars (ULXPs).

Application of theoretical models to XRPs in the super-Eddington regime represents an important step towards investigations of ULX systems. The high luminosities of some types of accretion-powered XRPs, and particularly the ULXPs, suggest radiation pressure plays a leading role in the accreting gas dynamics. Some sort of beaming may be required to preserve the structure of the surrounding accretion disk <cit.>. The luminosities of the ULXPs exceed L_ Edd, implying the accreting gas is unstable to ejection via the pressure of the intense radiation field <cit.>. The classical Eddington limit assumes full Thomson scattering, however, which is not applicable in XRPs because the electron scattering cross section for photons propagating parallel to the magnetic field is several orders of magnitude smaller than the Thomson value <cit.>.

Type II outbursts offer a unique laboratory to perform a detailed analysis of the physical processes occurring in the AC over a large range of  L_ X. For low L_ X, the accretion is sub-critical, meaning that L_ X≲ L_c, where L_c is the critical luminosity above which radiation pressure becomes dynamically important <cit.>. The AC is believed to transition to the super-critical regime at the peaks of the most extreme Type-II outbursts, with luminosity L_ X≳ L_c <cit.>. This can significantly alter both the dynamical structure of the AC, as well as the dominant mode of spectral formation <cit.>.

Physical simulation of the AC requires consideration of a variety of physical processes occurring in the magnetized, accreting plasma, as it passes through either a gas-dominated discontinuous shock (in low-luminosity sources) or a continuous (smooth) radiation-dominated shock (in high-luminosity sources such as the Type II outburst studied here). Seed photons produced via bremsstrahlung, cyclotron, and blackbody emission are reprocessed via bulk and thermal Comptonization before escaping through the walls and top of the column. The first model to include all of these physical features was developed by <cit.>, who applied their model to simulate the formation of the observed X-ray spectra for the sources Her X-1, Cen X-3, and LMC X-4.

The model of <cit.>, and also the similar model of <cit.>, proved successful at reproducing the observed X-ray pulsar spectra in super-critical sources, where radiation pressure plays a dominant role in controlling the dynamics of the accreting material. Additional progress was made by <cit.>, who developed a generalized conical geometry model with application to X-ray pulsars over a wide range of luminosities, and incorporates a more realistic bulk fluid velocity profile. However, these models do not treat the hydrodynamical structure of the column in detail; instead, they assume the velocity profile of the accreting gas is specified (known) in advance.

To date, the most self-consistent model for spectral formation in accretion-powered X-ray pulsars is the one developed by <cit.>, which is based on an iterative method that yields simultaneous, coupled solutions for the radiation-hydrodynamical structure of the accretion column, along with the radiation distribution escaping through the walls and top of the column. We focus on the application of this model to simulate the Type II outburst observed from (hereafter ) that occurred in 2019 <cit.>.

The model of <cit.> is a radiation-hydrodynamics simulation cast in a dipole geometry in which the bulk flow velocity is computed by solving a one-dimensional hydrodynamical equation coupled with ion and electron energy equations that fully account for the energy exchange processes between ions, electrons, and photons. It was developed primarily for the study of luminous pulsars. In this paper, we focus on the application of the  model to the analysis and interpretation of X-ray data obtained using NuSTAR and Insight-HXMT during the outburst.


The strong magnetic fields in accretion-powered XRPs cause the cyclotron energy levels of the electrons to be quantized; therefore, the electron scattering cross section is a complex function of the photon energy and propagation direction <cit.>.
The electron scattering cross section is also qualitatively different and depends on whether the photon is polarized in the ordinary or extraordinary mode. The extraordinary mode photons, in particular, have resonant interactions with the electrons which cause a sharp increase in the scattering cross section at the cyclotron centroid energy, E_ cyc, given by the standard relation


    E_ cyc(r) = 11.57 n B_12(r)  keV ,


where n corresponds to the Landau level from which the radiative de-excitation to the ground state occurs <cit.>, and the dipole magnetic field strength varies with the radius, r, according to


    B_12(r)=B_*/10^12  G(r/R_*)^-3 .


Here the stellar radius and the surface magnetic field strength at the polar cap are denoted by R_* and B_*, respectively. 

The interaction of the extraordinary mode photons with the electrons leads to the formation of a cyclotron resonant scattering feature (CRSF) at photon energy E_ cyc, which is often observed as an absorption-like feature, or simultaneous multiple features, in the spectra of XRPs. They are not true absorption features, but rather represent the scattering of radiation out of the beam due to collisions with electrons in the accreting plasma. The implications of cyclotron scattering for the formation of the dual CRSFs during the giant outburst of are discussed below.



§ OVERVIEW OF THE WWB17 MODEL
 

We investigated the radiation-hydrodynamics inside the AC of using the iterative model of WWB17, which is briefly described here. The complete process to build a self-consistent radiation and hydrodynamical pulsar AC model is described in detail in <cit.>, in which the accretion column aligns with the magnetic field dipole shape, and all of the physical quantities are functions of the radial coordinate, r, measured from the center of the NS, meaning they are averaged across the column cross-section at a given value of r.

The range of the dipole magnetic field radial coordinates provides the computational domain for solving a set of five coupled ordinary differential equations, which describe the spatial evolution of the bulk flow velocity, the total energy flux, and the electron, ion, and radiation sound speeds; therefore, the ion and electron temperatures are determined using separate energy equations. Coulomb collisions in high-luminosity sources rapidly equilibrate the ion and electron temperatures, but they are less efficient in low-luminosity sources, and the two temperatures may be quite different. Simultaneous solutions of the five coupled differential equations yield a complete description of the dynamical and thermal structure inside the AC.

The bulk flow velocity and electron temperature solutions are subsequently used to solve a second-order elliptic partial differential photon transport equation to obtain the radiation distribution inside the AC, which includes terms describing advection, diffusion, free-free (bremsstrahlung) emission and absorption, cyclotron emission and absorption, bulk and thermal Comptonization, and the escape of radiation through the column walls and top. Blackbody seed photons from the surface of the thermal mound are also included in the model. 

An upper boundary condition is imposed at the column top where photon transport makes a transition from diffusion to free streaming at the last scattering surface, and a zero total energy flux boundary condition is imposed at the stellar surface. The scattering of the photons by the electrons in the accreting plasma is modeled using the cross sections σ_|| and σ_⊥ for photons diffusing parallel or perpendicular to the magnetic field lines, respectively. The model also implements an angle-averaged cross section, σ̅, which is used to describe the Compton scattering interaction between the photons and electrons.

The solution for the radiation distribution inside the AC is used to compute the radial profile of the inverse-Compton temperature, T_IC, which plays a critical role in controlling the thermal balance between the gas and the radiation field. The complete set of six conservation equations, comprising five hydrodynamical relations and one radiation transport equation, must be solved iteratively because the photon transport equation requires hydrodynamical information (electron temperature and bulk velocity profiles), whereas the hydrodynamical solutions require knowledge of the inverse-Compton temperature profile. The iterative process is considered converged when both the electron temperature and inverse-Compton temperature profiles each change by less than 1% from the previous iteration. The result is a fully self-consistent model for the hydrodynamical and radiative structure of the AC.



§ ASTROPHYSICAL APPLICATION
 

We used the  model to analyze the formation of the phase-averaged X-ray spectra within 2 epochs during the giant outburst observed from in 2019. Data were collected by and Insight-HXMT during the low and high states, respectively, for both epochs as the NS was accreting at super-Eddington rates.





 §.§ 2019 Outburst of RX J0209: Data Collection


is a BeXRB pulsar (9.3 s period) located in the outer wing of the Small Magellanic Cloud (SMC). During its 2019 major outburst its observed L_ X exceeded 39 <cit.>. The outburst was monitored by all-sky X-ray instruments such as /GBM and /BAT, and pointed observations were made with . 
An X-ray light curve of the outburst that covers the brightest phase is plotted in Figure <ref>, where L_ X is estimated from GBM pulsed fluxes (12 – 50 keV) at a distance of D = 55kpc <cit.>.

The Nuclear Spectroscopic Telescope Array () mission is the first focusing high-energy X-ray telescope operating in the 3 – 79 keV energy range <cit.>. We analysed the DDT observation of (obsid: 90502352002) following <cit.>.
The Hard X-ray Modulation Telescope <cit.> collected data during 5 visits (proposal id: P0214059) over an ∼5 day period around the peak of the outburst <cit.>.
is equipped with the Low Energy (LE, 0.7 – 15 keV), Medium Energy (ME, 5 – 30 keV), and High Energy (HE, 20 – 250 keV) X-ray telescopes. The data were downloaded from the archive and analysed with hxmtsoft v2.04 using default options.

Based on the collected statistics, we limited our study to the following energy ranges for each detector; 2.0 – 8 keV for LE, 8 – 30 keV for ME; and 30 – 100 keV for HE <cit.>. We validated our analysis by comparing with results in the literature <cit.>. Specifically, we fitted the same spectral models and we compared the model parameters and residuals obtained with the previous results. Finally, we used the unfolded spectra in the  modeling process. 



 §.§ Observed Luminosity and Mass Flow Rate


The bolometric X-ray luminosity, L_ X, is related to the mass accretion rate, Ṁ, via the standard formula


    L_ X = G M_* Ṁ/R_* ,


where we assume all gravitational potential energy is converted into radiation, which is a reasonable approximation given the extreme density of the NS crust. We do not have complete knowledge of the bolometric luminosity, but we can obtain a reasonable approximation by using the energy bands corresponding to the specific detectors used to make the X-ray observations. For the low state of the RX J0209 outburst observed in 2019 using NICER and NuSTAR, the bolometric X-ray luminosity, L_ X, in the energy range 0.5 – 70.0 keV, was found to be L_ X=5.54 × 10^38  erg s^-1 <cit.>. For the high state of the 2019 outburst, observed using Insight-HXMT, the bolometric luminosity in the energy range 1.0 – 150.0 keV was found to be L_ X = 1.10 × 10^39  erg s^-1 <cit.>. The luminosity values are computed assuming isotropic emission with a source distance D = 55kpc.



 §.§ Cyclotron Resonance Scattering Feature (CRSF)


The formation of a CRSF in the observed X-ray spectrum via resonant cyclotron scattering is a challenging technical problem and beyond the scope of the present paper <cit.>. The location of the imprint radius provides valuable information about the strength of the magnetic field, and remains an active research topic. The hydrodynamic equations of the WWB17 model are based upon a dipole magnetic field geometry, and converging the model with a predetermined imprint radius has fundamental effects on the column height and structure.

The supercritical luminosities for both the low and high states of the RX J0209 outburst suggest the CRSF behavior will follow theory initially proposed by <cit.> in which a radiation-dominated shock forms at a radius scaling approximately linearly with the observed luminosity. We expect the imprint radius somewhere higher than the first photon scattering surface, which corresponds to a photon travelling one optical depth in the downward (parallel) direction from the top of the column. We also anticipate a negative correlation between the cyclotron centroid energy, E_cyc, and the source luminosity, L_ X <cit.>. The shock radius in the high state should be at a higher altitude (meaning E_cyc should be a lower value) when compared with the low state.

<cit.> theorized CRSF formation is concentrated along the column where the ratio of the extraordinary photon scattering length to the absorption length is at a minimum value. He found the relationship between the resonant scattering mean free path, ℓ_ scatt, and the absorption mean free path, ℓ_ abs, evaluated at the centroid energy, E_ cyc, for fully-ionized hydrogen gas is given by


    ℓ_scatt/ℓ_abs=3π/2αħ^3 c^4 n_e/E^3_cyc√(m_e/k_ B T_e) ,


where α is the fine structure constant, ħ is the reduced Planck constant, k_ B is Boltzmann's constant, c is the speed of light, and m_e, n_e, and T_e denote the electron mass, number density, and temperature, respectively. The ratio in Equation (<ref>) is a function of the radius, r, via the quantities E_ cyc and T_e, and the radius corresponding to the minimum value of this ratio is interpreted as the cyclotron imprint radius, r_ cyc.

More recent results by <cit.> place the CRSF within the radiative shock, at a radius higher in the column than the location suggested by <cit.>. <cit.> implemented a cylindrical geometry Monte Carlo radiative transfer code and demonstrated prominent CRSFs in the emergent spectrum are generated in the shock. They also noted the possibility of evenly spaced cyclotron lines due to resonant scattering in multiple harmonics. We follow the approach of <cit.> here.

The WWB17 model formation of the CRSF absorption-like feature in the X-ray spectrum is described in an approximate manner by multiplying the escaping phase-averaged spectrum by a 2D Gaussian distribution, A_ cyc(r,ϵ), which depends on the radius, r, and the photon energy, ϵ, given by

    A_ cyc(r,ϵ) = 1-d_cr/2 πσ_ cycσ_r
    e^-(ϵ - E_ cyc)^2/2 σ_ cyc^2
    e^-(r - r_ cyc)^2/2 σ_r^2 ,

where d_cr represents a combined strength parameter, and the standard deviations in the energy and radial directions are denoted by σ_cyc and σ_r, respectively. We note the cyclotron energy, E_ cyc, is a function of r via Equation (<ref>).



 §.§ Free Parameter Selection and Model Convergence


Table <ref> lists the values for the twelve free parameters used to model the outburst low and high states. It also includes values for three fixed quantities, three observational quantities, and several post-processing diagnostic quantities. 

We fix the NS mass and radius quantities using the canonical values M_* = 1.4 M_⊙ and R_* = 10km, respectively. The third fixed quantity is established by setting the perpendicular scattering cross-section equal to the Thomson limit, σ_⊥=σ__ T <cit.>. The previously published bolometric values of L_ X (at source distance D = 55kpc) are used to calculate the corresponding values of the mass accretion rate, Ṁ, using Equation (<ref>). Hence, we classify D, Ṁ, and L_ X as the three observational quantities. 

The next step in the WWB17 modeling process, after establishing the three fixed quantities and the three observational quantities, is to select the twelve free parameters. The first six free parameters are associated with the calculation of the continuum X-ray spectrum and the dynamical structure of the AC. These comprise (1) the electron scattering cross section for photons propagating parallel to the local magnetic field direction, relative to the Thomson limit, σ_||/σ__ T, (2) the angle-averaged electron scattering cross section, relative to the Thomson limit, σ̅/σ__ T, (3) the inner polar cap radius at the stellar surface, ℓ_1, (4) the outer polar cap radius at the stellar surface, ℓ_2, (5) the incident radiation Mach number at the top of the column, ℳ_ r0, and (6) the surface magnetic field strength, B_*, at the base of the column.

The final six free parameters are used to modify the spectrum during post-processing calculations, which occur after the radiative hydrodynamic equations are solved and temperatures are converged. There are five free parameters associated with dual CRSFs located at the fundamental and first harmonic cyclotron frequencies (n=1 and n=2, respectively), and one free parameter associated with interstellar absorption (N_ H).

Separate values for the free parameters in Table <ref> were obtained for the low and high states, respectively, with the exception of the surface magnetic field strength, B_*, which is a free parameter, but was constrained to have the same value in both states, as physically required. See <cit.> for a full discussion.

Cutaway views of the corresponding column structures are depicted in Figure <ref>. The profiles of the primary hydrodynamic and radiative variables in the high state are plotted in Figure <ref>. Similar results were obtained in the low state; therefore, the low state profiles are not included here. The plots include the relevant timescales and energy transfer rates for the electron heating (solid curves) and cooling (dotted curves) processes.

Finding phase-averaged spectra solutions with the WWB17 model requires connecting two software platforms in a manual fashion. The hydrodynamic equations are solved using Mathematica[https://www.wolfram.com/mathematica/], and the photon transport equation is solved using COMSOL[https://www.comsol.com/]. An automated process, such as the one implemented in XSPEC, would require creation of a tabulated model. Each iteration of a model takes 10 – 15 minutes; hence, this approach would require years of computation for a reasonably dense grid when considering all twelve free parameters. Therefore, automated fitting in XSPEC is not practical at this time. The model parameters were simply varied manually. The resulting theoretical spectra were then folded through the detector response matrix using XSPEC to make a quantitative comparison with the observed data in the detector plane. Using this method, we were able to quantitatively compare the WWB17 model with the observed data, yielding values for the reduced χ^2, as well as plots of the residuals. We find that χ^2_ red. = 1.17 in the low state and χ^2_ red. = 1.08 in the high state. These results are included in Figure <ref>.








lll
0pt
Model values for the low and high states.

Free Parameters     Low State     High State
Parallel cross-section σ_∥/σ__ T     4.44×10^-4     3.88 × 10^-4 

Angle-averaged cross-section σ̅/σ__ T     16.0×10^-4     9.50 × 10^-4 

Outer polar cap radius ℓ_2 (m)     640     640 

Inner polar cap radius ℓ_1 (m)     440     20 

Incident radiation Mach number ℳ_r0     1.66     1.40 

Surface magnetic field B_* (10^12G)     7.0     7.0 


CRSF radial σ_r (km)     4.1     4.1 


CRSF 1 energy σ_ cyc (keV)     11.0     8.25 


CRSF 1 strength d_cr (km keV)      95.0      52.0 


CRSF 2 energy σ_ cyc (keV)     16.0     16.5 


CRSF 2 strength d_cr ( km keV)     195.0      202.0 




N_ H (10^21cm^-2)     1.58     1.58 

Fixed Quantities     Low State     High State 

Perpendicular cross-section σ_⊥/σ__ T     1.0     1.0 

M_* (M_⊙)     1.4     1.4 

R_* (km)     10     10 

Observational Quantities†     Low State     High State 

L_ X (erg s^-1)     5.54×10^38     1.10×10^39 

Ṁ (g s^-1)     2.98 × 10^18     5.91 × 10^18 

Distance D (kpc)     55     55 

Diagnostic Quantities‡     Low State     High State 


L_ X^ tot (erg s^-1)     5.06×10^38     1.03×10^39 

L_ X^ brem (erg s^-1)     5.02×10^38     1.02×10^39 

L_ X^ cyc (erg s^-1)     2.52×10^35     4.32×10^34 

L_ X^ bb (erg s^-1)     3.46×10^36     7.59×10^36 

Sonic point radius (km)     14.91     15.34 

CRSF imprinting radius r_ cyc (km)     14.91     15.34 

CRSF 1 energy E_ cyc (keV)     24.44     22.45 

CRSF 2 energy E_ cyc (keV)     48.88     44.90 


Radius of max wall emission (km)     14.91     16.16 

First (parallel) scattering surface (km)     13.77     13.38 






L_ top (erg s^-1)     4.48×10^38     1.02×10^39 

L_ wall (erg s^-1)     5.78×10^37     1.13×10^37 

Thermal mound T_th (K)     2.90×10^8     3.00×10^8 



Surface impact velocity υ_*/c     8.07×10^-3     9.63×10^-3 

Incident ion Mach number ℳ_i0     203.0     183.6 


†The value of Ṁ is computed from L_ X using Equation (<ref>), see the text for the discussion of the bolometric luminosity L_ X.
‡Diagnostic luminosities are computed in the standard energy range 3-79 keV energy for both states for purposes of comparison.






§ DISCUSSION
 

The good agreement between the WWB17 theoretical model and the observed data obtained in the low and high states of the 2019 giant outburst allows us to make quantitative comparisons of the model parameters and physical profiles in the two states. Overall, the radiation hydrodynamics, energy exchange rates, and timescales are very similar in the two states, but there are some interesting differences.

Figure <ref> depicts the radiative hydrodynamic profiles along the full length of the AC for the high state converged model solution. The five panels in descending order depict the bulk velocity, temperatures, timescales, energy coupling rates, and pressure and pressure gradients. The timescales and energy coupling rates associated with gas heating processes are indicated by solid curves, whereas gas cooling processes are indicated by dotted curves. We segment the 8.2 km length of the column into four regions, which are labeled in the 4^ th panel of Figure <ref> (energy coupling rates). The regions separate hydrodynamic milestones as the bulk fluid enters the top of the column and eventually stagnates at the stellar surface and merges with the NS crust.

Region 1 begins where the accreting matter enters the column top with a velocity of 0.48 c. The high gas density ensures the electrons and ions equilibrate to the same temperature over the entire length of the column. Cyclotron cooling dominates all energy coupling rates in Region 1 on a timescale shorter than any heating rate, and therefore the gas temperature does not rapidly equilibrate (higher in the column) to the inverse-Compton temperature of the radiation field, T_ IC. The radius of maximum wall emission is at 16.2 km as photons continue transferring energy to the electrons via Compton scattering while approaching the radiation sonic point.

The plasma enters Region 2 after passing through the radiation sonic point while continuing to decelerate. The thermal Comptonization gas heating rate is a maximum at the top of Region 2, and starts dropping rapidly as photons and electrons approach the same temperature at the bottom of Region 2. The free-free absorption rate quickly rises and becomes the dominant gas heating mechanism. Electrons and photons finally reach thermal equilibrium at the inverse-Compton temperature, T_ IC, which marks the end of Region 2. The Compton scattering process transitions into inverse-Compton scattering downstream from this radius, as the accreting gas enters Region 3, where much of the remaining electron energy is transferred to the photons.

Free-free absorption continues to dominate in Region 3 as the electron temperature increases rapidly. Radiation pressure continues to decelerate the material as it transitions into the sinking regime discussed by <cit.>. Region 4 begins at the top of the thermal mound, where the parallel absorption optical depth, τ^abs_∥, equals unity. The electrons and protons decelerate to rest inside the thermal mound and merge with the star's crust.

We follow <cit.> and set the CRSF imprint radius equal to the radiation sonic point radius, and we allow for the possibility of dual CRSFs appearing at both the fundamental cyclotron frequency (n=1), as well as the first harmonic frequency (n=2). Multiple harmonics are reported in the literature for other sources. For example, <cit.> describe a model with a set of five harmonics for the 1999 giant outburst of 4U 0115+63. In our simulations, the centroid energy for the fundamental frequency was calculated to be 24.44 keV and 22.45 keV for the low and high states, respectively. The centroid energy for the first harmonic was set to twice the fundamental value, corresponding to energies 48.88 keV and 44.90 keV, respectively.

One of the most important free parameters is the surface magnetic field strength at the polar cap, for which we obtained B_* = 7.0 ×10^12 G, a value that can be tested against observational constraints. The best observational estimates of the surface B-field strength rely on two methods. These include the application of torque models that can explain the spin evolution of the NS during the 2019 outburst, and the identification of critical transitions in either the Hardness Ratios (HR) or the pulse profile shape that are believed to be related to changes in column morphology.

Torque modeling of the data, assuming a distance of D = 55kpc, reveals a magnetic field strength of 2 – 6×10^12G at the NS poles, with uncertainties mainly due to the adopted bolometric correction and the details of the specific torque model <cit.>. On the other hand, critical transitions in the pulse profile shape <cit.> and HRs <cit.> suggest a higher B-field strength of ∼2×10^13 G. The differences between the two estimates could be explained by the presence of a multipolar B-field <cit.>.

An additional important result of our analysis is an improved understanding of how the structure of the accretion column evolves with increasing luminosity during the outburst. The outer column radii in both states is found to be 640 m. The unchanged outer polar cap radius is expected because the inner edge of the radiation-dominated accretion disk should not vary significantly, if at all, with the increased mass flow rate <cit.>. The pressure balance between the disk and magnetic field in this inner region remains essentially constant. However, the column is significantly more hollow in the low state (440 m inner polar cap radius), as compared with the high state (20 m inner polar cap radius), which is nearly completely filled. See Figure <ref>.

More puzzling is the transition from a mostly hollow column in the low state to a nearly filled column in the high state. The inner polar cap radius obtained in the low state is ℓ_1 = 440m, whereas in the high state we obtained  ℓ_1 = 20m. This question can't be fully resolved without incorporating a detailed model for the structure of the extended pulsar magnetosphere and the surrounding accretion disk, including the binary orbit and the properties of the companion star, which is beyond the scope of this paper.

We believe that the changes in the inner polar cap radius, ℓ_1, between the low and high states are due to the super-Eddington accretion rates involved here. The high luminosities during the outburst will undoubtedly cause some heating of the accretion disk, which may cause it to puff up <cit.>. This change in the accretion disk geometry will cause a large variation in the radius where the matter is first picked up and entrained onto the magnetic field. In particular, the outer-most pick-up radius is likely to move outwards as the disk puffs up in the high state. Keeping in mind that the outer-most pick-up radius corresponds to the inner wall of the accretion column where it connects with the stellar surface, we conclude an increase in the accretion rate will tend to produce a more filled column in the high state, as we obtained here.



§ CONCLUSION


In this study we modeled the spectral formation occurring in the binary X-ray pulsar RX J0209.6-7427 during the 2019 super-Eddington Type II outburst. We focused on analysis of the X-ray spectra detected during low and high states, between which the accretion rate changed by a factor of ∼ 2. We employed the previously published self-consistent WWB17 model for the radiative and dynamical processes occurring in the AC of the highly magnetized NS, and we find that the simulated spectra closely match the observed spectra during the low and high states, observed with NuSTAR and Insight-HXMT, respectively.

We review our primary findings here and focus on a comparison of the free parameters and the radiative and hydrodynamic profiles inside the AC. This reveals a number of interesting changes as the outburst transitions from the low state to the high state:


    
  * Bremsstrahlung photons provide the primary source of seed photons contributing to the phase-averaged spectrum. The cyclotron and blackbody seed photons make a negligible contribution. The luminosities of the various respective reprocessed seed photon components are indicated in Table <ref>. It is clear that bremsstrahlung provides over 99% of the power in the observed spectra for both the low and high states.
    
    
  * The luminosity radiated from the column top is significantly higher than the wall luminosity for both states, and the relative contribution of the top component increases with luminosity (see Figure <ref> and Table <ref>). This has immediate implications for the shape of the resulting pulse profile. For example, we predict the pulsar profiles of the ULXPs should be roughly sinusoidal, as observed in many cases <cit.>.

    
    
    
  * The angle-averaged and parallel electron scattering cross sections, σ̅ and σ_||, respectively, are larger in the low state. The enhanced cross sections in the low state are due to the stronger cyclotron absorption parameter, d_cr, at the fundamental frequency (n=1) which amplifies the effect of the cyclotron resonance, hence increasing the scattering cross sections (see Table <ref>). The reduced cross sections in the high state cause a larger fraction of the photons to be beamed along the magnetic field lines, leading to a larger column top luminosity. Our findings are consistent with the results of <cit.> and <cit.>.
    
    
  * The CRSF energy, E_ cyc, is 2 keV lower and the radiation sonic point is 0.43 km higher in the high state when compared with the column structure in the low state. This expected behavior is in agreement with the ideas first proposed by <cit.>. The theoretical model developed by <cit.> predicts a negative correlation between centroid energy and emission altitude according to
    
    E_ cyc/E_*=(R_*+h/R_*)^-3 ,

    where E_* is the surface value for the cyclotron energy according to Equation (<ref>). The h parameter corresponds to the emission altitude at which photon diffusion through the column walls is a maximum. Setting the value of h equal to the CRSF imprint altitude (which is also the sonic point altitude) establishes that Equation (<ref>) is essentially a combination of Equations (<ref>) and (<ref>).
    
    
  * The super-Eddington luminosities in Type II outbursts do not disrupt the accretion column structure <cit.> because the scattering cross section for photons propagating parallel to the magnetic field is sub-Thomson by several orders of magnitude <cit.>. There can be an indirect effect on the AC structure, however, if the super-Eddington flux heats the surrounding accretion disk, thereby causing changes in the accretion stream feeding into the column. <cit.>.




Our results have important implications for the application of AC radiative models to the analysis of major outbursts from BeXRBs and ULXPs. In particular, ULXPs share spectral similarities with BeXRBs <cit.>, which have so far mainly been explored using empirical and phenomenological models. Next-generation X-ray telescopes, such as the High Energy X-ray Probe <cit.>, will have higher effective areas extending over larger energy ranges than current detectors, which will facilitate much more detailed studies of super-Eddington bursts <cit.>.


We are grateful to the anonymous referee, who made several important comments and suggestions that led to substantial improvements in the manuscript. GV acknowledges support by H.F.R.I. through the project ASTRAPE (Project ID 7802). The authors would also like to thank Dr. Kenneth Wolfram for his insights and constructive conversations.



aasjournal

