We study multiple scattering of elastic waves with a Monte Carlo method. We take into account the mode conversions and the polarization of the S waves. Some important physical parameters relevant to the description of the polarization are recalled, such as the definition and properties of the elastic Stokes vector. We briefly derive the scattering and Mueller matrices, as well as the differential and total scattering cross sections for one spherical inclusion embedded in a homogeneous matrix. The results of the single-scattering problem are used as a building block for multiple scattering. A Monte Carlo method to simulate the propagation of full elastic waves is presented. We pay a special attention to the convergence toward the diffusive regime which exhibits the equilibration of the P and S energy densities. Our simulations show the shear energy to become very rapidly dominant in the coda and the S to P energy density ratio to tend to 10.4 for a Poisson solid, as predicted by the equipartition theorem. However, the typical timescale and length scale to reach equipartition heavily depend on the scattering parameter kPa, where kP is the P wave number and a is the sphere radius. For Rayleigh scattering (kPa≪1) we find a smooth evolution of energy density with time and a slow convergence toward the equilibration, mainly because of the large difference between the P and S scattering mean free paths in this case. On the other hand, for Rayleigh-Gans scattering (kPa~1.2,1.6) a peak of energy associated with the forward scattered waves is observed, followed by a slow decay according to the diffusion approximation. We find that after only a few mean free times, equipartition is reached in spite of the strong anisotropy of the scattering in this regime. As the scattering parameter kPa increases, we find that equipartition is again delayed because the transport mean free paths become quite large. We find that a large source-station distance favors a rapid equilibration. This effect is seen to be very pronounced for Rayleigh scatterers. When a source of P waves is considered, the equipartition time can be twice as long as compared with a shear source. The time evolution of the EP/ES ratio could be used as a marker for the different scattering mechanisms. ¿ 2000 American Geophysical Union