By viewing $\{x_1, \ldots, x_n\}$ and $\{y_1, \ldots, y_n\}$ as two separate samples, $\sum r_i^2 = (n - 1)S_x^2 + (n - 1)S_y^2$, where $S_x^2 = (n - 1)^{-1}\sum_{i = 1}^n (x_i - \bar{x})^2$ is the well-known sample variance. So your problem reduces to how to re-express the sample variance as the sum of squared paired differences.
In non-parametric statistics, this formulation is called U-statistic. For sample variance, it is easy to verify that the sample variance is a U-statistic of order-2: \begin{align} S_x^2 = \frac{1}{n - 1}\sum_{i = 1}^n (x_i - \bar{x})^2 = \frac{1}{n(n - 1)}\sum_{1 \leq i \neq j \leq n}\frac{(x_i - x_j)^2}{2}. \end{align}\begin{align} S_x^2 = \frac{1}{n - 1}\sum_{i = 1}^n (x_i - \bar{x})^2 = \frac{1}{n(n - 1)}\sum_{1 \leq i \neq j \leq n}\frac{(x_i - x_j)^2}{2}. \tag{1} \end{align}
If you prefer increasing sum indices (as opposed to distinct sum indices), $(1)$ can be clearly rewritten as: \begin{align} S_x^2 = \frac{2}{n(n - 1)}\sum_{1 \leq i < j \leq n}\frac{(x_i - x_j)^2}{2}. \tag{2} \end{align}
Therefore, the quantity of your interest can be written as \begin{align} \sum_{i = 1}^nr_i^2 = \frac{1}{2n}\sum_{1 \leq i \neq j \leq n}[(x_i - x_j)^2 + (y_i - y_j)^2]. \end{align}\begin{align} \sum_{i = 1}^nr_i^2 &= \frac{1}{2n}\sum_{1 \leq i \neq j \leq n}[(x_i - x_j)^2 + (y_i - y_j)^2] \\ &= \frac{1}{n}\sum_{1 \leq i < j \leq n}[(x_i - x_j)^2 + (y_i - y_j)^2]. \end{align}
Proof of (1). Since $x_i - x_j = 0$ when $i = j$, it follows that \begin{align} & \sum_{1 \leq i \neq j \leq n}(x_i - x_j)^2 \\ =& \sum_{i = 1}^n\sum_{j = 1}^n(x_i - x_j)^2 \\ =& \sum_{i = 1}^n\sum_{j = 1}^n[(x_i - \bar{x}) - (x_j - \bar{x})]^2 \\ =& \sum_{i = 1}^n\sum_{j = 1}^n[(x_i - \bar{x})^2 - 2(x_i - \bar{x})(x_j - \bar{x}) + (x_j - \bar{x})^2] \\ =& n\sum_{i = 1}^n(x_i - \bar{x})^2 + 0 + n\sum_{j = 1}^n(x_j - \bar{x})^2 \\ =& 2n(n - 1)S_x^2. \end{align}