# Reduced density matrix of the edge spin-1/2 in AKLT spin chain

I am trying to understand the paper titled, "Entanglement in a Valence-Bond-Solid State" by Fan, Korepin, and Roychowdhury (https://arxiv.org/abs/quant-ph/0406067).

I was able to understand the theorem which states that the reduced density matrix of some contiguous spin chain block is independent of the location and the total length of the spin chain.

Now, the complication is to tackle the spin-1/2's at the two edges. The reduced density matrix is,

$$\rho_{\hat{L}} = \frac{1}{3^{L}}\sum_{\alpha_1,\ldots \alpha_L=1}^{3} \left( I_{\bar{0}}\otimes \sigma_{ \alpha_L }\ldots \sigma_{ \alpha_1 } \right) | \Psi^{-} \rangle \langle \Psi^{-} | \left( I_{\bar{0}}\otimes \sigma_{ \alpha_L }\ldots \sigma_{ \alpha_1 } \right)^{\dagger}$$ where $L$ is the length of the chain. I don't understand the underlying algebra going from above form to,

$$\rho_{\hat{L}} = \frac{1}{4}\left(1 -p(L) \right) I + p(L) | \Psi^{-} \rangle \langle \Psi^{-}|$$

where $p(L) = (-1/3)^L$. I thought maybe a brute force method might work but I wasn't able to execute it. Other than that, I am not sure how to proceed. Can anyone point me in the right direction?

Note: I do know an alternate approach using the VBS state and the coherent spin representation but I am skeptical if this approach was used by Fan et al.

It is straightforward to see that $\hat\rho_L=(I\otimes\mathcal E)^N(|\Psi^-\rangle\langle\Psi^-|)$, where $$\mathcal E(\rho) = \tfrac13 \sum_\alpha \sigma_\alpha \rho \sigma_\alpha^\dagger\ .$$ It is also straightforward to see that $$\mathcal E(\rho) = \tfrac 43 I - \tfrac13\rho$$ since $\mathcal E(\rho)+\tfrac13\rho$ is the sum over all Paulis including the identity, which maps any $\rho$ to the identity.
This implies that for any bipartite state $\hat\rho$ with reduced density matrix $I$, $$(I\otimes \mathcal E)(\hat\rho) = \tfrac 43 I - \tfrac 13 \hat\rho$$ (which has again reduced density matrix $I$), and in particular $(I\otimes \mathcal E)(I) = I$, and thus $$(I\otimes \mathcal E)^N(\hat\rho) = (1-(-\tfrac 13)^N) I + (-\tfrac 13)^N \hat\rho\ ,$$ where the first term can most easily be inferred from the normalization condition ($\mathcal E$ is trace preserving).
Setting $\hat \rho = |\Psi^-\rangle\langle\Psi^-|$, this yields the result.