Inverse of geometric matrix series and its properties

Question

The question begins with the rather simple situation of a geometric series of matrices, namely, let $A \in \mathbb{R}^{N\times N}$ be a real matrix ($A$ can be nonsymmetric) of full rank and $C \in \mathbb{R}^{N}$ a real vector. Suppose that

$$ B = \sum_{j=0}^{\infty} A^j C\, C^\top (A^\top)^j $$

is convergent (we may assume that the largest singular value $\sigma_\max(A) < 1$, for example). We also assume that $\text{rank}( A C \vert A^2 C \vert \ldots \vert A^N C ) = N $.

One can re-write $A^j C = \sum_{i=1}^N c_i \lambda_i \mathbf{v}_i$ for $\{\lambda_1, \ldots, \lambda_N \}$ the distinct eigenvalues of $A$, $\{\mathbf{v}_1, \ldots, \mathbf{v}_N \}$ the associated eigenvectors and coefficients $\{c_1, \ldots, c_N\}$ such that $C = \sum_{i=1}^N c_i \mathbf{v}_i$. Then the above infinite sum simplifies to

$$ B = \sum_{i,j = 1}^N \frac{c_i c_j^H}{1 - \lambda_i \lambda_j^H}\, \mathbf{v}_i \mathbf{v}_j^H $$

The objects of interest are now the following products:

$$ \gamma_k := C^\top (A^\top)^k B^{-1} A^k C \ \in \mathbb{R}^+ , \qquad k \geq 0$$

Question: what can we say in theory of $\gamma_k$ as a function of $k$?

Observations that I have gathered, but have not gotten me very far:

1. By power iteration we generally expect that $A^k C / \|A^k C\| \to \mathbf{v}_1$, the largest eigenvector of $A$ (here for simplicity I assume there aren't two conjugate complex largest eigenvalues, $\lambda_1 \in \mathbb{R}$). Therefore $$ \lim_{k \to \infty} \frac{\gamma_k}{\lambda_1^2} = \mathbf{v}_1' B^{-1} \mathbf{v}_1 > 0$$ The convergence can be slow, of course.

2. For $k \to \infty$, from (1.) it is easy to see that $ \gamma_k \to 0$.

3. For $k$ small - where "small" here depends on $A$, its spectral radius, and $C$ to a lesser extent, but e.g. $0 \leq k \leq 5$ - numerically one finds $$ \gamma_k \approxeq 1 $$ I do not have an intuition of why this would be the case.

4. Matrix $B$ can be proven invertible under the given assumptions above, but it is often very badly determined.

I'd assume this has already come up in some literature, but I could not find anything myself that goes beyond the first part on the geometric matrix series.

Any suggestion is much welcome.

Edit 1: Here is an example plot of $\gamma_k$ for $N = 30$, $A$ and $C$ with elements each randomly (i.i.d.) drawn from a uniform distribution $U([-0.5, 0.5])$ but with $A$ re-scaled to $\overset{\sim}{A} = A / \lambda_\max(A) * 0.9$: