  The Athena MHD code has been extended to integrates the motion of particles
coupled with the gas via aerodynamic drag, in order to study the dynamics of
gas and solids in protoplanetary disks and the formation of planetesimals. Our
particle-gas hybrid scheme is based on a second order predictor-corrector
method. Careful treatment of the momentum feedback on the gas guarantees exact
conservation. The hybrid scheme is stable and convergent in most regimes
relevant to protoplanetary disks. We describe a semi-implicit integrator
generalized from the leap-frog approach. In the absence of drag force, it
preserves the geometric properties of a particle orbit. We also present a
fully-implicit integrator that is unconditionally stable for all regimes of
particle-gas coupling. Using our hybrid code, we study the numerical
convergence of the non-linear saturated state of the streaming instability. We
find that gas flow properties are well converged with modest grid resolution
(128 cells per pressure length \eta r for dimensionless stopping time
tau_s=0.1), and equal number of particles and grid cells. On the other hand,
particle clumping properties converge only at higher resolutions, and finer
resolution leads to stronger clumping before convergence is reached. Finally,
we find that measurement of particle transport properties resulted from the
streaming instability may be subject to error of about 20%.
