Flory characteristic ratio of hindered rotation chain

This is an exercise on Rubinstein's textbook of Polymer Physics.

1. Definition of Flory characteristic ratio:
   $$C_\infty:=\lim_{n\to\infty}\frac{\langle \mathbf{R}^2 \rangle}{nb^2}$$
   where $\langle \mathbf{R}^2 \rangle$ is the mean square end-to-end vector, $b$ is Kuhn length of the polymer and $n$ is corresponding chain length.
2. Calculation of $\langle \mathbf{R}^2\rangle$:
   Let $\mathbf{r}_i$ be the bond vector between $i$ and $(i-1)$th particle, and $\mathbf{r}_1=\mathbf{R}_1-\mathbf{R}_0\equiv0$; for an $n$-bead chain,:
   $$\begin{align}\langle \mathbf{R}^2 \rangle&=\left\langle\left(\sum_{i=1}^n \mathbf{r}_i^2\right)\cdot\left(\sum_{j=1}^n \mathbf{r}_j^2\right)\right\rangle\\
   &=\sum_i\sum_j\langle \mathbf{r}_i\cdot\mathbf{r}_j\rangle\end{align}$$
3. Calculation of $\langle\mathbf{r}_i\cdot\mathbf{r}_j\rangle$:
   Now consider a local coordinate frame that $\mathbf{r}_i=(b,0,0)^T$, then $\mathbf{r}_{i+1}$ could be considered as rotating $\mathbf{r}_i$ $\theta$ by $z$-axis for bond angle and $\phi_{i}$ by $x$-axis for torsion. The rotation matrix is thus
   $$\mathbf{T}_i:=\left(
   \begin{array}{ccc}
   \cos (\theta ) & -\sin (\theta ) & 0 \\
   \sin (\theta ) \cos (\phi_i) & \cos (\theta ) \cos (\phi_i) & -\sin (\phi_i) \\
   \sin (\theta ) \sin (\phi_i ) & \cos (\theta ) \sin (\phi_i ) & \cos (\phi_i ) \\
   \end{array}\right)$$
   or $\pi-\theta$, $\pi - \phi$ according to the definitions of bond angle $\theta$ and torsion angle $\phi$. In every reference coordinate frame $R_i$, the bond vector is $(b,0,0)^T$, and in $R_i$, $\mathbf{r}_{i+1}=\mathbf{T}_i\cdot\mathbf{r}_i=(b \cos (\theta ),b \sin (\theta ) \cos (\phi_i ),b \sin (\theta ) \sin (\phi_i ))^T$. Now consider a simple case of $\mathbf{r}_i\cdot\mathbf{r}_{i+2}$, the $\mathbf{r}_{i+2}$ is calculated in the $(i+1)$th reference frame $R_{i+1}$ where $\mathbf{r}_{i+1}$ is $(b,0,0)^T$ and $\mathbf{r}_{i+2}=\mathbf{T}_{i+1}\cdot\mathbf{r}_{i+1}$, here we must 'rotate' one of the vector into the other frame to do the inner product, i.e. rotate $\mathbf{r}_i=(b,0,0)^T$ from $R_i$ to $R_{i+1}$ or vice versa. Now note that $\mathbf{r}_{i+1}$ in the $R_i$ frame is $(b \cos (\theta ),b \sin (\theta ) \cos (\phi_i ),b \sin (\theta ) \sin (\phi_i ))^T$ and $(b,0,0)^T$ in the $R_{i+1}$ frame, therefore if we want to transform some vector in $R_i$ into $R_{i+1}$, we need a transform matrix to transform the $(b \cos (\theta ),b \sin (\theta ) \cos (\phi_i ),b \sin (\theta ) \sin (\phi_i ))^T$ into $(b,0,0)^T$ which is $\mathbf{T}_i^{-1}$, or $\mathbf{T}_i^T$, because in frame $R_i$, we obtain $\mathbf{r}_{i+1}$ by $\mathbf{T}_i\cdot(b,0,0)^T$ and $\mathbf{T}_i$ is a rotation matrix, which is orthonormal. By transform $\mathbf{r}_i$ into $R_{i+1}$, the inner product $\mathbf{r}_i\cdot\mathbf{r}_{i+2}$ is
   $$\langle\mathbf{T}_i^T\mathbf{r}_i,\mathbf{T}_{i+1}\mathbf{r}_{i+1}\rangle=$$ $$\langle\mathbf{r}_i,\mathbf{T}_i\mathbf{T}_{i+1}\mathbf{r}_{i+1}\rangle=(b,0,0)\mathbf{T}_i\mathbf{T}_{i+1}(b,0,0)^T$$
   By induction, $\mathbf{r}_i\cdot\mathbf{r}_j=(b,0,0)\mathbf{T}_i\mathbf{T}_{i+1}\mathbf{T}_{i+2}\cdots\mathbf{T}_{j-1}(b,0,0)^T$, since $\theta$ for every bonds are same and $\phi_{1,2,\dots,n}$ are independently distributed, let $\mathbf{T}$ be the average transform matrix of ${\mathbf{T}_i}$:
   $$\begin{align}\mathbf{T}&:=\frac{\int_{-\pi}^{\pi}\mathbf{T}_i(\phi)\exp{(-U(\phi)/k_BT)}\mathrm{d}\phi}{\int_{-\pi}^{\pi}\exp{(-U(\phi)/k_BT)}\mathrm{d}\phi}\\
   &=\left(
   \begin{array}{ccc}
   \cos (\theta ) & -\sin (\theta ) & 0 \\
   \sin (\theta ) \langle\cos (\phi) \rangle & \cos (\theta )\langle \cos (\phi)\rangle &0 \\
   0& 0& \langle\cos (\phi)\rangle \\
   \end{array}\right)\end{align}$$
   We therefore have $\langle\mathbf{r}_i\cdot\mathbf{r}_j\rangle=(b,0,0)\mathbf{T}^{|j-i|}(b,0,0^T)=b^2(\mathbf{T}^{|j-i|})_{11}$. Here every $\sin(\phi)$ term vanishes in the integral because $U$ is even.
4. Now the Flory characteristic ratio
   $$\begin{align}C_\infty&=\lim_{n\to\infty}\frac{1}{nb^2}\sum_i\sum_j \mathbf{r}_i\cdot\mathbf{r}_j\\
   &=\lim_{n\to\infty}\frac{1}{n}(\sum_i\sum_j T^{|j-i|})_{11}\\
   &=\lim_{n\to\infty}\left(-\frac{1}{n}\frac{-2 \mathbf{T}^{n+1}+\mathbf{T}^2 n+2 \mathbf{T}-n\mathbf{I}}{(\mathbf{I}-\mathbf{T})^2}\right)_{11}\\
   &=\left(\frac{\mathbf{I}+\mathbf{T}}{\mathbf{I}-\mathbf{T}}\right)_{11}\\
   &=\frac{1+\cos(\theta)}{1-\cos(\theta)}\frac{1-\langle\cos(\phi)\rangle}{1+\langle\cos(\phi)\rangle}\end{align}$$
   Here we have $\frac{\bullet}{n}=0$ and $\mathbf{T}^n=0$ at $n\to\infty$, because $\mathbf{T}$ is constant and correlation between 2 beads goes to $0$ as distance goes to infinity.