  A major uncertainty in the structure and dynamics of magnetized, radiation
pressure dominated neutron star accretion columns in X-ray pulsars and
pulsating ultraluminous X-ray sources is that they are thought to be subject to
the photon bubble instability. We present the results of two dimensional
radiation relativistic magnetohydrodynamic simulations of a non-accreting,
static atmosphere to study the development of this instability assuming
isotropic Thomson scattering in the slow diffusion regime that is relevant to
neutron star accretion columns. Photon bubbles generally grow faster toward
shorter wavelengths, until a maximum growth rate is achieved at the radiation
viscosity length scale, which is generally quite small and requires high
numerical resolution to simulate. We confirm the consistency between our
simulation results and linear theory in detail, and show that the nonlinear
evolution inevitably leads to collapse of the atmosphere with the higher
resolution simulation collapsing faster due to the presence of shorter length
scale nonlinear structures. At least in static atmospheres with horizontally
periodic boundary conditions, this resolution dependence may make simulations
of the nonlinear dynamics of photon bubble instability in neutron star
accretion columns challenging. It remains to be seen whether these difficulties
will persist upon inclusion of an accretion flow through the top and
magnetically-confined horizontal boundaries through which photons can escape.
Our results here provide a foundation for such future work.
