I consider a 2 state Markov chain: $X = \{1,2\}$, transitions are $M(i,j)$ and the matrix has a unique stationary distribution $\pi$:
$$ \pi(1) = \frac{M(2,1)}{2-M(1,1) - M(2,2)} \\ \pi(2) = \frac{M(1,2)}{2-M(1,1) - M(2,2)}. $$
I try to compute $$ \sum_{i,j\in X} \pi(i) (M^k(i,j) -\pi(j)) \mu_i \mu_j $$ where $\mu_1, \mu_2$ are reals. This is supposed to be $$ \pi(1) \pi(2) (\mu_1 - \mu_2)^2 (1-M(1,2) - M(2,1))^k $$but why?
Calculating the powers of $M$ (see here: How to compute the powers of $2\times2$ Markov matrices) shows that the matrix whose $(i,j)^{\rm th}$ entry is $\pi(i)(M^k(i,j)-\pi(j))$ can be written $$A:=\underbrace{\pi(1)\pi(2)(1-M(1,2)-M(2,1))^k}_{\rm scalar}\pmatrix{1&-1\\-1&1}.$$
The expression you want is \begin{eqnarray*} \mu^T A \mu &=&\pi(1)\pi(2)(1-M(1,2)-M(2,1))^k \pmatrix{\mu_1&\mu_2}\pmatrix{1&-1\\-1&1}\pmatrix{\mu_1\\ \mu_2}\\ &=&\pi(1)\pi(2)(1-M(1,2)-M(2,1))^k (\mu_1-\mu_2)^2. \end{eqnarray*}