Assume we have a series of measurements $x_i$ of, say, some physical quantity with (unknown) real value $x_r$. In a physics book I found the statement that, if we could exclude systematic errors in our measurement, we have \begin{equation} (1)\quad x_r=\lim_{n\to\infty}\frac 1n\sum_{i=1}^n x_i. \end{equation} There is also some kind of "proof", starting with definitions \begin{equation} \bar{x}_n:=\frac 1n\sum_{i=1}^n x_i,\quad e_i:=x_r-x_i,\quad \varepsilon:=x_r-\bar{x}, \end{equation} that is, \begin{equation} \varepsilon=\frac 1n\sum_{i=1}^n(x_r-x_i)=\frac 1n\sum_{i=1}^n e_i. \end{equation}
Next step is the calculation \begin{equation} \varepsilon^2=\frac{1}{n^2}\left(\sum_i e_i\right)^2=\frac{1}{n^2}\sum_i e_i^2+\frac{1}{n^2}\sum_i\sum_{j\neq i}e_ie_j\approx\frac{1}{n^2}\sum_i e_i^2 \end{equation} for $n\to\infty$, using that $e_i$ and $e_j$ are statistically independent. I can make no sense of this conclusion, and then, how does this yield the above statement (1)?
This looks an awful lot like the Law of Large Numbers to me.
I'll assume that each of the measurements are independent, identically distributed and (as you say that we can "exclude systematic errors"), each has a mean of $x_r$. So $E(x_i)=x_r$ for each $i$, $E(e_i)=x_r-E(x_i)=0$ for each $i$, and $E(\bar{x}_n)=E\left(\frac{1}{n}\sum_i{x_i} \right)=\frac{1}{n}\sum_iE(x_i)=x_r$ for all $n$.
I'll also edit your book's notation a little: I'll use $\varepsilon_n:=x_r-\bar{x}_n$ instead of just $\varepsilon$ as it seems that the book's use of $\varepsilon$ still depends on $n$.
We have that $E(\varepsilon_n)=x_r-E(\bar{x}_n)=0$ for all $n$. Because of this, $$\mathrm{Var}(\varepsilon_n)=E(\varepsilon_n^2)-E(\varepsilon_n)^2=E(\varepsilon_n^2)-0=E(\varepsilon_n^2)$$ ... so your final line of working suggests that the book wants to be working with the variance of $\varepsilon_n$.
Furthermore, if the $x_i$ are independent, then the $e_i$ are also independent, which means that their covariances are zero. So, for all $i\neq j$: $$0=\mathrm{Cov}(e_i,e_j)=E(e_ie_j)-E(e_i)E(e_j)=E(e_ie_j)-0=E(e_ie_j)$$ ... which suggests why the book ignores the $\frac{1}{n^2}\sum_i\sum_{j\neq i}e_ie_j$ term.
Putting all of this together (and finding the epectation of your final line of working) we get:
$$ \begin{align} \mathrm{Var}(\varepsilon_n) &= E(\varepsilon_n^2) \\ &= \ldots \\ &= \frac{1}{n^2}\sum_i E(e_i^2) \\ &= \frac{1}{n^2}\sum_i \mathrm{Var}(e_i) \\ &= \frac{1}{n^2}\times n \times \mathrm{Var}(e_1) \\ &= \frac{1}{n} \mathrm{Var}(e_1) \end{align}$$ where we've used the fact that if the $e_i$ are identically distributed then they all have the same variance.
Intuitively now, you can see that if $n$ gets large then the variance of your error term $\varepsilon_n$ becomes small. More rigorously, we could use Chebyshev's inequality to show that $\varepsilon_n$ converges in probability to zero and hence that $\bar{x}_n$ converges in probability to $x_r$ (your statement (1)).