I suffered over this proof quite a bit in my pre-arithmetic youth and wrote up
a number of summaries. A condensed and extremely superficial version is given here, mostly to indulge my own sense of nostalgia.
If by chance someone finds it at all useful, of course I will be delighted. I apologize that I don't say anything about physical
intuition (because I have none), and for repeating parts of the previous nice answers.
It's been years since I've thought about these matters, so I will forgo
all attempts at even a semblance of analytic rigor.
The usual setting has a compact Riemannian manifold $M$, two hermitian bundles $E^+$ and $E^-$, and a linear
operator
$$P:H^+\rightarrow H^-,$$
where $H^{\pm}:=L^2(E^{\pm})$.
With suitable assumptions (ellipticity), $ker(P)$ and $coker(P)$
have finite dimension, and the number of interest is the index:
$$Ind(P)=dim(ker(P))-dim(coker(P)).$$
This can also be expressed as
$$dim(ker(P))-dim(ker(P^{*})),$$ where
$$P^{*}:H^-\rightarrow H^+$$
is the Hilbert space adjoint. A straightforward generalization of the Hodge theorem allows us also to write this in terms of Laplacians
$\Delta^+=P^* P$ and $\Delta^-=PP^*$
as
$$dim(ker(\Delta^+))-dim(ker(\Delta^-)).$$
Things get a bit more tricky when we try to identify the index with the expression ('supertrace,' so-called)
$$Tr(e^{-t\Delta^+})-Tr(e^{-t\Delta^-}).$$
The operator
$$e^{-t\Delta^{\pm}}:H^{\pm}\rightarrow H^{\pm}$$
sends a section $f$ to the solution of Heat equation
$$\frac{\partial}{\partial t} F(t,x)+\Delta^{\pm}F(t,x)=0$$
($x$ denoting a point of $M$) at time $t$ with intial condition $F(0,x)=f(x).$
One important part of this is that there are discrete Hilbert direct sum decompositions
$$H^+=\oplus_{\lambda} H^+(\lambda)$$
and $$H^-=\oplus_{\mu} H^-(\mu)$$
in terms of finite-dimensional eigenspaces for the Laplacians. And then, the identities
$$\Delta^-P=PP^{*}P=P\Delta^+$$
and
$$\Delta^+P^{*}=P^{*}PP^{*}=P^{*}\Delta^-$$
show that the (supersymmetry) operators $P$ and $P^{*}$ can be used to define isomorphisms between all non-zero eigenspaces of the two Laplacians with
a correspondence for eigenvalues as well.
Thus, once you believe that the exponential operators are trace class,
it's easy to see that the only contributions to the trace are from the kernels of the plus and minus Laplacians.
But there, the heat evolution operators are clearly the identity, allowing us to identify the supertrace with the index.
To summarize up to here, we have
$$Ind(P)=Tr(e^{-t\Delta^+})-Tr(e^{-t\Delta^-}).$$
This identity also makes it obvious that the supertrace is in fact independent of $t>0$.
The proofs under discussion all have to do with identifying this supertrace in terms of local expressions that
relate naturally to characteristic classes. The beginning of this process involves first writing the operator
$e^{-t\Delta^+}$ in terms of an integral kernel
$$K^+_t(x,y)=\sum_i e^{-t\lambda_i } \phi^+_i(x)\otimes \phi^+_i(y)$$
where the $\phi^+_i$ make up an orthonormal basis of eigenvectors for the Laplacian.
That is,
$$e^{-t\Delta^+}f=\int_M K^+_t(x,y)f(y)dvol(y)=\sum_i e^{-t\lambda_i } \int_M \phi^+_i(x) \langle \phi^+_i(y),f(y)\rangle dvol(y).$$
Formally, this identity is obvious, and the real work consists of the global analysis necessary to justify the formal computation.
Obviously, there is a parallel discussion for $\Delta^-$. Now, by an infinite-dimensional version of the formula
that expresses the trace of a matrix as a sum of diagonals, we get that
$$Tr(e^{-t\Delta^+})=\int_M Tr(K^+_t(x,x))dvol(x)=\int_M \sum_ie^{-t\lambda_i}||\phi^+_i(x)||^2 dvol(x)$$
and similarly for $Tr(e^{-t\Delta^-})$. One needs therefore, techniques to evaluate the density
$$\sum_ie^{-t\lambda_i}||\phi^+_i(x)||^2 dvol(x)-\sum_ie^{-t\mu_i}||\phi^-_i(x)||^2 dvol(x).$$
More analysis gives an asymptotic expansion for the plus and minus densities of the form
$$ a^{ \pm }{-d/2}(x) t^{-d/2}+a^{ \pm }{d/2+1}(x) t^{-d/2+1}+\cdots $$
where $d$ is the dimension of $M$.
Up to here was completely general, but then the proof begins to involve special cases, or
at least, broad division into classes of cases. But note that even for the special cases mentioned in the original question,
one would essentially carry out the procedure outlined above for a specific operator $P$.
The breakthrough in this line of thinking
came from Patodi's incredibly complicated computations for the operator $d+d^*$
going from even to odd differential forms,
where one saw that the
$a^{+}_i(x)$
and
$a^{-}_i(x)$
canceled each other out locally, that is, for each point $x$, for all the
terms with negative $i$.
Therefore, one could take a limit
$$\lim_{t\rightarrow 0}[\sum_ie^{-t\lambda_i}||\phi^+_i(x)||^2 dvol(x)-\sum_ie^{-t\mu_i}||\phi^-_i(x)||^2 dvol(x)],$$
that he identified with the Euler form. This important computation set the pattern that recurred in all other versions of
the heat kernel approach to index theorems. One proves the existence of an analogous
limit as $t\rightarrow 0$ and identifies it. The identification
as a precise differential form representative for a characteristic class is referred to sometimes as a {\em local index theorem}, a statement
more refined than the topological formula for the global index. There is even a beautiful version of a local {\em families index theorem}
that relates eventually to deep work in arithmetic intersection theory and Vojta's proof of the Mordell conjecture.
As I understand it, Gilkey's contribution was an invariant theory
argument that tremendously simplified the calculation and allowed a differential form representative for
the $\hat{A}$ genus to emerge naturally
in the case of the Dirac operator. And then, I believe there is a $K$-theory argument that deduces the index theorem for a general elliptic operator
from the one for the twisted Dirac operator.
Experts can correct me if I'm wrong, but
from a purely mathematical point of view, essentially all the work on the heat kernel proof was done at this point.
Subsequent {\em interpretations} of the proof, supersymmetry, path integrals, loop spaces, etc., were tremendously
influential to many areas of mathematics and physics, but the mathematical core of the index theorem itself appears to have remained largely unchanged for almost forty years.
