Radiation MHD Simulations of Neutron Star Accretion Columns in X-ray Pulsars and Ultra-Luminous X-ray Sources

The dynamics of magnetized, radiation pressure dominated plasmas has long been a problem of central importance to the brightest objects observed in the universe, from massive stars to quasars. These plasmas are notoriously vulnerable to the formation of strong density inhomogeneities, simply because pressure continuity in this regime only requires continuity of temperature, not density. Moreover, in the presence of a magnetic field and stratification due to a gravitational field, such plasmas are prone to photon bubble instabilities, which drive the creation of such density inhomogeneities. This may allow the plasma to radiate at much higher luminosities than one would expect from photon diffusion through a more homogeneous medium. Nowhere in the universe should these effects be more important than in high luminosity accretion onto neutron stars, both in accretion-powered X-ray pulsars and in the recently discovered pulsating ultraluminous X-ray sources. The latter are emitting (isotropic) luminosities as high as 100 times the critical value at which gravity balances radiation pressure (the so-called Eddington limit). The central problem in such sources is how an accretion column on the surface of the neutron star can sustain such high levels of radiation pressure, and whether such columns can be stably confined by the magnetic field of the neutron star. We propose to conduct a suite of radiation magnetohydrodynamic simulations of these accretion columns, ranging from the relatively low luminosity regime of accreting X-ray pulsars to the high luminosity regime of ultra-luminous X-ray sources. These simulations will be performed with our radiation magnetohydrodynamic code Athena++, which incorporates full angle-dependent, but frequency-integrated, radiation transport. We have already run some preliminary 2D simulations with simplified opacities that exhibit the development of photon bubble instabilities. As part of this proposal, we will extend this into 3D and characterize the emergent spectra and viewing-angle dependent variability of the emergent luminosity. We will also incorporate special relativistic hydrodynamics and the more complex opacities associated with the high magnetic fields of the neutron star. Our results will be useful for interpreting spectra and variability data (particularly X-ray pulse profiles) obtained by NASA high energy astrophysics missions, such as RXTE and NuSTAR. Our simulations will also help answer a broader astrophysics question of whether there truly exists a maximum luminosity for an accretion powered source, and what is its value. Because of the inherent, time-dependent and three dimensional nature of this type of flow, this question can only be credibly answered with numerical simulations aimed at interpreting observational data.