Why is the Thomson Problem so hard to crack? I read about the Thomson Problem in a Wiki Article (https://en.wikipedia.org/wiki/Thomson_problem)

The objective of the Thomson problem is to determine the minimum electrostatic potential energy configuration of N electrons constrained to the surface of a unit sphere that repel each other with a force given by Coulomb's law. The physicist J. J. Thomson posed the problem in 1904[1] after proposing an atomic model, later called the plum pudding model, based on his knowledge of the existence of negatively charged electrons within neutrally-charged atoms.


Minimum energy configurations have been rigorously identified in only a handful of cases.

I was very surprised when I read this because the solutions for the lower cases are very easy. What makes it difficult for scientists to generalize it to $N$ cases? Why is a computer not able to solve this for general $N$?
Is the Wikipedia Article dubious, especially since it lists solved cases for up to $N=400$ (which are not at all "handful")
Then the article also lists this strange looking conjecture,

According to a conjecture, if $m = n+2$, $p$ is the polyhedron formed by the convex hull of $m$ points, $q$ is the number of quadrilateral faces of $p$, then the solution for $m$ electrons is $f(m) = 0^n +3n -q$

Please, is the Thomson problem really tough?
If yes, Why?!
If no, How has it been solved for arbitrary value of $N$?
 A: This is not an answer, but is too long for a comment, and is in fact in response to a previous comment which stated that all stable 2D configurations will have charges residing on the circumference of the disc.  This is simply not true.  In fact a rather suprising result occurs when n (the number of charges) reaches 12.
The following argument follows a paper that, if I remember correctly, appeared in Nature in the 80’s.  I can’t find the paper, and I know that the arguments below are not exactly those of that article (althought must be very similar - from memory).  In this area of research, the following is a very well known result.  What it also shows is that the 2D problem is a great distance from ‘very easy' (as claimed in a comment from the OP)  
It is observed via numerical calculations that the minimisation of 
$$W\left( n \right)=\sum\limits_{i<j}{\frac{1}{\left| {{r}_{i}}-{{r}_{j}} \right|}}$$
over the surface of a disc, at some critical $n={{n}_{c}}$, has configurations of charges some of which exist on the surface/middle of the disc.  Prior to this, all minimal configurations have charge, equally spaced upon the circumference.  Therefore to determine ${{n}_{c}}$ consider first the energy of the system when charges are equally distributed over the circumference of the disc – we find:
$$W\left( n \right)=\frac{n}{2}\sum\limits_{m=1}^{n-1}{\frac{1}{2\sin \left( \frac{2\pi }{n}m/2 \right)}}=\frac{n}{4}\sum\limits_{m=1}^{n-1}{\frac{1}{\sin \left( \frac{\pi }{n}m \right)}}$$
We want to show $W\left( n \right)>W\left( n-1 \right)+n-1$ which is the energy of a system with one charge at the center of the disc.  First from the well known identity (Mellin Transform).
$$\frac{1}{\sin \left( \tfrac{1}{n}\pi s \right)}=\frac{n}{\pi }\int\limits_{0}^{\infty }{\frac{{{u}^{s-1}}}{{{u}^{n}}+1}du}$$
We have 
$$\frac{n}{4}\sum\limits_{m=1}^{n-1}{\frac{1}{\sin \left( \frac{\pi }{n}m \right)}}=\frac{{{n}^{2}}}{4\pi }\int\limits_{0}^{\infty }{\frac{{{u}^{n-1}}-1}{\left( {{u}^{n}}+1 \right)\left( u-1 \right)}du}$$
Let $u={{e}^{ax}}\to du=a{{e}^{ax}}dx$ for a>0 and real.
$$\frac{n}{4}\sum\limits_{m=1}^{n-1}{\frac{1}{\sin \left( \frac{\pi }{n}m \right)}}=\frac{{{n}^{2}}a}{4\pi }\int\limits_{-\infty }^{\infty }{\frac{{{e}^{axn}}-{{e}^{ax}}}{\left( {{e}^{anx}}+1 \right)\left( {{e}^{ax}}-1 \right)}dx}=\frac{{{n}^{2}}a}{8\pi }\int\limits_{-\infty }^{\infty }{\left( \coth \left( \frac{ax}{2} \right)\tanh \left( \frac{nax}{2} \right)-1 \right)dx}$$
Let $a=2$
$$\frac{W\left( n \right)}{{{n}^{2}}}=\frac{1}{4\pi }\int\limits_{-\infty }^{\infty }{\left( \coth \left( x \right)\tanh \left( nx \right)-1 \right)dx}$$
And so
$$\frac{W\left( n \right)}{{{n}^{2}}}-\frac{W\left( n-1 \right)}{{{\left( n-1 \right)}^{2}}}=\frac{1}{2\pi }\int\limits_{0}^{\infty }{\left( \frac{\tanh \left( nx \right)-\tanh \left( \left( n-1 \right)x \right)}{\tanh \left( x \right)} \right)dx}$$
Note 
$$\frac{\tanh \left( nx \right)-\tanh \left( \left( n-1 \right)x \right)}{\tanh \left( x \right)}=\cosh \left( x \right)\text{sech}\left( nx \right)\text{sech}\left( \left( n-1 \right)x \right)$$
So consider 
$$I=\frac{1}{2\pi }\int\limits_{0}^{\infty }{\cosh \left( x \right)\text{sech}\left( nx \right)\text{sech}\left( \left( n-1 \right)x \right)dx}$$
Now $\text{sech}\left( z \right)\tilde{\ }2{{e}^{-z}}$ so $\text{sech}\left( nx \right)\text{sech}\left( \left( n-1 \right)x \right)\tilde{\ }4{{e}^{-\left( 2n-1 \right)x}}$.  Now consider
$$\text{sech}\left( nx \right)\text{sech}\left( \left( n-1 \right)x \right)-\text{sec}{{\text{h}}^{\text{2}}}\left( anx \right)\tilde{\ }4\left( {{e}^{-\left( 2n-1 \right)x}}-{{e}^{-2anx}} \right)$$ 
for some a.   If $a>1-\frac{1}{2n}$ then
$\frac{1}{2\pi }\int\limits_{0}^{\infty }{\cosh \left( x \right)\text{sech}\left( nx \right)\text{sech}\left( \left( n-1 \right)x \right)dx}$ > $\frac{1}{2\pi }\int\limits_{0}^{\infty }{\cosh \left( x \right)\text{sec}{{\text{h}}^{\text{2}}}\left( anx \right)dx}$ 
so let $a=1$ and we have 
$$\frac{W\left( n \right)}{{{n}^{2}}}-\frac{W\left( n-1 \right)}{{{\left( n-1 \right)}^{2}}}>\frac{1}{2\pi }\int\limits_{0}^{\infty }{\cosh \left( x \right)\text{sec}{{\text{h}}^{\text{2}}}\left( nx \right)dx}$$
But 
$$\frac{1}{2\pi }\int\limits_{0}^{\infty }{\frac{\cosh \left( x \right)}{{{\cosh }^{2}}\left( nx \right)}dx}=\frac{\csc \left( \frac{\pi }{2n} \right)}{4{{n}^{2}}}$$
which we may obtain from properties of the beta function for example.  So
$$\frac{W\left( n \right)}{{{n}^{2}}}-\frac{W\left( n-1 \right)}{{{\left( n-1 \right)}^{2}}}>\frac{\csc \left( \frac{\pi }{2n} \right)}{4{{n}^{2}}}$$
Write this as
$$W\left( n \right)-W\left( n-1 \right)>\frac{\csc \left( \pi /2n \right)}{4}+\frac{\left( 2n-1 \right)W\left( n-1 \right)}{{{\left( n-1 \right)}^{2}}}$$
If $W\left( n \right)-W\left( n-1 \right)>n-1$ we need
$$\frac{\csc \left( \pi /2n \right)}{4}+\frac{\left( 2n-1 \right)W\left( n-1 \right)}{{{\left( n-1 \right)}^{2}}}>n-1$$
First let this be an equality and so we have 
$$\frac{W\left( n-1 \right)}{{{\left( n-1 \right)}^{2}}}=\frac{4\left( n-1 \right)-\csc \left( \pi /2n \right)}{\left( 2n-1 \right)4}$$
Letting n become very large will give us the lowest possible bound.  Noting $\csc \left( \pi /2n \right)\tilde{\ }2n/\pi $therefore yields
$$\frac{W\left( n \right)}{{{n}^{2}}}\ge \frac{1}{2}-\frac{1}{4\pi }$$
The argument from here goes something like this (my notes get blurry from this point)
Numerical computation shows $W\left( n \right)-W\left( n-1 \right)-n+1<0$ for $n<12$, and $>0$ for n=12 and n=13.  Also numerical computation shows$\frac{W\left( 13 \right)}{{{13}^{2}}}-\frac{1}{2}+\frac{1}{4\pi }>0$ which implies $W\left( n \right)-W\left( n-1 \right)>n-1$ for all $n\ge 14$.   Hence we have shown that if the system contains 11 charges or less, then they will settle on the circumference.  Add one more, and it must necessarily be elsewhere. 
Can someone provide a reference for this result perhaps?
Edit: The paper i was thinking of is the following
A.A.BEREZIN, An unexpected result in classical electrostatics, Nature
315 (1985) 104.
This reference came from the following paper by A. Worley arxiv:physics/0609231v1 
A: The configuration space $X$ of $N$ points distributed on the surface of a $2$-sphere is the $N$-fold product $(S^{2})^{N}$ with the "diagonals" removed, i.e., tuples $(x_{1}, \dots, x_{N})$ with $x_{i} = x_{j}$ for some $i \neq j$. The Coulomb energy may be viewed as a function $E:X \to \mathbf{R}$ that grows without bound at the ends (loosely, "the boundary") of the space.
There are at least four difficult aspects to minimizing $E$ on $X$:

*

*("Large domain") The domain of $E$ is $2N$-dimensional. (Consider the problem of numerically minimizing a "high-frequency" function of $2N$ variables over a cube. If you merely bisect each edge of the cube, you have to inspect $2^{2N} > 10^{0.6N}$ cubelets. This is, to put it mildly, computationally infeasible.)


*("Complicated topology") The domain of $E$ is non-compact, and has $\frac{1}{2}N(N - 1)$ ends; there's no simple, useful analytic description of $X$.


*("Complicated critical point structure") The energy function has many unstable critical points and local minima. (For instance, equally space $N$ points around some equator, or equally space $(N - 2)$ points around some equator and put the remaining particles at the poles, etc.)
To develop a sense of how delicate critical points might be, it's worth considering the Figure-8 paper of Alain Chenciner and Richard Montgomery. They proved analytically that there exists a stable periodic figure-8 solution of the $3$-body problem in celestial mechanics. The "energy well" in which the solution sits is about $12$ orders of magnitude shallower than the energy of the configuration itself. One isn't likely to detect such phenomena numerically/experimentally. Again, that's three particles, albeit in a dynamical setting.


*("Lack of recursive structure") A locally-minimizing configuration for $N$ particles gives no general strategy for finding a minimizing configuration for $(N + 1)$ particles.
A: Consider a 2d-analogue of this problem.  It should be "easier", right?  We pack circles into a circle, all of the same diameter, which we try to maximize. This has been studied, and is itself a very hard problem.  Even for 14 circles, it is unknown what the largest possible diameter is.  Not so easy at all.
