We present axisymmetric two-temperature general relativistic radiation magnetohydrodynamic (GRRMHD) simulations of the inner region of the accretion flow onto the supermassive black hole M87. We address uncertainties from previous modeling efforts through inclusion of models for (1) self-consistent dissipative and Coulomb electron heating (2) radiation transport (3) frequency-dependent synchrotron emission, self-absorption, and Compton scattering. We adopt a distance $D=16.7$ Mpc, an observer angle $theta = 20^{circ}$, and consider black hole masses $M/M_{odot} = (3.3times10^{9}, 6.2times10^{9})$ and spins $a_{star} = (0.5, 0.9375)$ in a four-simulation suite. For each $(M, a_{star})$, we identify the accretion rate that recovers the 230 GHz flux from VLBI measurements. We report on disk thermodynamics at these accretion rates ($dot{M}/dot{M}_{mathrm{Edd}} sim 10^{-5}$). The disk remains geometrically thick; cooling does not lead to a thin disk component. While electron heating is dominated by Coulomb rather than dissipation for $r gtrsim 10 GM/c^2$, the accretion disk remains two-temperature. Radiative cooling of electrons is not negligible, especially for $r lesssim 10 GM/c^2$. The Compton $y$ parameter is of order unity. We then compare derived and observed or inferred spectra, mm images, and jet powers. Simulations with $M/M_{odot} = 3.3times10^{9}$ are in conflict with observations. These simulations produce mm images that are too small, while the low-spin simulation also overproduces X-rays. For $M/M_{odot} = 6.2times10^{9}$, both simulations agree with constraints on radio/IR/X-ray fluxes and mm image sizes. Simulation jet power is a factor $10^2-10^3$ below inferred values, a possible consequence of the modest net magnetic flux in our models.
Download