We present evidence for a duality between Jackiw-Teitelboim gravity minimally coupled to a free massive scalar field and a single-trace two-matrix model. One matrix is the Hamiltonian H of a holographic disorder-averaged quantum mechanics, while the other matrix is the light operator \cal O dual to the bulk scalar field. The single-boundary observables of interest are thermal correlation functions of \cal O. We study the matching of the genus zero one- and two-boundary expectation values in the matrix model to the disk and cylinder Euclidean path integrals. The non-Gaussian statistics of the matrix elements of \cal O correspond to a generalization of the ETH ansatz.
We describe multiple ways to construct double-scaled matrix models that reproduce the gravitational disk correlators. One method involves imposing an operator equation obeyed by H and \cal O as a constraint on the two matrices. Separately, we design a model that reproduces certain double-scaled SYK correlators that may be scaled once more to obtain the disk correlators.
We show that in any single-trace, two-matrix model, the genus zero two-boundary expectation value, with up to one \cal O insertion on each boundary, can be computed directly from all of the genus zero one-boundary correlators. Applied to the models of interest, we find that these cylinder observables depend on the details of the double-scaling limit. To the extent we have checked, it is possible to reproduce the gravitational double-trumpet, which is UV divergent, from a systematic classification of matrix model `t Hooft diagrams. The UV divergence indicates that the matrix integral saddle of interest is perturbatively unstable. A non-perturbative treatment of the matrix models discussed in this work is left for future investigations.
