Coherence as an (incomplete) measure of the number of correlated spins Consider a system of $N$ spin $1/2$ particles in one dimension. Working in the conventional shared $\hat{\sigma}_i^z$ basis for the system's Hilbert space (spin up/spin down basis), a coherence term of order $q$ is defined as the non-zero (off-diagonal) matrix element $|\sigma_1,...,\sigma_N \rangle \langle \sigma_1',..., \sigma_N'|$ in the density operator, such that $\sum_{i=1}^N (S_i'-S_i)=\frac{1}{2}\sum_{i=1}^N (\sigma_i'-\sigma_i) = q$ with $\sigma_i \in \{-1, 1\}$. The weight of such an order $q$ coherence term in the density operator is then naturally called the $q^{\mathrm{th}}$ (multiple quantum) coherence intensity $I_q$ (denoted as MQC intensity for short). 
Now in the supplementary information section of this paper (Appendix C - page 11), it is stated that:

MQC intensities are an incomplete measure of many-spin
correlations since a signal in $I_q$ indicates there are at least $|q|$ spins present in the correlations. On the other hand, a
correlation with $m$ spins can in principle give rise to all $I_q$ with $q = −m, −m + 1, · · · , m$.

I have a problem with these claims regarding the relationship between correlations and coherences. For example, consider the simple case of $m=2$ uncorrelated spins. By the definition of uncorrelated systems, this means that the total density operator of the two spins can be written as:
$$\hat{\rho} = \hat{\rho}_1 \otimes \hat{\rho_2}$$
where $\hat{\rho}_1$ and $\hat{\rho}_2$ are the reduced density operators of the individual spins. As an example, let the individual (pure) states be $|\psi_1 \rangle = |\psi_2 \rangle = \frac{1}{\sqrt{2}}(|\uparrow \rangle + |\downarrow \rangle)$. We have:
$$\hat{\rho}_1 = \hat{\rho}_2 = \frac{1}{2}(|\uparrow \rangle \langle \uparrow| + |\uparrow \rangle \langle \downarrow|+|\downarrow \rangle \langle \uparrow|+|\downarrow \rangle \langle \downarrow|)$$
The total density operator is then:
$$\hat{\rho} = \hat{\rho}_1 \otimes \hat{\rho_2}=\frac{1}{2}(|\uparrow \rangle \langle \uparrow| + |\uparrow \rangle \langle \downarrow|+|\downarrow \rangle \langle \uparrow|+|\downarrow \rangle \langle \downarrow|) \otimes \frac{1}{2}(|\uparrow \rangle \langle \uparrow| + |\uparrow \rangle \langle \downarrow|+|\downarrow \rangle \langle \uparrow|+|\downarrow \rangle \langle \downarrow|)$$
$$\hat{\rho}= \frac{1}{4}\Big(|\uparrow \uparrow \rangle \langle \uparrow \uparrow|+|\uparrow \uparrow \rangle \langle \uparrow \downarrow|+|\uparrow \uparrow \rangle \langle \downarrow \uparrow|+|\uparrow \uparrow \rangle \langle \downarrow \downarrow|+|\uparrow \downarrow \rangle \langle \uparrow \uparrow|+|\uparrow \downarrow \rangle \langle \uparrow \downarrow|+|\uparrow \downarrow \rangle \langle \downarrow \uparrow|+|\uparrow \downarrow \rangle \langle \downarrow \downarrow|+|\downarrow \uparrow \rangle \langle \uparrow \uparrow|+|\downarrow \uparrow \rangle \langle \uparrow \downarrow|+|\downarrow \uparrow \rangle \langle \downarrow \uparrow|+|\downarrow \uparrow \rangle \langle \downarrow \downarrow|+ |\downarrow \downarrow \rangle \langle \uparrow \uparrow|+|\downarrow \downarrow \rangle \langle \uparrow \downarrow|+|\downarrow \downarrow \rangle \langle \downarrow \uparrow|+|\downarrow \downarrow \rangle \langle \downarrow \downarrow|\Big)$$
Clearly, this state contains all coherence terms $q \in \{-2,-1,0,1,2\}$, even though it is an uncorrelated state for the two spins by definition. This is in direct contradiction with the statement quoted from the paper. 
I feel like the paper is referring to a different kind of correlation than what I am thinking of. What am I missing?
 A: So I found the answer myself after digging through other parts of the paper. It looks like the mentioned statement is only true in the context of the paper, not in general. In the main text, it is shown that for the initial density operator, which is the Gibbs distribution (at high temperature) $\hat{\rho}(0) \simeq \frac1{2^N}(\hat{1} - \sum_{i=1}^{N} \epsilon_i \hat{\sigma}_i^z)$, the time-evolved density operator under an arbitrary unitary transformation has the form:
$$\hat{\rho} = \frac{1}{2^N}(\hat1 - \sum_{\{\alpha_k\}}b_{\{\alpha_k\}} \bigotimes_{k=1}^N \hat{\sigma}^{\alpha_k})$$
(I'm using a different notation than the paper to make my point a little bit more clear.) Here $\alpha_k \in \{0,x,y,z\}$ with $\hat{\sigma}^0 \equiv \hat 1$. The coefficients $b_{\{\alpha_k\}}$ depend on the specific unitary evolution in question. Also, the sum does not include all the $\alpha_k$ being zero, which would give a $\bigotimes_{k=1}^N \hat{\sigma}^{\alpha_k} =\hat1 \otimes \hat1 ...\otimes \hat 1 = \hat 1$. You can see from this that each $\bigotimes_{k=1}^N \hat{\sigma}^{\alpha_k}$ term has a nonzero contribution to the corresponding correlation function $\langle \sigma_1^{\alpha_1}...\sigma_2^{\alpha_2}\rangle$. For example, consider $N=4$ spins, the two spin correlation $\langle \sigma_1^{x}\sigma_2^{z}\rangle$ is:
$$\langle \sigma_1^{x}\sigma_2^{z}\rangle = \mathrm{Tr}(\hat \sigma_1^{x} \hat \sigma_2^{z} \hat \rho) = - b_{x,z,0,0} \mathrm{Tr}\Big((\hat {\sigma}^x)^2 \otimes (\hat {\sigma}^z)^2\otimes \hat 1 \otimes \hat 1\Big) = - b_{x,z,0,0}$$
where I used the fact that $\mathrm{Tr}(\hat A \otimes \hat B) \equiv \mathrm{Tr}(\hat A) \times \mathrm{Tr}(\hat B)$. The other terms in the sum are zero because $\mathrm{Tr}(\hat \sigma^{\alpha_k}) \equiv 0$ for all $\alpha_k \neq 0$. 
Thus, the coefficients $b_{\{\alpha_k\}}$ are simply the correlations between the corresponding components of the particles' spins. Now regarding the coherences, you can see that an $m$ spin correlation (so a term such as $\hat \sigma^{\alpha_1} \otimes ...\hat \sigma^{\alpha_m}\otimes \hat 1\otimes ... \hat 1 $ in the density operator), can in principle give rise to any off-diagonal matrix element $\langle \{\sigma_k\}| \hat \sigma^{\alpha_1} \otimes ...\hat \sigma^{\alpha_m}\otimes \hat 1\otimes ... \hat 1 |\{\sigma'_k\} \rangle$ with coherence between $-m$ and $m$, because $\hat \sigma^x$ and $\hat \sigma^y$ operators basically flip spins, and the $\hat \sigma^z$ operator doesn't. For example a two spin correlation in the $N=4$ case can have the term  $\hat \sigma^x \otimes \hat \sigma^x \otimes \hat 1 \otimes \hat 1$ which gives a $q=m = 2$ and $q= - m = -2$ coherence; or a $\hat \sigma^z \otimes \hat \sigma^z \otimes \hat 1 \otimes \hat 1$ term which gives $q=0$ (writing the other possible terms shows that they also give $q=1$ and $q=-1$). 
In conclusion, the statement given in the paper is true for the specific system with the specific state studied in the paper, but not in general.
A: I suspect they're talking about something simpler. It sounds like they're just stating the fact that for any density matrix $\lvert\rho_{jk}\rvert^2\leq\rho_{jj}\rho_{kk}$. 
If your actual system of $M$ spins is represented by $\rho_{M}$, but you write the density matrix for $N>M$ spins, then you'll be writing $\rho_{M}\otimes\frac{1}{2^{N-M}}\mathbb{1}$; you're free to write the density matrix for the rest of the universe, but it is in a completely unknown state. If you go off-diagonal for any of the extra $N-M$ spins, you'll find a zero in the matrix.
In other words, any matrix element that is not zero can involve off-diagonal elements of only some subset of the $M$ spins that are actually in your system. But note that the converse is not true: for example, if one of the $M$ spins is in the eigenstate $\lvert\uparrow\rangle_j$, then any off-diagonal element involving $\lvert\downarrow\rangle_j$ or $\langle\downarrow\rvert_j$.
