<s>[INST] $\lim_{N\rightarrow\infty}\left< \psi{\left(N_k+1\right)} - \psi{\left(N_j+1\right)} \right >$ 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
[1] Bounding $\sum\limits_{k=1}^n{ \binom{n}{k} \, \frac{\left(-1\right)^k}{k} \left[ p_1^k - p_2^k \right] } $
[2] Identity and bounding of ${ \sum\limits_{k=1}^N{ \binom{N}{k}\,\dfrac{p^k \, \left( - 1 \right)^{k} }{k} } } $ when $0<p<1$?
[3] Proofs of $\lim\limits_{n \to \infty} \left(H_n - 2^{-n} \sum\limits_{k=1}^n \binom{n}{k} H_k\right) = \log 2$
[4] Covariance of product of two functions of two binomial distributions
 [/INST] My answer is correct. Problem solved.
</s>