We develop a class of matrix models which implement and formalize the `eigenstate thermalization hypothesis' (ETH) and point out that in general these models must contain non-Gaussian corrections, already in order to correctly capture thermal mean-field theory, or to capture non-trivial OTOCs as well as their higher-order generalizations. We develop the framework of these `ETH matrix models', and put it in the context of recent studies in statistical physics incorporating higher statistical moments into the ETH ansatz. We then use the `ETH matrix model' in order to develop a matrix-integral description of JT gravity coupled to a single scalar field in the bulk. This particular example takes the form of a double-scaled ETH matrix model with non-Gaussian couplings matching disk correlators and the density of states of the gravitational theory. Having defined the model from the disk data, we present evidence that the model correctly captures the JT+matter theory with multiple boundaries, and conjecturally at higher genus. This is a shorter companion paper to the work [1], serving both as a guide to the much more extensive material presented there, as well as developing its underpinning in statistical physics.
