We describe the algorithm, implementation and numerical tests of a multifluid dust module in the Athena++ magnetohydrodynamic (MHD) code. The module can accommodate an arbitrary number of dust species interacting with the gas via aerodynamic drag (characterized by the stopping time), with a number of numerical solvers. In particular, we describe two second-order accurate, two-stage, fully-implicit solvers that are stable in stiff regimes including short stopping time and high dust mass loading, and they are paired with the second-order explicit van-Leer and Runge-Kutta gas dynamics solvers in Athena++, respectively. Moreover, we formulate a consistent treatment of dust concentration diffusion with dust back-reaction, which incorporates momentum diffusion and ensures Galilean invariance. The new formulation and stiff drag solvers are implemented to be compatible with most existing features of Athena++, including different coordinate systems, mesh refinement, shearing-box and orbital advection. We present a large suite of test problems, including the streaming instability in linear and nonlinear regimes, as well as local and global setting, which demonstrate that the code achieves the desired performance. This module will be particularly useful for studies of dust dynamics and planet formation in protoplanetary disks.
Comment: 32 pages, 14 figures, ApJS resubmitted after addressing referee's comments