I like to think of the covariance as a symmetric bilinear form. Indeed, for $X,Y,Z\in L^2(\Omega)$ and $\lambda \in \mathbb{R}$
- $\mathrm{Cov}(X,Y)=\mathbb{E}[(X-\mathbb{E}[X])(Y-\mathbb{E}[Y])]=\mathbb{E}[(Y-\mathbb{E}[Y])(X-\mathbb{E}[X])]=\mathrm{Cov}(Y,X)$.
- $\mathrm{Cov}(X+\lambda Y,Z)=\mathbb{E}[(X+\lambda Y-\mathbb{E}[X+\lambda Y])(Z-\mathbb{E}[Z])]=\mathbb{E}[(X-\mathbb{E}[X])(Z-\mathbb{E}[Z])+\lambda(Y-\mathbb{E}[Y])(Z-\mathbb{E}[Z])]=\mathrm{Cov}(X,Z)+\lambda\mathrm{Cov}(Y,Z)$
by linearity of the expected value. Obviously $\mathrm{Var}(X)=\mathrm{Cov}(X,X)$. Using all those facts leads to
$$ \mathrm{Var}\left(\sum_{j=1}^nX_j\right)=\mathrm{Cov}\left(\sum_{i=1}^nX_i,\sum_{j=1}^nX_j \right)=\sum_{i=1}^n\mathrm{Cov}\left(X_i,\sum_{j=1}^nX_j\right)=\sum_{i=1}^n\sum_{j=1}^n\mathrm{Cov}(X_i,X_j).$$