We know that covariance between samples $(X_1,\ldots,X_n)$ and $(Y_1,\ldots,Y_n)$ is
$$\frac{1}{n-1} \sum_{i=1}^n (X_i - \overline{X})(Y_i-\overline{Y}). \qquad (*)$$
But I didn't find the derivation of this formula and don't understand why it doesn't look like the following $\frac{1}{(n-1)^2} \sum_{i=1}^n \sum_{j=1}^n (X_i - \overline{X})(Y_j-\overline{Y})$ (i.e. why the sum has only $n$ terms instead of $n^2$ terms).
It looks like the expression (*) is correct only if samples $(X_1,\ldots,X_n)$ and $(Y_1,\ldots,Y_n)$ satisfy the following condition: $P(X=X_i,Y=Y_j) = 0, \, ~~\forall i \ne j~~$ (here $X$ and $Y$ are random variables which have corresponding population distributions). Can you confirm my guess?
You want to compute covariance of two r.v. $X$ and $Y$ taking into account pairs of observations $(x_1,y_1) \dots (x_n,y_n)$ for both variables. Cross terms (i.e. $(x_i,y_j)$ with $i \neq j$) do not correspond to observations and they don't appear when computing the sample covariance, but it does not necessarily mean that $P(X=x_i,Y=y_j)=0$.