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.