Question
What is \begin{align} \lim_{N\rightarrow\infty}\left< \psi{\left(N_k+1\right)} - \psi{\left(N_j+1\right)} \right >? \end{align}
Here, $N_k$ and $N_j$ are multinomially distributed random variable with expected values $\left<N_k\right> = N\,p_k$ and $\left<N_j\right> = N\,p_j$, respectively; $0<p_k, p_j <1$; and $\psi$ is the digamma function.
Context
This question implicitly solves the questions in [1] and [2]. The solution is based on a related question [3]; and specifically makes use of of Srivatsan's second method therein. At one point below, a binomial identity in [4] is used.
Approach
To solve this problem, we generalize Srivatsan's second method that is derived for a related problem (see [3]).
\begin{align} \left< \psi{\left(N_k+1\right)} \right > & = - \gamma + \sum\limits_{N_k=0}^{N}{\left[\binom{N}{N_k}\, p_k^{N_k}\,\left( 1 - p_k\right)^{N-N_k} \, \sum\limits_{t=1}^{N_k}{\dfrac{1}{t}}\right]} \\ & = - \gamma + \sum\limits_{N_1=0}^{N}{\left[\binom{N}{N_k}\, p_k^{N_k}\,\left( 1 - p_k\right)^{N-N_k} \, \sum\limits_{t=1}^{N_k}{ \int_{0}^{1}{x^{t-1}} }\right]}\,dx \\ & = - \gamma + \int_{0}^{1} \sum\limits_{N_1=0}^{N}{\left[\binom{N}{N_k} \, p_k^{N_k}\,\left( 1 - p_k\right)^{N-N_k} \, { \dfrac{ x^{N_1} -1 }{x - 1} }\right]}\,dx \\ & = - \gamma + \int_{0}^{1} { { \dfrac{\sum\limits_{N_1=0}^{N}\binom{N}{N_k} \, p_k^{N_k}\,\left( 1 - p_k\right)^{N-N_k} \, x^{N_1} - \sum\limits_{N_1=0}^{N}\binom{N}{N_k} \, p_k^{N_k}\,\left( 1 - p_k\right)^{N-N_k} \, 1 }{x - 1} } }\,dx \\ & = - \gamma + \int_{0}^{1} { { \dfrac{ \left(1 - p_k + x\,p_k\right)^N - 1 }{x - 1} } }\,dx \end{align} In the step above, we use an identity found in [4]. Furhter, $\gamma$ is the Euler-Mascheroni constant.
Next, make the substitution $y = 1 - p_k + x\,p_k$ so $\dfrac{dy}{p_k} = dx$ \begin{align} \left< \psi{\left(N_k+1\right)} \right > & = - \gamma + \int_{1-p_k}^{1} { { \dfrac{ y^N - 1 }{ \dfrac{y - 1 + p_k}{p_k} - 1} } }\, \dfrac{dy}{p_k} \\ & = - \gamma + \int_{1-p_k}^{1} { { \dfrac{ y^N - 1 }{ y - 1 + p_k - p_k} } }\, dy \\ & = - \gamma + \int_{1-p_k}^{1} { { \dfrac{ y^N - 1 }{y - 1} } }\, dy \\ & = - \gamma + \int_{1-p_k}^{1} { \sum\limits_{k=0}^{N-1}{y^k} }\, dy \\ & = - \gamma + \left. \sum\limits_{k=1}^{N}{\dfrac{y^k}{k}} \right|_{1-p_k}^{1} \\ & = - \gamma + H_N - \sum\limits_{k=1}^{N}{\dfrac{\left(1 - p_k\right)^k}{k}} \end{align} where $H_{N}$ is the $N^{\textrm{th}}$ harmonic number.
If one notes that
\begin{equation}
- \ln (1 - z) = z + \frac{z^2}{2} + \frac{z^3}{3}\cdots = \sum\limits_{k=1}^{\infty}{\dfrac{ z ^k}{k}}
\end{equation}
is valid for $\left|z\right| \leq 1$, then the summation $\sum\limits_{k=1}^{N}{\dfrac{\left(1 - p_k\right)^k}{k}}$ is understood to be the $n$th partial sum of the Taylor series expansion of $ -\ln{\left(1-y\right)}$ evaluated at $y= 1- p_k$. Therfore, with respect to the question at hand, we write
\begin{align}
\lim_{N\rightarrow \infty }\left< \psi{\left(N_k+1\right)} - \psi{\left(N_j+1\right)} \right >
&
= \left[- \gamma + \lim_{N\rightarrow \infty }\left\{H_N\right\} - \left(- \log{\left(1- \left[1 - p_k\right]\right)}\right) \right] -\left[- \gamma + \lim_{N\rightarrow \infty }\left\{H_N\right\} - \left(- \log{\left(1- \left[1 - p_j\right]\right)}\right) \right]
\\
&
= \log{\left( p_k \right)} - \log{\left( p_j \right)}
\end{align}
Solution
Finally, we find \begin{align} \lim_{N\rightarrow \infty }\left< \psi{\left(N_k+1\right)} - \psi{\left(N_j+1\right)} \right > = \log{\left( \dfrac{p_k}{p_j} \right)}. \end{align}
Citations
[4] Covariance of product of two functions of two binomial distributions
