The complete answer is that, $x_i$ maximizes the given function, if and only if $\sum_i x_i=0$. Here is a proof base on achille hui's great comment:
$$
\sum\limits_{i<j}\|x_i - x_j\|^2 = \frac12\sum\limits_{i,j}\|x_i - x_j\|^2 = \frac12\big({\small 4\sum\limits_j \|x_j\|^2 + 4\sum\limits_i \|x_i\|^2 - 2\langle \sum\limits_i x_i , \sum\limits_j x_j}\rangle\big) $$
$$=16-\|\sum_i x_i\|^2.
$$
In particular, any regular tetrahedron whose vertices lie on the sphere satisfies this condition, but there are other examples, such as a square on the equator.
Edit 1:
It is worth noting that the proof above does not use in any essential way the dimension of the sphere (here $2$), or the number of points (here $4$).
Indeed, if we want to place $N$ points $x_1,\dots,x_N$ on the $n$-dimensional sphere $\mathbb{S}^n \subseteq \mathbb{R}^{n+1}$, maximizing $\sum\limits_{i<j}\|x_i - x_j\|^2$, then the same calculation gives:
$$
\sum\limits_{i<j}\|x_i - x_j\|^2 = \frac12\sum\limits_{i,j}\|x_i - x_j\|^2 =\frac12\sum\limits_{i,j}\|x_i\|^2+\|x_j\|^2-2 \langle x_i , x_j\rangle
$$
Since $$\sum\limits_{i,j}\|x_i\|^2=\sum_i \sum_j \|x_i\|^2=\sum_i N \|x_i\|^2=N \sum_i \|x_i\|^2=N^2
$$
and
$$
\sum\limits_{i,j} \langle x_i , x_j\rangle=\langle \sum\limits_i x_i , \sum\limits_j x_j\rangle,
$$
we get
$$
\sum\limits_{i<j}\|x_i - x_j\|^2=\frac12 \big( 2N^2 - 2\langle \sum\limits_i x_i , \sum\limits_j x_j\rangle\big)
=N^2-\|\sum_i x_i\|^2.
$$
Edit 2: (Generalization)
Here we again focus only on four points:
$f : (0,4] \to (0,\infty)$ be any monotone increasing concave function.
We use $\sum_{i<j} |x_i - x_j|^2
=16 - \left|\sum_i x_i\right|^2.$
Since $f$ is concave and monotone increasing,
$$\sum_{i<j}f\left(|x_i - x_j|^2\right) \le 6f\left(\frac16\sum_{i<j}|x_i-x_j|^2\right) = 6f\left(\frac83-\frac16|\sum_ix_i|^2\right)\le 6f\left(\frac83\right).$$
In particular, $E_f=\sum_{i<j}f\left(|x_i - x_j|^2\right)$ is maximized whenever $|x_i-x_j|$ are all equal, i.e. $x_i$ is a tetrahedron. (This automatically implies $\sum_i x_i=0$).
When $f$ is strictly concave, the tetrahedron is the unique maximizer.
In particular this works for the map $x \to \sqrt{x}$.
Let
- $f : (0,4] \to (0,\infty)$ be any monotonic decreasing convex function.
- $p_1,\ldots, p_4$ be any $4$ distinct points on $S^2$.
- $q_1,\ldots, q_4$ be $4$ points on $S^2$ forming a regular tetrahedron.
Notice
$$\begin{align}\sum_{i<j} |p_i - p_j|^2
&= \frac12\sum_{i,j}|p_i - p_j|^2
= \frac12\left(4\sum_j |p_j|^2 + 4\sum_i|p_i|^2 - 2\sum_{ij} p_i\cdot p_j\right)\\
&= 16 - \left|\sum_i p_i\right|^2
\end{align}$$
Since $f$ is convex.
$$\sum_{i<j}f\left(|p_i - p_j|^2\right) \ge 6f\left(\frac16\sum_{i<j}|p_i-p_j|^2\right) = 6f\left[\frac83-\frac16\left|\sum_ip_i\right|^2\right]$$
Since $f$ is monotonic decreasing,
$$\sum_{i<j}f\left(|p_i - p_j|^2\right) \ge 6f\left(\frac83\right)$$
Notice $\displaystyle\;|q_i - q_j|^2 = \frac83\;$ for all $i \ne j$, we have
$$\sum_{i<j}f\left(|p_i - p_j|^2\right) \ge \sum_{i<j}f\left(|q_i - q_j|^2\right)$$
Now the map $\displaystyle\;x \mapsto \frac{1}{\sqrt{x}}\;$ is monotonic decreasing and convex on $(0,4]$. This leads to
$$\sum_{i<j} \frac{1}{|p_i - p_j|} \ge \sum_{i<j} \frac{1}{|q_i - q_j|}$$
As a result, the electrostatic potential energy is minimized when $p_1,\ldots,p_4$ are vertices of a regular tetrahedron.
Best Answer
Nope.
Given any two points $p,q \in S^2$ separated by an angle $\theta \in (0,\pi]$, let $c = \cos\theta$ and $s = \sin\theta$. We have
$$\frac{1}{|p-q|} = \frac{1}{\sqrt{2(1-c)}}$$ Fix $q$ and perturb $p$ by a small amount to $p' = p + \delta p\in S^2$ so that $c \to c + \delta c$, we have
$$\delta\left[\frac{1}{|p - q|}\right] \stackrel{def}{=} \frac{1}{|p' -q|} - \frac{1}{|p-q|} = \frac1{\sqrt{2}}\left(\frac{\delta c}{2\sqrt{1-c}^3} + \frac{3(\delta c)^2}{8\sqrt{1-c}^5} \right) + O((\delta c)^3)\tag{*1} $$
Let $a \in S^2$ so that $a \perp p$ and $q = c p + s a$, Let $b = p \times a$.
For any $0 < \epsilon \ll \theta$, consider perturbations of $p$ to a point $p'$ at angle $\epsilon$ from it.
ie. pick a $\phi$ from $[0,2\pi)$ and set
$$p' = P(\epsilon,\phi) \stackrel{def}{=} \cos\epsilon p + \sin\epsilon (a \cos\phi + b \sin\phi)$$
The corresponding $\delta c$ equals to $c(\cos\epsilon - 1) + s\sin\epsilon\cos\phi$. Substitute this into $(*1)$ and only keeping terms up to $\epsilon^2$, we obtain
$$\delta\left[\frac{1}{|p-q|}\right]_{p'=P(\epsilon,\phi)} = \frac1{\sqrt{2}}\left[\frac{2s\epsilon\cos\phi - c\epsilon^2}{4\sqrt{1-c}^3} + \frac{3(s\epsilon\cos\phi)^2}{8\sqrt{1-c}^5} \right] + O(\epsilon^3)$$ Taking angular average over $\phi$, we obtain
$$\begin{align} \frac{1}{2\pi} \int_0^{2\pi} \delta\left[\frac{1}{|p-q|}\right]_{p'=P_(\epsilon,\phi)} d\phi & = \frac{\epsilon^2}{\sqrt{2}}\left[-\frac{c}{4\sqrt{1-c}^3} + \frac{3s^2}{16\sqrt{1-c}^5} \right] + O(\epsilon^3)\\ &= \frac{\epsilon^2}{\sqrt{2}}\frac{3-c}{16\sqrt{1-c}^3} + O(\epsilon^3) \end{align} $$ Since $\displaystyle\;\frac{3-c}{16\sqrt{1-c}^3} > 0\;$, we find for all sufficient small $\epsilon$,
$$\frac{1}{2\pi}\int_0^{2\pi} \frac{d\phi}{|P(\epsilon,\phi) - q|} > \frac{1}{|p-q|}\tag{*2}$$ For any $4$ points $p_1,\ldots,p_4$, replace $p$ by $p_1$ and $q$ by $p_2,p_3,p_4$ and sum over the pairs.
$(*2)$ tell us for all sufficient small $\epsilon$, the electrostatic potential energy
$$E(p_1,p_2,p_3,p_4) \stackrel{def}{=} \sum_{i<j} \frac{1}{|p_i - p_j|}$$ satisfies
$$\frac{1}{2\pi} \int_0^{2\pi}E( P(\epsilon,\phi),p_2, p_3, p_4 ) > E(p_1,p_2,p_3,p_4)$$
For such $\epsilon$, we can always find a point $p'$ at an angle from $p_1$ with $$E(p_1',p_2,p_3,p_4) > E(p_1,p_2,p_3,p_4)$$
This means $p_1,p_2,p_3,p_4$ cannot be a local maximum for $E(\cdots)$.
Since this is true for any $4$ points on $S^2$, $E(\cdots)$ doesn't have any local maximum.