This paper is concerned with the ensemble statistics of the energy density of a random system subjected to a harmonic point load. Both the mean and variance of the reverberant energy density ͑i.e., the total response minus the direct field͒ are investigated. It is shown that the ensemble average of the reverberant energy density is not spatially homogeneous, but rather the value at the drive point is between two and three times the spatially averaged value, depending on the modal overlap factor. This result is closely analogous to established results regarding the case of transient excitation. Expressions are also derived for the relative variance of the reverberant energy density both remote from the drive point and at the drive point, and a number of anomalies are found in existing results. A comparison is made with simulation results for a randomized plate, and this comparison highlights the importance of the ensemble size considered in the simulations. The present analysis is based on a random point process model of the system natural frequencies, and both Gaussian orthogonal ensemble ͑GOE͒ statistics and Poisson statistics are considered.