We looked into this with Greta and Damir in [this paper][1].  Unfortunately, other than the Harris-Willenbring identity we didn't find much.  For our bounds we even had to invent a new (not too difficult) LR-identity which follows from the *skew Cauchy identity* (it's Lemma 5.1 in the paper):
$$
\sum_{\lambda\vdash n} \left( c^{\lambda}_{\mu\nu}\right)^2 \, = \, \sum_{\alpha, \beta,
\gamma, \delta}  c^{\mu}_{\alpha\gamma}  c^{\mu}_{\alpha \delta} 
c^{\nu}_{\beta\gamma} c^{\nu}_{\beta \delta}
$$


  [1]: https://arxiv.org/abs/1804.04693