I have a few suggestions, but I'm not sure about their efficiency. I'm sure that if this is the only thing that you want to compute my suggestion is good enough.
What I do when I have to compute something like this is to use a lot of replacement rules and a dot product. First, define
P[k_,i_,j_]:=delta[i,j] - k[i] k[j]/k^2
For the P[p-k,i,j] term just use a replacement rule (one of many):
P[p - k, i, j]//. {(a_ - b_)[c_] -> a[c] - b[c]}
gives you the desired output.
Now you only need to use replacement rules, for example
k[i] delta[i,j]-> k[j]
is achieved with
rule={p_[a_] delta[a_, b_] -> p[b]}
for k^2
rule2={p_[i_] p_[i_]->p^2}
The dimension of your system
rule3={delta[i_,i_]->Dim}
You have to make sure to add all the replacement, I just added a few, you need also the case
k[i] delta[j,i]-> k[j]
but now you know how. Finally, you can either keep the use of
k[i] p[i]
or define a dot product
k[i] p[i]->dot[p,k]
I go for the second, always. If you follow my advise, define P as
P[k_,i_,j_]:=delta[i,j] dot[k,k] - k[i] k[j]
now, when you face your expression, keep in mind that you have to divide for (k^2 (k-p)^2). Why do I use the dot? Because (p-k)^2 is not p^2+k^2-2pk. Mathematica will take
p[i]p[i] k[j]k[j]
to be the same as
p[i]k[i] p[j] k[j]
which is wrong. Now it is easier to simplify.
Saludos!