We propose a material point method (MPM) to model the evolving multi-body contacts due to crack growth and fragmentation of thermo-elastic bodies. By representing particle interface with an implicit function, we adopt the gradient partition techniques introduced by Homel and Herbold (2017) to identify the separation between a pair of distinct material surfaces. This treatment allows us to replicate the frictional heating of the evolving interfaces and predict the energy dissipation more precisely in the fragmentation process. By storing the temperature at material points, the resultant MPM model captures the thermal advection–diffusion in a Lagrangian frame during the fragmentation, which in return affects the structural heating and dissipation across the frictional interfaces. The resultant model is capable of replicating the crack growth and fragmentation without requiring dynamic adaptation of data structures or insertion of interface elements. A staggered algorithm is adopted to integrate the displacement and temperature sequentially. Numerical experiments are employed to validate the diffusion between the thermal contact, the multi-body contact interactions and demonstrate how these thermo-mechanical processes affect the path-dependent behaviors of the multi-body systems.