  Numerical simulation of numerous mm/cm-sized particles embedded in a gaseous
disk has become an important tool in the study of planet formation and in
understanding the dust distribution in observed protoplanetary disks. However,
the mutual drag force between the gas and the particles can become so stiff,
particularly because of small particles and/or strong local solid
concentration, that an explicit integration of this system is computationally
formidable. In this work, we consider the integration of the mutual drag force
in a system of Eulerian gas and Lagrangian solid particles. Despite the
entanglement between the gas and the particles under the particle-mesh
construct, we are able to devise a numerical algorithm that effectively
decomposes the globally coupled system of equations for the mutual drag force
and makes it possible to integrate this system on a cell-by-cell basis, which
considerably reduces the computational task required. We use an analytical
solution for the temporal evolution of each cell to relieve the time-step
constraint posed by the mutual drag force as well as to achieve the highest
degree of accuracy. To validate our algorithm, we use an extensive suite of
benchmarks with known solutions in one, two, and three dimensions, including
the linear growth and the nonlinear saturation of the streaming instability. We
demonstrate numerical convergence and satisfactory consistency in all cases.
Our algorithm can for example be applied to model the evolution of the
streaming instability with mm/cm-sized pebbles at high mass loading, which has
important consequences for the formation scenarios of planetesimals.
