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[%]