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