  A numerical method to implement a linearized Coulomb collision operator in
the two-weight $\delta f$ Monte Carlo method for multi-ion-species neoclassical
transport simulation is developed. The conservation properties and the
adjointness property of the operator in the collisions between two particle
species with different temperatures are verified. The linearized operator in a
$\delta f$ Monte Carlo code is benchmarked with other two kinetic simulations,
a $\delta f$ continuum gyrokinetic code with the same linearized collision
operator and a full-f PIC code with Nanbu collision operator. The benchmark
simulations of the equilibration process of plasma flow and temperature
fluctuation among several particle species show very good agreement between
$\delta f$ Monte Carlo code and the other two codes. An error in the H-theorem
in the two-weight $\delta f$ Monte Carlo method is found, which is caused by
the weight spreading phenomenon inherent in the two-weight $\delta f$ method.
It is demonstrated that the weight averaging method serves to restoring the
H-theorem without causing side effect.
