How to understand exactly why gravity darkening happens on rotating stars?
The argument goes something like this.
Hydrostatic equilibrium means that the local pressure gradient is proportional to the local density multiplied by a latitude-dependent local gravity. If the pressure just depends on density and temperature, this means that those quantities will also just depend on latitude and therefore will be constant along an equipotential surface. i.e. pressure, temperature and density are functions of the effective gravitational potential $\phi$.
For stars with radiative outer envelopes, the heat flux is proportional to the temperature gradient, multiplied by some stuff (like inverse opacity) that just depends on density and temperature.
But
$$\nabla T(\phi) = \frac{dT}{d\phi}\nabla \phi = f(\phi)g_{\rm eff}$$
If we now say that at the surface $\sigma T_{\rm eff}^4$ equals the radiative flux, then we recover the Von Zeipel gravity darkening law that $T_{\rm eff}$ is proportional to $g_{\rm eff}^{1/4}$.
The missing step in this argument is to show the $f(\phi)$ is constant. Given that the photosphere is defined as where the optical depth is some fixed value (usually 1 or 2/3), and may be assumed to depend only on temperature and density, then this also lies on an equipotential. But $dT/d\phi$ also only depends on $\phi$ and so must also be constant along an equipotential.
For more detail, although missing the last paragraph above(!), see https://www.astro.umd.edu/~jph/Stellar_Rotation.pdf .
The situation is much more complex for stars with convective envelopes or differential rotation and I think can only be tackled through detailed modelling.