The question addressed in this paper is whether it is possible to model in a fully self-consistent way a thermal plasma ion population by a finite number of test particles of a Monte Carlo model. This Monte Carlo model for the guiding centre drift motion has been developed to study collisional ion transport in an axisymmetric tokamak equilibrium. The model calculates the drift motion of test particles interacting collisionally with ion and electron populations representing the plasma. Momentum conservation in ion–ion collisions is enforced. Calculations with the model for a JET ohmic equilibrium in limiter configuration have shown an 8% particle loss rate during an energy confinement time and a corresponding energy loss rate. To describe steady state plasmas the Monte Carlo model has been extended with a recycling scheme in order to conserve particles, momentum and energy. The recycling of a lost ion as a neutral follows a physics based model while the recycling of the lost momentum and energy is enforced by an ad hoc prescription. The test particles are intended to model the plasma they interact with. To test the self-consistency of this model, test particle profiles of density, temperature, drift velocity, heat flux, etc are accumulated by snapshot techniques during the entire motion of each test particle. Extensive calculations are made on a JET equilibrium with two sets of assumed plasma profiles, peaked and hollow. The results show that the accumulated density and temperature profiles agree approximately with the assumed plasma profiles, demonstrating a degree of self-consistency. A discussion of this result is presented. From the recycling of the lost test particle energy the recycling power function is accumulated. This function is required to maintain a steady state for the test particle population. It agrees with the corresponding function required to maintain a steady state for the plasma ion population, which is obtained from the divergence of the ion heat flux derived from the test particle heat flux. The approximate agreement of these profiles completes the self-consistency of the calculations.