$\newcommand{\S}{\mathbb{S}^2}$Let$$M=\{(x_1,x_2,x_3,x_4) \in \mathbb{S}^2 \times \mathbb{S}^2 \times \mathbb{S}^2 \times \mathbb{S}^2 \, |\,\, \text{ all the } x_i \, \text{ are distinct}\} $$
Let $E:M \to \mathbb{R}$ be defined by $$E(x_1,x_2,x_3,x_4)=\sum_{i < j}\frac{1}{\| x_i - x_j \|},$$
where $\| x_i - x_j \|$ denotes the Euclidean distance in $\mathbb{R}^3$.
Question: Let $p:=(x_1,x_2,x_3,x_4)$ be the configuration of a planar square lying on the equator of $\S$. Is $p$ a local minimum of $E$?
(Of course, it is not a strict minima since one can rotate).
Here is an attempt:
Let $\beta_i(t)$ be a path in $\mathbb{S}^2$, $\beta_i(0)=x_i, \dot \beta(0)=w_i \in T_{x_i}\S$.
Consider the path $$\alpha(t)=(\beta_1(t),\beta_2(t),\beta_3(t),\beta_4(t)).$$ Using $$ \begin{align} &\frac{d}{dt}| \beta_i(t) - \beta_j(t) |^{-1}=\frac{d}{dt}(| \beta_i(t) - \beta_j(t) |^2)^{-\frac{1}{2}}\\&=| \beta_i(t) - \beta_j(t) |^{-3}\big(\langle \dot \beta_i(t), \beta_j(t)\rangle+\langle \beta_i(t), \dot \beta_j(t)\rangle\big), \tag{1} \end{align} $$ we get $$ \frac{d}{dt}E(\alpha(t))=\sum_{i<j}| \beta_i(t) - \beta_j(t) |^{-3}\big(\langle \dot \beta_i(t), \beta_j(t)\rangle+\langle \beta_i(t), \dot \beta_j(t)\rangle\big). \tag{2} $$ In particular, denoting the length of the square's edge by $a$, and assuming that $x_1,x_2,x_3,x_4$ are the square's vertices arranged in a cyclic order we have $$ dE_p(0,0,0,w)=\sum_{i=1}^3 | x_i - x_4 |^{-3}\langle x_i,w\rangle=a^{-3} \langle x_1+x_3,w\rangle+(a\sqrt 2)^{-3} \langle x_2,w \rangle=0, $$ where we used the fact that $x_3=-x_1$, and $x_2=-x_4$, so $\langle x_2,w \rangle=-\langle x_4,w \rangle=0$ as $w \in T_{x_4}\S$.
Differentiating equation $(2)$ again, we get $$ \frac{d^2}{dt^2}| \beta_i - \beta_j |^{-1}=| \beta_i - \beta_j |^{-3}\bigg(3| \beta_i - \beta_j |^{-2}\big(\langle \dot \beta_i, \beta_j\rangle+\langle \beta_i, \dot \beta_j\rangle\big)^2+\langle \ddot \beta_i, \beta_j\rangle+\langle \beta_i, \ddot \beta_j\rangle+2 \langle \dot\beta_i, \dot \beta_j\rangle\bigg)\tag{2}. $$ Now, consider first all the $i<j$ such that $j=i+1 \text{mod} 4$, i.e. $i-j$ is an edge of the square, or equivalently $d_{ij}=|x_i-x_j|=a$.
If we choose e.g. $j=4$, then the two neighbors are $i=1,3$, and so combining terms $1-4,3-4$ we get $$ \langle x_1+x_3, \ddot \beta_4(0)\rangle=0, $$ and similarly for the other two edges $1-2,2-3$. Thus $$ \frac{d^2}{dt^2}|_{t=0}E(\alpha(t))=a^{-3}\bigg( 3a^{-2}\sum_{i=1}^4\big(\langle x_i,\dot \beta_{i+1}\rangle+\langle x_{i+1},\dot \beta_{i}\rangle\big)^2 +2\langle \dot \beta_{i}, \dot \beta_{i+1}\rangle\bigg)+A, $$ where $A$ is the part correspondong to $i-j$ equal $1-3$, $2-4$ (the diagonals).
Consider the pair $2-4$:
$\langle x_2,\dot \beta_{4}\rangle=-\langle x_4,\dot \beta_{4}\rangle=0$, so the first summand $3(a\sqrt 2)^{-2}...$ vanishes.
Thus we are left with $$A=(a\sqrt 2)^{-3}\bigg( \langle \ddot \beta_2, \beta_4\rangle+\langle \beta_2, \ddot \beta_4\rangle+2 \langle \dot\beta_2, \dot \beta_4\rangle \bigg). $$ Since $$ \langle \ddot \beta_2, \beta_4\rangle=-\langle \ddot \beta_2, \beta_2\rangle=|\dot \beta_2|^2, $$ we get $$ A=(a\sqrt 2)^{-3}\bigg( |\dot \beta_2+\dot \beta_4|^2+ |\dot \beta_1+\dot \beta_3|^2 \bigg). $$ Thus, up to a factor of $a^{-3}$, we have $$ \frac{d^2}{dt^2}|_{t=0}E(\alpha(t))=3a^{-2}\sum_{i=1}^4\big(\langle x_i,\dot \beta_{i+1}\rangle+\langle x_{i+1},\dot \beta_{i}\rangle\big)^2 +2\langle \dot \beta_{i}, \dot \beta_{i+1}\rangle+ $$ $$ \sqrt 2^{-3}\bigg( |\dot \beta_2+\dot \beta_4|^2+ |\dot \beta_1+\dot \beta_3|^2 \bigg). $$
Since $a=\sqrt 2$, we get $$ \frac{d^2}{dt^2}|_{t=0}E(\alpha(t))=3/2\sum_{i=1}^4\big(\langle x_i,\dot \beta_{i+1}\rangle+\langle x_{i+1},\dot \beta_{i}\rangle\big)^2 +2\langle \dot \beta_{i}, \dot \beta_{i+1}\rangle+ $$ $$ \sqrt 2^{-3}\bigg( |\dot \beta_2+\dot \beta_4|^2+ |\dot \beta_1+\dot \beta_3|^2 \bigg). $$
Is the last quantity $\ge 0$?
(If I am not mistaken in my computations so far...).
As suggested in the comments, the answer is negative.
The variation that moves one opposite pair toward the north pole and the other opposite pair toward the south pole, all at equal speeds, decreases the energy at second order.
Indeed, we take $$ x_1=(-\frac{1}{\sqrt 2},\frac{1}{\sqrt 2},0), x_2=(\frac{1}{\sqrt 2},\frac{1}{\sqrt 2},0), x_3=(\frac{1}{\sqrt 2},-\frac{1}{\sqrt 2},0),x_4=(-\frac{1}{\sqrt 2},-\frac{1}{\sqrt 2},0), $$ then the given variation corresponds to $$ \dot \beta_1=\dot \beta_3=(0,0,-1),\dot \beta_2=\dot \beta_4=(0,0,1). $$ Then the formula $$ \frac{d^2}{dt^2}|_{t=0}E(\alpha(t))=3/2\sum_{i=1}^4\big(\langle x_i,\dot \beta_{i+1}\rangle+\langle x_{i+1},\dot \beta_{i}\rangle\big)^2 +2\langle \dot \beta_{i}, \dot \beta_{i+1}\rangle+ $$ $$ \sqrt 2^{-3}\bigg( |\dot \beta_2+\dot \beta_4|^2+ |\dot \beta_1+\dot \beta_3|^2 \bigg), $$ gives in this situation $$ 2 \sum_{i=1}^4 \langle \dot \beta_{i}, \dot \beta_{i+1}\rangle+\sqrt 2^{-3}(4+4). $$ Since all $\langle \dot \beta_{i}, \dot \beta_{i+1}\rangle=-1$, we obtain $$ \frac{d^2}{dt^2}|_{t=0}E(\alpha(t))=8(\frac{1}{\sqrt 2^{3}}-1)<0 $$ as required.