This is only a hint on how to approach this problem, where I try to point out a number of questions that have to be clarified before a complete answer can be given. We assume that $\nu$ is a probability measure.
First, as an operator, $L^t$ needs a domain. We assume initially that the domain is the Hilbert space $H = L^2(M, \nu)$, but this will be taken up for revision at the end.
Lets also split the operator into two parts. Let $k_t(x,y) = e^{-||x-y||^2/4t}$. The first term in the operator can be written as
$$f \mapsto \theta^t f$$
where $\theta^t = \int_M k_t(x,y) d\nu(y)$. Thus it is simply a scalar multiplication of the identity operator. The second term is
$$f \mapsto I^t(f) := \int_M f(y) k_t(\cdot, y) d\nu(y),$$
which is an integral operator. Since $\nu$ is a probability measure and $|k_t(x,y)| \leq 1$ it follows that $I^t(f) \in H$ and $I^t$ is a bounded operator. As I understand the question, empirical versions of these two operators are formed by replacing the measure $\nu$ by an empirical measure $\varepsilon_n$ obtained by sampling $n$ points independently from $\nu$.
Focusing on the integral operator $I^t$ and its empirical approximation $I^t_n$ we would like to know if $I^t_n$ is a good approximation of $I^t$ and in which ways. If we fix $f \in H$ and $x \in M$ then
$$I_n^t(f)(x) = \frac{1}{n} \sum_{i=1}^n f(x_i) k(x, x_i),$$
and since the random variables $f(x_i) k(x,x_i)$ for fixed $f$ and $x$ are i.i.d. with mean $I^t(f)(x)$ it follows from the usual law of large numbers that
$$I^t_n(f)(x) \rightarrow I^t(f)(x) \quad \text{a.s.}$$
for $n \to \infty$. Since $f(x_i)^2 k(x_i, x)^2 \leq f(x_i)^2$ and $f \in H$ is square integrable, the random variables have finite variance and the standard CLT can be used to assess the accuracy.
However, we might want to say something about the approximation uniformly over $f$ and / or $x$. If we fix $f \in H$ does it hold that
$$||I^t_n(f) - I^t(f)||_2 \rightarrow 0 \quad \text{a.s.}$$
for $n \to \infty$? Moreover, since $I_n^t$ and $I^t$ are bounded operators on $H$ we could also ask if
$$||I^t_n - I^t|| \rightarrow 0 \quad \text{a.s.}$$
for $n \to \infty$ with $||\cdot||$ the operator norm on the Banach space of bounded operators?
In both cases above we ask for a law of large numbers to hold in a certain vector space $-$ the Hilbert space $H$ or the Banach space of bounded operators. A solid reference is Probability in Banach Spaces by Ledoux and Talagrand. Central limit results stating that $\sqrt{n}(I^t_n - I^t)$ or $\sqrt{n}(I^t_n(f) - I^t(f))$ are approximately Gaussian processes can also be obtained in some cases.
To obtain good results it might be necessary to restrict attention to a different domain than $H$. One possibility is a Sobolev space for nice choices of $M$ where the functions have some degree of smoothness.
