Problem
I have a multi-type branching process with $S$ types, where an individual of type $i$ generates $i$ offsprings. The probability that an offspring is of type $j$ is $p_j$, regardless of its parent. Assume we start the process with one individual of type $i_0$, I am interested in determining the variance of the process at the $n$-th generation. That is, by defining $Z^{(n)}_j$ the number of individuals of type $j$ at the $n$-th generation, I am interested in finding $\mathbb{E}\left[Z^{(n)}_jZ^{(n)}_i|Z^{(0)}_k=\delta_{ki_0}\right]$. Finding $\mathbb{E}\left[Z^{(n)}_j|Z^{(0)}_k=\delta_{ki_0}\right]$ is relatively straightforward.
Harris [1] gives a formula for it (page 37), but it is unclear when the expectations are conditioned on the initial conditions or not. Athreya [2] gives the general recipe to find it and cites Harris, but as I tried to work it out I get inconsistent answers.
My attempt (summary)
I find that the variance of the total population at the $n$th generation is: \begin{align} \text{Var}\left[\sum_j Z^{(n)}_j|Z_k^{(0)}=\delta_{ik}\right]&=(n-1)i \mu^{n}(\sigma^2+\mu^2-\mu)+i\mu{n-1}-i^2\mu^{2(n-1)} \end{align} where I defined $\mu=\sum_i ip_i$ and $\sigma^2=\sum_i i^2 p_i - \mu^2$. I believe this to be wrong, as I expect the variance to grow as $n$ increases, but clearly this expression is decreasing (with the leading order $\propto \mu^{2n}$). I do no see where my mistake is. Heuristically I expect that $\mathbb{E}\left[Z^{(n)}_jZ^{(n)}_i|Z^{(0)}_k=\delta_{ki_0}\right]$ should be proportional to $\mu^{2(n-1)}$.
My attempt (full derivation)
$\underline{\text{First Moment}}$
Given an individual of type $i$, the distribution of its offsprings is a multinomial distribution with parameters $(i,\{p_1,p_2,\dots,p_k\})$. The probability generating function of the multinomial is:
\begin{equation} f_i(\mathbf{s})=\left(\sum_j p_j s_j\right)^i. \end{equation}
Following Harris or Athreya, starting with one individual of type $i$, the average number of individuals of type $j$ after one generation is: \begin{equation} \mathbb{E}[Z_{j}^{(1)}|Z_k^{(0)}=\delta_{ik}]=\frac{\partial f_i(\mathbf{1})}{\partial s_j}:=M_{ij} \end{equation} where $\mathbf{1}=\{1,1,\dots,1\}$.
Using a central result in the theory of branching processes, we can represent the $n$-th generation as:
\begin{equation} \mathbb{E}[Z_{j}^{(n)}|Z_k^{(0)}=\delta_{ik}]=\frac{\partial f^{(n)}_i(\mathbf{1})}{\partial s_j}, \end{equation} where \begin{equation} \mathbf{f}^{(n)}(\mathbf{s})=\mathbf{f}\left(\mathbf{f}^{(n-1)}(\mathbf{s})\right)=\overbrace{\mathbf{f}\left(\mathbf{f}\left(\mathbf{f}\left(\dots\right)\right)\right)}^{n\text{ times}} \end{equation}.
Using the chain rule, we can write an expression for $\mathbb{E}[Z_{j}^{(n)}]$:
\begin{align} \mathbb{E}[Z_{j}^{(n)}|Z_k^{(0)}&=\delta_{ik}]=\left.\frac{\partial f^{(n)}_i(\mathbf{s})}{\partial s_j}\right|_{\mathbf{s}=\mathbf{1}}\\ &=\sum_k \left.\frac{\partial f_i}{\partial s_k}\frac{\partial f^{(n-1)}_k}{\partial s_j}\right|_{\mathbf{s}=\mathbf{1}}\\ &=\sum_{k_1,k_2}\left. \frac{\partial f_i}{\partial s_{k_1}}\frac{\partial f_{k_1}}{\partial s_{k_2}}\frac{\partial f^{(n-2)}_{k_2}}{\partial s_j}\right|_{\mathbf{s}=\mathbf{1}}\\ &=\sum_{k_1,k_2,\dots k_{n-1}}\left. \frac{\partial f_i}{\partial s_{k_1}}\frac{\partial f_{k_1}}{\partial s_{k_2}}\dots\frac{\partial f^{}_{k_{n-1}}}{\partial s_j}\right|_{\mathbf{s}=\mathbf{1}}\\ &=\left[M^n\right]_{ij} \end{align} We can compute explicitly $\left[M^n\right]_{ij}$ in our case, since $M_{ij}=ip_j$, we have \begin{equation} \left[M^n\right]_{ij}=ip_{k_1}k_1p_{k_2}\dots k_{n-1}p_j=ip_j\mu^{n-1} \end{equation} where we defined $\mu=\sum_i i p_i$.
$\underline{\text{Second Moment}}$
Following Athreya [2] (p185 Eq.(3)), we can express the second moments as:
\begin{equation} \mathbb{E}\left[Z^{(1)}_jZ^{(1)}_l-\delta_{jl}Z^{(1)}_j|Z^{(0)}_k=\delta_{ik}\right]=\left.\frac{\partial^2 f_i(\mathbf{s})}{\partial s_j\partial s_l}\right|_{\mathbf{s}=\mathbf{1}}:=V^{(l)}_{ij} \end{equation}
So I apply the second derivative on the previous result:
\begin{align} \left.\frac{\partial^2 f^{(n)}_i(\mathbf{s})}{\partial s_j\partial s_l}\right|_{\mathbf{s}=\mathbf{1}}&=\left.\frac{\partial}{\partial s_l}\left(\sum_{k_1,k_2,\dots k_{n-1}} \frac{\partial f_i}{\partial s_{k_1}}\frac{\partial f_{k_1}}{\partial s_{k_2}}\dots\frac{\partial f^{}_{k_{n-1}}}{\partial s_j}\right)\right|_{\mathbf{s}=\mathbf{1}}\\ &=\sum_{k_1,k_2,\dots k_{n-1}} \left[\left.\frac{\partial^2 f_i}{\partial s_{k_1}s_{l}}\frac{\partial f_{k_1}}{\partial s_{k_2}}\dots\frac{\partial f^{}_{k_{n-1}}}{\partial s_j}\right|_{\mathbf{s}=\mathbf{1}}+\left.\frac{\partial f_i}{\partial s_{k_1}}\frac{\partial^2 f_{k_1}}{\partial s_{k_2}\partial s_l}\dots\frac{\partial f^{}_{k_{n-1}}}{\partial s_j}\right|_{\mathbf{s}=\mathbf{1}}+ \dots \\ + \left.\frac{\partial f_i}{\partial s_{k_1}}\frac{\partial f_{k_1}}{\partial s_{k_2}}\dots\frac{\partial^2f^{}_{k_{n-1}}}{\partial s_j\partial s_l}\right|_{\mathbf{s}=\mathbf{1}}\right]\\ &=\sum_{k_1,k_2,\dots k_{n-1}} \left(V^{(l)}_{i{k_1}} M_{{k_1}k_2}\dots M_{k_{n-1}j} + M_{i{k_1}}V^{(l)}_{{k_1}{k_2}}M_{k_2k_3}\dots M_{k_{n-1}j}+\dots + M_{i{k_1}}M_{k_1k_2}\dots M_{k_{n-2}k_{n-1}}V^{(l)}_{k_{n-1}j}\right)\\ &=\left[V^{(l)}M^{n-1}\right]_{ij}+\left[M V^{(l)}M^{n-2}\right]_{ij}+\left[M^2V^{(l)}M^{n-3}\right]_{ij}+\dots \left[M^{n-1}V^{(l)}\right]_{ij}\\ &=\sum_{k=0}^{n-1} \left[M^kV^{(l)}M^{n-1-k}\right]_{ij} \end{align}
Therefore the variance of the total population at the $n$-th generation (given that we start with a known type) can be expressed as:
\begin{align} \text{Var}\left[\sum_j Z^{(n)}_j\right]&=\sum_{jl}\text{Cov}\left[ Z^{(n)}_jZ^{(n)}_l\right]\\ &=\mathbb{E}[Z^{(n)}_jZ^{(n)}_l|Z^{(0)}_k=\delta_{ik}]-\mathbb{E}[Z^{(n)}_j|Z^{(0)}_k=\delta_{ik}]\mathbb{E}[Z^{(n)}_l|Z^{(0)}_k=\delta_{ik}]\\ &=\sum_{jl}\sum_{m=0}^{n-1} \left[M^kV^{(l)}M^{n-1-m}\right]_{ij}+\sum_j \left[M^n\right]_{ij}-\sum_{jl}\left[M^n\right]_{ij}\left[M^n\right]_{il} \end{align}
Similarly to $M_{ij}$, we can explicitly compute $V^{(l)}_{ij}$: \begin{align} V^{(l)}_{ij}=i(i-1)p_jp_l \end{align}
Therefore we have: \begin{align} \sum_{k=0}^{n-1} \left[M^kV^{(l)}M^{n-1-k}\right]_{ij}&=\sum_{k=0}^{n-1}\sum_{ab} ip_a\mu^{k-1}a(a-1)p_lp_b bp_j \mu^{n-1-(k-1)}\\ &=\sum_{k=0}i p_j\mu^{n}(\sigma^2+\mu^2-\mu)\\ &=(n-1)i p_j\mu^{n}(\sigma^2+\mu^2-\mu) \end{align} where we defined $\sigma^2=\sum_i i^2 p_i-\mu^2$.
Unless I made an error above, here is the inconsistency: I would have expected the variance to grow as a factor of $\propto\mu^{2n}$ as $n$ grows, since the mean number of individuals grows as $\propto\mu^{n}$. Writing it explicitly:
\begin{align} \text{Var}\left[\sum_j Z^{(n)}_j\right]&=(n-1)i \mu^{n}(\sigma^2+\mu^2-\mu)+i\mu{n-1}-i^2\mu^{2(n-1)} \end{align}
[1] Harris, T. E. (1963). The theory of branching processes (Vol. 6). Berlin: Springer.
[2] Athreya, K. B., Ney, P. E., & Ney, P. E. (2004). Branching processes. Courier Corporation.