Implementation of the Shearing Box Approximation in Athena

TL;DR

This paper details the implementation of the shearing box approximation in the Athena MHD code, introducing algorithms for energy conservation, orbital advection, and boundary condition modifications to improve simulation accuracy of accretion disk dynamics.

Contribution

It presents new numerical methods for implementing the shearing box approximation in Athena, including energy-conserving schemes and an orbital advection algorithm for MHD.

Findings

Energy conservation for epicyclic oscillations to round-off error.

Orbital advection reduces truncation error and improves efficiency.

Conservation of magnetic flux and other quantities to better than 0.03% over hundreds of orbits.

Abstract

We describe the implementation of the shearing box approximation for the study of the dynamics of accretion disks in the Athena magnetohydrodynamics (MHD) code. Second-order Crank-Nicholson time differencing is used for the Coriolis and tidal gravity source terms that appear in the momentum equation for accuracy and stability. We show this approach conserves energy for epicyclic oscillations in hydrodynamic flows to round-off error. In the energy equation, the tidal gravity source terms are differenced as the gradient of an effective potential in a way which guarantees that total energy (including the gravitational potential energy) is also conserved to round-off error. We introduce an orbital advection algorithm for MHD based on constrained transport to preserve the divergence-free constraint on the magnetic field. This algorithm removes the orbital velocity from the time step constraint, and makes the truncation error more uniform in radial position. Modifications to the shearing box boundary conditions applied at the radial boundaries are necessary to conserve the total vertical magnetic flux. In principle similar corrections are also required to conserve mass, momentum and energy, however in practice we find the orbital advection method conserves these quantities to better than 0.03% over hundreds of orbits. The algorithms have been applied to studies of the nonlinear regime of the MRI in very wide (up to 32 scale heights) horizontal domains.