The physics of a magnetically confined, radiation pressure supported column of plasma plays a defining role in understanding the observations of accretion-powered X-ray pulsars, including pulsating ultraluminous X-ray sources (ULXs). Near the neutron star accretor, the accretion flow is constrained by the strong magnetic field to fall along the magnetic field lines. At a sufficiently high accretion rate, the inflow is shocked above the stellar surface and forms a columnar structure below, radiating most of accretion power via sideways emission in a so-called ‘fan-beam’ pattern. The misalignment of the anisotropic radiation emission with respect to the neutron star spin axis results in the observed pulsations. We perform radiative relativistic MHD simulations to study the nonlinear dynamics of the accretion column. The column structure is extremely dynamical and exhibits kHz quasi-periodic oscillations. The existence of the photon bubble instability is identified in simulated accretion columns but proved to be not responsible for triggering the oscillatory behaviors. Instead, the oscillations originate from the inability of the system to resupply heat to locally balance the sideways cooling. The column structure is very sensitive to the shock geometry, which directly determines the cooling efficiency. The time-averaged column structures from the simulations can be approximately reproduced by a 1D stationary model provided one corrects for the actual 2D shape of the time-averaged column. I will also discuss how reduction of scattering opacity by the magnetic field can alter the column structure and variability.