Mis understanding of the Heat Method
When solving the equation $(M - tL_C)u = \delta_\gamma$, you effectively have to invert the operator:
$$
u = (M - tL_C)^{-1} \delta_\gamma
$$
Note that while the individual operators $M$ and $L_C$ are only local, containing information about the individual vertices and edges of the mesh, the inverse operator is decidedly not local. Inversion is a global operation that incorporates information across the whole matrix being inverted, which means information across the whole mesh. (Note that $M$ has a trivial inverse since it's diagonal, but $L_C$ definitely doesn't.)
In practice, you might not literally compute the inverse matrix but rather use Gaussian elimination or some such to solve the $u$ for a given $\delta_\gamma$. In that case, it's the solving process that integrates information across the whole mesh. Effectively $L_C$ is encoding a bunch of constraints between neighboring vertices, and the solver has to satisfy all the constraints across the whole mesh at once.
(Note that it's not a trivially parallelizable job. There are methods for using parallelism to accelerate solving large linear systems, but they're going to be multi-pass methods that transfer information up and down between different &quot;scales&quot; of the problem. Think like parallel prefix scan, or FFT—that sort of thing.)
BTW, another perspective on this equation is that it is solving for the fixed point of $u$ in
$$
u = M^{-1}(\delta_\gamma + tL_Cu)
$$
This is yet another way of algebraically rewriting that same equation. In this form it looks a lot like the rendering equation $L = L_e + \int L \, f_{\text{brdf}}$! It has a similar structure, where $\delta_\gamma$ is like the emitted radiance and $L_C$ is like the scattering at surfaces. Just as in the rendering equation, you're looking for a global equilibrium solution. (And if rendering was discretized, we could similarly solve it by turning the scattering operator into a matrix and throwing the whole thing into a linear solver.)