To get arbitrarily many formal variables, you can use Array. But with those variables, your function definition won't work because of the Apply statement. So I modified your definition as follows (I reduced the point number for testing purposes):
pts = Apply[{2 \[Pi] #1, ArcCos[2 #2 - 1]} &, RandomReal[1, {10, 2}], 1]; Clear[energy]; Clear[a]; vars = Array[a, {Length[pts], 2}]; energy[p_] := Module[{cart}, cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]], Cos[#[[1]]]} &, p]; Total[Outer[Exp[-Norm[#1 - #2]] &, cart, cart, 1], 2]]; FindMinimum[energy[vars], Transpose[{Flatten@vars, Flatten@pts}]]
{32.2548, {a[1, 1] -> 1.93787, a[1, 2] -> 1.72361, a[2, 1] -> 1.11355, a[2, 2] -> 0.893035, a[3, 1] -> 6.21077, a[3, 2] -> 2.1405, a[4, 1] -> 3.06917, a[4, 2] -> 2.14062, a[5, 1] -> 1.06997, a[5, 2] -> 2.50937, a[6, 1] -> 4.21367, a[6, 2] -> 1.69561, a[7, 1] -> 5.07748, a[7, 2] -> 2.48594, a[8, 1] -> 4.31041, a[8, 2] -> 0.111206, a[9, 1] -> 4.25016, a[9, 2] -> 3.31368, a[10, 1] -> 5.11923, a[10, 2] -> 0.955784}}
The form of the array passed to energy matches the $N\times2$ dimension list that is expected by the line creating the cart variable. In the FindMinimum statement the dummy variables and initial conditions are specified as a single list of pairs by using Flatten on both.
There is the usual wrinkle that the minimization may need to be tweaked for precision, but that's a different issue.
Finally, to get the minimizing point list, you have to do
vars/.Last[%] Edit
Depending on the function to be optimized, it's sometimes faster to avoid the use of derivatives by specifying the initial conditions for FindMinimum in the form of three numbers:
FindMinimum[energy[vars], Transpose[{Flatten@vars, Flatten@pts, Flatten@pts - .1, Flatten@pts + .1}]] Edit 2
I did get a significant speed-up with this for your example, but the performance depends on the random starting points (and on the choice of bracket width) so I can't say anything definitive. That seems like a topic for a different question.
Edit 3
Though I didn't look at the speed issue in detail, forcing FindMinimum to work with numerical derivatives may be the worst option here. That will happen if you define your function energy only for numerical arguments, as in
energy[p : {{_?NumericQ, _?NumericQ} ..}] := followed by either your own or my initial definition above. So although that's a common advice people give in these applications, it is not going to be the fastest approach here.
Edit 4
I just had another idea on how to improve the speed of my solution: the use of Norm might make it harder to estimate the Hessian for this function. And indeed, when I got rid of Norm there was a significant speed gain (note that the initial solution above is already faster than the _?NumericQ approach even when the latter is compiled while mine is not). I think this is worth adding here because Norm seems like a natural thing to use in pair potentials, even if the energy expression becomes more complicated than the one in this question.
So here is the new version, with Norm replaced by Sqrt[(#1 - #2).(#1 - #2)]. Observe that I have now put back the original particle number of 100 because on my laptop this takes less than 8 seconds to evaluate!
pts = Apply[{2 \[Pi] #1, ArcCos[2 #2 - 1]} &, RandomReal[1, {100, 2}], 1]; Clear[energy]; Clear[a]; vars = Array[a, {Length[pts], 2}]; energy[p_] := Module[{cart}, cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]], Cos[#[[1]]]} &, p]; Total[Outer[Exp[-Sqrt[(#1 - #2).(#1 - #2)]] &, cart, cart, 1], 2]]; FindMinimum[energy[vars], Transpose[{Flatten@vars, Flatten@pts}]]