$\newcommand{\cov}{\operatorname{cov}}\newcommand{\E}{\operatorname{E}}$I'm going to assume that where you wrote "multinormal" you meant "multinomial" and not multivariate normal.
The probability that you wrote above should be
$$
P(N_1=n_1,N_2=n_2,\ldots,N_r=n_r)=\frac{n!}{n_1!n_2!\cdots n_r!} p_1^{n_1} p_2^{n_2} \cdots p_r^{n_r}.
$$
You have $n$ independent trials and at each trial you choose one of $r$ alternatives, with probabilities $p_1,\ldots,p_r$. The random variable $N_k$ is how many times the $k$th alternative was chosen.
For $k=1,\ldots,r$, let $M_k$ be the number of times the $k$th alternative was chosen on the first trial, so that $M_k$ must be either $0$ or $1$.
A consequence of the fact that the $n$ trials are independent and identically distributed is that $\cov(N_k,N_j) = n\cov(M_k,M_j)$. Assuming $k\ne j$, we can now write
\begin{align}
& \cov(M_k,M_j) = \E(M_k M_j) - \E(M_k)\E(M_j) \\[8pt]
= {} & \Pr(M_k=1=M_j) - \Pr(M_k=1)\Pr(M_j=1) \\[8pt]
= {} & 0 - p_kp_j.
\end{align}
The correlation is the covariance divided by the product of the two standard deviations. The standard deviation is the square root of the variance. The variance is
\begin{align}
& \operatorname{var}(M_k) = \E(M_k^2) - (\E(M_k))^2 \\[8pt]
= {} & \Pr(M_k=1) - (\Pr(M_k=1))^2 = p_k - p_k^2 = p_k(1-p_k).
\end{align}
This may be familiar to you from seeing it stated that the variance of a Bernoulli distribution is $pq$, where $p+q=1$. The random variable $M_k$ is Bernoulli-distributed.