Okay so after a lot of complicated math (Euler-Maclaurin formula, as well as some umbral calculus), I can answer my own question kinda.
So suppose you want $C,p,q,r$ s.t.
$$
C = \lim_{n\to\infty} \frac{\prod_{k=1}^n k^{p(k)}}{n^{q(n)}e^{r(n)}}
$$
Instead of looking at $a_{p,k}$, it turns out Bernoulli polynomials occur pretty naturally as $q$'s
So we want to estimate $w(n) = \prod_{k=1}^n k^{p(k)}$
Taking logs, we get $\ln w(n) = \sum_{k=1}^n p(k) \ln(k)$
Letting $\phi(x) = p(x) \ln x$, Euler Maclaurin then says
$$
\ln w(n) = \phi(1) + \int_1^n \phi(x) dx + \sum_{k=1}^p \frac{B_k^+}{k!} \left(\phi^{(k-1)}(n) - \phi^{(k-1)}(1)\right) + R_p
$$
Now note we want something of the form $w(n) = c_nn^{q(n)}e^{r(n)}$, and so the $q(n)$ term arises out of the $\ln n$ terms of $\ln w(n)$.
Now by integration by parts, $$\begin{align*}
&\int_1^n p(x) \ln(x) dx \\
&= p^{(-1)}(n) \ln n - \int_1^n \frac{p^{(-1)}(x)}{x} dx
\end{align*}$$
Note we assume the constant term of $p^{(-1)}$ to be zero, so that the right integral does not produce any $\ln$'s
Also, if $f(x) = h(x) \ln(x)$ for some polynomial $h$, we have
$$f'(x) = h'(x) \ln(x) + h(x)/x$$
Noticing that $\phi^{(|p|)}$ is the last derivative with any $\ln n$ terms, we obtain
$$q(n) = \sum_{n=0}^{|p|+1} \frac{B_k^+}{k!} p^{(k-1)}(n)$$
Defining $P = p^{(-1)}$, and using the umbral calculus linear operator defined by $L(z^k) = B_k^-$, we obtain
$$
\begin{align*}
&q(n) \\
&=\sum_{n=0}^{|P|} \frac{B_k(1)}{k!} P^{(k)}(n) \\
&=\sum \frac{L((z+1)^k)}{k!} P^{(k)}(n) \\
&=L \left( \sum \frac{(z+1)^k}{k!} P^{(k)}(n)\right) \\
&=L(P(z+n+1)) \\
&=\sum_{k=0}^{|P|} P_k L((z+n+1)^k) \\
&=\sum_{k=0}^{|P|} P_k B_k(n+1)
\end{align*}
$$
This implies that if $q(n) = B_k(n+1)$, we obtain $P(n) = n^k$, thus $p(n) = kn^{k-1}$