We study an ensemble of diverse (inhomogeneous) thermostatically controlled loads aggregated to provide the demand response (DR) services in a district-level energy system. Each load in the ensemble is assumed to be equipped with a random number generator switching heating/cooling on or off with a Poisson rate, r, when the load leaves the comfort zone. Ensemble diversity is modeled through inhomogeneity/disorder in the deterministic dynamics of loads. Approached from the standpoint of statistical physics, the ensemble represents a non-equilibrium system driven away from its natural steady state by the DR. The ability of the ensemble to recover by mixing faster to the steady state after its DR’s use is advantageous. The trade-off between the level of the aggregator’s control, commanding the devices to lower the rate r, and the phase-space-oscillatory deterministic dynamics is analyzed. Then, we study the effect of the load diversity, investigating four different disorder probability distributions (DPDs) ranging from the case of the Gaussian DPD to the case of the uniform with finite support DPD. We show that stronger regularity of the DPD results in faster mixing, which is similar to the Landau damping in plasma physics. Our theoretical analysis is supported by extensive numerical validation.