A local time-stepping algorithm has been developed to improve the numerical ef ciency of Lagrangian particlebased Monte Carlo methods for obtaining the steady-state solutions of the modeled probability density function equations of turbulent reacting ows. On each step in the pseudo-time-marching algorithm, the properties of each particle are advanced by a time step, the magnitude of which depends on the particle's spatial location. This algorithm has been incorporated into the consistent hybrid nite volume/particle method. The performance of the local time-stepping method is evaluated in terms of numerical ef ciency and accuracy through application to a nonreacting bluff-body ow. For this test case, it is found that local time stepping can accelerate the global convergence of the hybrid method by as much as an order of magnitude, depending on the grid stretching. Additionally, local time stepping is found to improve signi cantly the robustness of the hybrid method mainly due to the accelerated convection of error waves out of the computational domain. The method is very simple to implement, and the small increase in CPU time per step (typically 3% %) is a negligible penalty compared to the substantial reduction in the number of time steps required to reach convergence. NomenclatureC u = Courant number based on particle velocity C Ä = Courant number based on turbulent frequency D b = bluff-body diameter F = mass-density function for uniform time stepping O F = mass-density function for the local time stepping Q f = mass-weighted probability density function M 2 = total number of grid cells m = particle mass q = mean particle-mass density R j = jet radius r = radial distance s = a global time variable U = particle velocity V = sample space variable for U W = isotropic Wiener process X = particle position x = axial distance x = sample space variable for X ± = Dirac delta function ¹ = numerical weight » x = grid-stretching factor in x direction » y = grid-stretching factor in y direction Ä = conditional mean turbulent frequency h i = mean eld 0 = rms uctuating quantity Q = mass-weighted eld