Estimate on sum of squares of multinomial coefficients

Douglas already commented that the asymptotics for fixed $p$ and $l\to \infty$ shoudl follow from standard methods. One gets $$a_{\ell}^p\approx (p+1)^{2\ell+\frac{p+1}{2}}(4\pi \ell)^{-\frac{p}{2}}.$$ See theorem 4 in "Counting Abelian squares", by Richmond and Shallit. Notice that these numbers appear also in combinatorics when considering abelian squares, or more generally abelian powers, on a fixed alphabet.

For the asymptotics that you're interested in, at least in the unweighted case, one can say $$a _{\ell} ^p=\sum _{j=0} ^{\ell} \binom{p}{j}\sum _{a _1+ \cdots +a _j = \ell \atop a _i \geq 1} \binom{\ell}{a _1,a _2,\dots,a _j}^2$$ which makes it clear that $a _{\ell}^p$ is a polynomial in $p$ of fixed degree $\ell$. The coefficient of $\binom{p}{\ell}$ is $(\ell!)^2$, and the coefficient of $\binom{p}{\ell -1}$ is $\frac{\ell-1}{4}(\ell!)^2$, so you have $$a _{\ell}^p =\ell!p^{\ell}-\ell!\frac{\ell(\ell-1)}{4}p^{\ell-1}+O(p^{\ell-2}).$$