We want to calculate next sum in different point in limit of large $N, N_f$. Initial expression is expansion around $h=0$ by definition. And the answer is known . This is ($N_f/N= \kappa$) $$ \lim_{N \to \infty} \frac{1}{N^2} Z = \kappa^2 \log \frac{1}{1- h^2} \ . $$ We start from $$ Z= \sum_{r=0}^{N N_f} \ h^{2 l} \sum_{\tau \vdash r }s_{\tau}(1^{N_f})s_{\tau}(1^{N_f}) = \prod_{n=1}^{N} \sum_{l_n=0}^{N_f} \ \ e_{l_n} (H) \det_{1\le i,j\le N} e_{l_i-i+j} (H)= $$ where $s_{\sigma}(1^{N_f})$ is Schur function and $\sigma \vdash r$ run over partition(with restriction $λ_1\leq N$ for $λ=(λ_1,..,λ_l)$), or alternatively we are using $e_{l_n} (H)$ - elementary symmetric polynomials, where $H=\left(\underbrace{h,\ldots,h}_{N_f \text{ entries}}\right)$. (See Ira Gessel comments Cauchy identity, with sum restricted over partitions with first part $\leq n$ ). We need to reexpand this around $h=1$.

Lets do it: $$ = \prod_{n=1}^{N} \sum_{l_n=0}^{N_f} \ h^{2 l_n} \ \binom{N_f}{l_n} \det_{1\le i,j\le N} \binom{N_f}{l_i-i+j} =|h^{2 l_n} =1+2 l_n \log h + \cdots | $$ $$ = \det_{1\le i,j\le N} \binom{N_f^2}{N_f-i+j} + \log h \sum_{k=1}^N \det_{1\le i,j\le N} (N_f+ k-j) \binom{N_f^2}{N_f -i+j} + $$ $$ {} + (\log h)^2 \left \{ \sum_{k=1}^N \det_{1\le i,j\le N} \frac{(N_f+ k-j)(N_f^2+(1-N_f)(k-j))}{2 N_f -1} \binom{N_f^2}{N_f-i+j} \right. + $$ $$ \left. +\sum_{k, l=1}^N \det_{1\le i,j\le N} (N_f+ k-j) (N_f+ l-j) \binom{N_f^2}{N_f-i+j} \right \} + $$ $$ {} + (\log h)^3\left \{\sum_{k=1}^N \det_{1\le i,j\le N}\!\!\!\! \frac{ (N_f+k-j)^2 \left(N_f(N_f+1)+ (2- N_f)(k-j)\right)}{3(2 N_f-1)} \!\!\binom{N_f^2}{N_f-i+j}+ \cdots \right \}+ O( (\log h)^4) . $$ Using $$ \sum_{l=0}^{N_f} \ \binom{N_f}{l} \binom{N_f}{l- i+j} = \binom{2 N_f}{N_f-i+ j} $$ we came to $$ \det_{1\le i, j \le N} \ \binom{2 N_f}{N_f-i +j} = s_{N^{N_f}} \left(1^{2N_f}\right)= \det_{1\le i,j \le N} e_{N_f -i+j} \left(\underbrace{1,\ldots,1}_{2N_f \text{ entries}}\right) = $$ $$ =\frac{ G[N+ 2 N_f+1] G[N+1] G[N_f+1]^2 } { G[2 N_f+1] G[N+ N_f+1]^2 }\ = Z_0 \ . $$ (See https://math.stackexchange.com/questions/2681139/determinant-of-a-matrix-with-binomial-coefficient-entries?rq=1). First correction is simple $$ \sum_{k=1}^N \det_{1\le i,j\le N} (N_f+ k-j) \binom{N_f^2}{N_f-i+j} = N_f N \, Z_0 $$ because of $$ \sum_{k=1}^N \det_{1\le i,j\le N} ( k-j) \binom{N_f^2}{N_f-i+j} = \epsilon^{j_1 .. j_N} [(1- j_1)+\cdots +(N- j_N)] $$ $$ \times\binom{N_f^2}{N_f-1+j_1} \cdots \binom{N_f^2}{N_f-N+j_N} = 0 \ . $$ But with next terms we get in trouble (looking at $- 2k\cdot j$ terms) $$ \sum_{k=1}^N \det_{1\le i,j\le N} ( k-j)^2 \binom{N_f^2}{N_f-i+j} = \epsilon^{j_1 \cdots j_N} \left \{ [(1+ j_1^2)+ \cdots +(N^2+ j_N^2)] + \right. $$ $$ \left. {} - 2( j_1+ 2 j_2+ \cdots+ N j_N) \right \} \binom{N_f^2}{N_f-1+j_1} \cdots \binom{N_f^2}{N_f-N+j_N}= 2 \frac{N(N+1)(2 N+1)}{6} Z_0 + \text{???} $$ How to take determinants of this type?

Finally we have first terms of expansion: $$ Z= Z_0(1+ N N_f \log h +...) $$ It is very interesting to obtain all terms.