Given $N$ spins, which eigenspaces of total angular momentum $S^2$ are linked by flipping a single spin?

Let's write the $N$-particle space $H^{otimes N}$ as $H^{otimes N-1}otimes H$ so that we have the basis $$ lvert bar{s},m_{bar{s}}; m_Nrangle$$ where $bar{s}$ is the total spin for the $N-1$ particles and $m_N$ the $z$-spin of the $N$-th particle (since there can be more than one copy of any given $bar{s}$ representation, we would also have to have a marker for that duplicity. I'm omitting it here for convenience of notation as I don't think it changes anything about the argument [Edit: confirmed]). Expressing a state of definite total spin $s$ in terms of this is a standard application of Clebsch-Gordan coefficients: $$ lvert s,m_srangle = sum_{bar{s}}sum_{m_{bar{s}}} sum_{m_N}C^{sm_s}_{bar{s}m_{bar{s}}frac{1}{2}m_N}lvert bar{s},m_{bar{s}}; m_Nrangle,tag{1}$$ where the coefficients are only non-zero for $lvert bar{s}-frac{1}{2}rvertleq sleq bar{s}+frac{1}{2}$ and $m_s = m_{bar{s}} + m_N = m_{bar{s}} pmfrac{1}{2}$.

The question is now when $langle s_1, m_{s_1}lvert sigma^x_N lvert s_2,m_{s_2}rangle$ is non-zero. In the expansion (1), the action of $sigma^x_N$ is just to flip $m_N$ but doesn't touch the $bar{s}$ part. Since the $ lvert bar{s},m_{bar{s}}; m_Nrangle$ are an orthonormal basis, that means this overlap can only be non-zero when the two $lvert s_i,m_{s_i}rangle$ have at least one non-zero $bar{s}$ in common. This can only happen when that $bar{s}$ fulfills both $lvert bar{s}-frac{1}{2}rvertleq s_1leq bar{s}+frac{1}{2}$ and $lvert bar{s}-frac{1}{2}rvertleq s_2leq bar{s}+frac{1}{2}$. Write $s_2 = s_1 + x$, then $$ lvert bar{s}-frac{1}{2}rvert leq s_1leq bar{s}+frac{1}{2} - x,$$ which means $x$ is at most 1, since $lvert bar{s}-frac{1}{2}rvert$ and $bar{s} + frac{1}{2}$ differ at most by 1. This is what we wanted to show.

It is indeed the case that if the two $psi_i$ are in respective spin sectors $s_i$ (i.e., that $S^2|psi_irangle = s_i(s_i+1)|psi_irangle$ for $i=1,2$), then $langle psi_1 | sigma_n^x | psi_2 rangle$ can be non-zero only if $|s_1-s_2|=0$ or $1$. The key idea is directly from ACuriousMind's answer, which I encourage you to read first. What follows is just a lengthy expansion of that answer to explicitly handle the multiplicity of identical irreducible representations. I expect there is a much more elegant way to do this, so alternative answers are encourage.

---

For any Hilbert space $mathcal{H}_{N} = (mathcal{H}_{1})^{otimes N} cong (mathbb{C}^2)^{otimes N} cong mathbb{C}^{2^N}$ of $N$ spin-1/2 particles, we can use the structure theorem for finite-dimensional von Neumann algebras to decompose it as begin{align} mathcal{H}_{N} cong bigoplus_{s={0,1/2}}^{N/2}left(mathcal{S}^{(N)}_s otimes mathcal{E}^{(N)}_sright) end{align} where $mathcal{S}^{(N)}_s cong {mathbb{C}}^{2s+1}$ is an irreducible representation of the SU(2) group generated by ${S^x,S^y,S^z}$ and is spanned by the $2s+1$ orthonormal vectors $|s,mrangle$ for $-sle m le s$. The direct sum starts at $s=0$ or $s=1/2$, depending on whether $N$ is even or odd, respectively, and the dimension of $mathcal{E}^{(N)}_s$ — call it $D(N,s)$ — is just the multiplicity associated with how many copies of the spin-$s$ irrep there are. Any state $|psirangle in mathcal{H}_{N}$ can therefore be written as begin{align} |psirangle &= sum_{s={0,1/2}}^{N/2} sum_{m=-s}^{s} sum_{e=1}^{D(N,s)} alpha(s,m,e) |s:m:erangle end{align} with $alpha(s,m,e)in mathbb{C}$ and where, using an arbitrary orthonormal basis ${|s;erangle}$ of $mathcal{E}^{(N)}_{s}$, we have defined the orthonormal basis begin{align} |s:m:erangle := |s,mrangleotimes|s;erangle in mathcal{S}^{(N)}_{s} otimes mathcal{E}^{(N)}_{s}. end{align} For each fixed value of $s$ and $e$, and for $-sle m le s$, this spans a distinct spin-$s$ irrep of SU(2) in $mathcal{H}_N$. Our assumption that the two $|psi_irangle$ live in sectors with a fixed value $s=s_i$ means that, for them, the respective coefficients $alpha_i(s,m,e)$ vanish except when $s=s_i$, i.e., the sum over $s$ becomes trivial.

We can alternatively decompose any $|psirangleinmathcal{H}_N$ with respect to the tensor structure $mathcal{H}_{N}= mathcal{H}_{N-1}otimes mathcal{H}_{1}$ to get begin{align} |psirangle &= sum_{bar{s}={0,1/2}}^{(N-1)/2} sum_{bar{m}=-bar{s}}^{bar{s}} sum_{f=1}^{D(N-1,bar{s})} sum_{hat{m}=-1/2}^{1/2} beta(bar{s},bar{m},f,hat{m}) |bar{s}:bar{m}:frangle_{N-1} otimes |hat{m}rangle_1 end{align} where $beta(bar{s},bar{m},f,hat{m}) in mathbb{C}$. Here we have used the structure theorem on $mathcal{H}_{N-1}$ but not $mathcal{H}_{1}$ since it is trivial for the latter. Parameters associated with $mathcal{H}_{N-1}$ and $mathcal{H}_{1}$ are indicated with bars ($bar{s}$ and $bar{m}$) and hats ($hat{m}$) respectively. Analogously to before, we have defined the orthonormal basis begin{align} |bar{s}:bar{m}:frangle_{N-1} := |bar{s},bar{m}rangleotimes|bar{s};frangle in mathcal{S}^{(N-1)}_{bar{s}} otimes mathcal{E}^{(N-1)}_{bar{s}}, end{align} which, for each fixed value of $bar{s}$ and $f$, and for $-bar{s}le bar{m} le bar{s}$, spans a spin-$bar{s}$ irrep of SU(2) in $mathcal{H}_{N-1}$. This can be combined with the orthonormal basis $|hat{m}rangle in mathcal{H}_1$ for $hat{m} in{-1/2,1/2}$, which (trivially) spans the lone spin-$1/2$ irrep of SU(2) in $mathcal{H}_{1}$, to form a tensor-product basis for the joint Hilbert space of the pair of spins. In particular, we can perform the Clebsch–Gordan decomposition begin{align} |s:m:bar{s}:frangle:= sum_{bar{m}=-bar{s}}^{bar{s}} sum_{hat{m}=-1/2}^{1/2} C^{s,m}_{bar{s},bar{m};1/2,hat{m}}|bar{s}:bar{m}:frangle_{N-1}otimes|hat{m}rangle_1 end{align} with each allowed fixed choice of $s$, $bar{s}$, and $f$ spanning a distinct spin-$s$ irrep of SU(2) in $mathcal{H}_N$ for $-s le m le s$.

Recall that in our first decomposition of $|psirangle$ we chose an arbitrary orthonormal basis ${|s;erangle}$ for $mathcal{E}^{(N)}_{s}$; all that has happened is that, by fixing a value of $bar{s}$ and $f$ in the state above, we have effectively picked out one of the basis vectors. In other words, there is a choice of basis for $mathcal{E}^{(N)}_{s}$ and a function $e(bar{s},f)$ such that $|s:m:e(bar{s},f)rangle=|s:m:bar{s}:frangle$.

We then have begin{align} beta(bar{s},bar{m},f,hat{m}) &= Big[{}_{N-1} langlebar{s}:bar{m}:f|otimes{}_{1}langlehat{m}|Big]|psirangle &= sum_{s={0,1/2}}^{N/2} sum_{m=-s}^{s} sum_{e=1}^{D(N,s)} alpha(s,m,e) Big[{}_{N-1} langle bar{s}:bar{m}:f | otimes {}_{1}langlehat{m}|Big] |s:m:erangle &= sum_{s={0,1/2}}^{N/2} sum_{m=-s}^{s} alpha(s,m,e(bar{s},f)) C^{s,m}_{bar{s},bar{m};1/2,hat{m}} end{align} where to get the first (second) equality we used our first (second) decompositions of $|psirangle$. To get the third equality we expanded in terms of Clebsch–Gordan coefficients and used the mapping $eto e(bar{s},f)$. Now, following ACuriousMind, we just make the key observation that the Clebsch–Gordan coefficient $C^{s,m}_{bar{s},bar{m};1/2,hat{m}}$ vanishes unless $bar{s}= spm 1/2$. This means that for $|psi_irangle$, for which $alpha_i(s,m,e)$ vanishes unless $s=s_i$, we can conclude that $beta_i(bar{s},bar{m},f,hat{m})$ vanishes unless $bar{s} = s_i pm 1/2$.

Finally, using the second decomposition for $|psi_1rangle$ and $|psi_2rangle$, we see that for any local operator on the $N$th qubit, say $sigma_N^x$, we have begin{align} langlepsi_1|sigma_N^x|psi_2rangle &= sum_{bar{s}={0,1/2}}^{(N-1)/2} sum_{bar{m}=-bar{s}}^{bar{s}} sum_{f=1}^{D(N-1,bar{s})} sum_{hat{m}_1=-1/2}^{1/2} sum_{hat{m}_2=-1/2}^{1/2} beta_1(bar{s},bar{m},f,hat{m}_1)^* beta_2(bar{s},bar{m},f,hat{m}_2) {}_1 langle hat{m}_1|sigma_N^x |hat{m}_2rangle_1 end{align} where we have used the orthonormality of $|bar{s}:bar{m}:frangle_{N-1}$. Having observed that $beta_i(bar{s},bar{m},f,hat{m})$ vanishes unless $bar{s}=s_i pm 1/2$, all terms in the above sum vanish except when $s_1 = s_2-1$, $s_2$, or $s_2+1$.