First, let's assume $a=1$; for other values we can scale a solution with $\sqrt{a}$.

So we want to minimize $H=\sum_{i,j} (1-\|x_i-x_j\|^2)^2$.

I globally optimized the problem numerically for $n=4$ and $n=5$ and obtained as solutions

for $n=4$, the square with side length $\sqrt{\frac{2}{3}}\approx0.8165$, giving $H_4=\frac{2}{3}$

for $n=5$, the regular pentagon with side length $\sqrt{\frac{5-\sqrt{5}}{6}}\approx0.6787$, giving $H_5=\frac{5}{3}$

So in these cases the "generalizations of the equilateral triangle" are the regular polygons. This might be an indication that the regular polygons could be solutions to the optimization problem for other $n$. In any case the regular polygon yields an upper bound. To find the value of $H$ for the regular polygon one just has to find the optimal side length and then plug in the values. Perhaps surprisingly, $H^{C_n}$ takes a particular simple form:
$$H^{C_n}=\frac{n(n-3)}{6}.$$
So we know $H^n_{min}\leq H^{C_n}$ and this bound is sharp for $n\in \{3,4,5\}$.

For larger $n$, there are different solutions that obtain the same minimum as the regular polygon. For example distributing the points evenly on the vertices of the equilateral triangle gives $$H=\left\lfloor\frac{\binom{n-1}{2}}{3}\right\rfloor,$$ which is the same as $H^{C_n}$ above if $n$ is a multiple of $3$.

Let me give you two pictures of solutions that are not the regular $n$-gons, but obtain the same minimum:

for $n=6$:

for $n=7$:

I think minimizing this function is interesting also in dimensions other than $d=2$, even $d=1$, but mainly $d>2$. Maybe one can exploit the symmetry in the polynomial $H$ in some way. I don't know the references and the connection to explosives that Ryan Budney mentions in his comment.