Here's a simple solution using `Nearest`. I've written it assuming the sites are in a list `ps`, rather than embedded inside a function `p`; I think this is easier to generate and manipulate. You may want to modify the code as appropriate.

    num = 5000;
    ps = RandomInteger[{1, 100}, {num, 2}];

    nf = Nearest[Table[ps[[i]] -> i, {i, Length[ps]}]]

    upToNthNearestSites[k_, 0] := {k} (* the "zeroth" nearest neighbour, i.e. itself *)
    upToNthNearestSites[k_, n_] := Module[{pk, near, d},
      pk = ps[[k]];
      near = upToNthNearestSites[k, n - 1]; (* get the nearest neighbours up to order n-1 *)
      near = nf[pk, Length[near] + 1]; (* get one more; this is one of the nth nearest *)
      d = N@EuclideanDistance[pk, ps[[Last@near]]]; (* distance to the nth nearest *)
      nf[pk, {Infinity, d}] (* the solution is all the sites up to that distance *)
      ]
    nthNearestSites[k_, n_] := Module[{pk, near0, near, d},
      pk = ps[[k]];
      near0 = upToNthNearestSites[k, n - 1];
      near = nf[pk, Length[near0] + 1];
      d = N@EuclideanDistance[pk, ps[[Last@near]]];
      near = nf[pk, {Infinity, d}];
      Complement[near, near0] (* same as above except remove neighbours closer than n *)
      ]

The nearest neighbours to the $k$th site are given by `nthNearestSites[k, 1]`, the second nearest by `nthNearestSites[k, 2]`, and so on. On my machine, even with a *million* points, the initial construction of `nf` takes a little over a second, and after that `nthNearestSites[1, 2]` takes negligible time.

P.S.

1. You *could* just define `nthNearestSites[k_, n_] := Complement[upToNthNearestSites[k, n], upToNthNearestSites[k, n - 1]]`, but that would end up evaluating the $(n-1)$th neighbours twice (as well as the $(n-2)$th, the $(n-3)$th, and so on). In the implementation above, it makes exactly $2n$ calls to the `NearestFunction`.

2. I'm not too happy about having to put the `N` around `EuclideanDistance`. Unfortunately, `NearestFunction` doesn't accept something like $\sqrt2$ as the search radius.