Let $A\in \mathbb{S}^{N\times N}$ be a symmetric, real and stable matrix, i.e., $\rho(A)<1$, where $\rho(A)$ stands for the spectral radius of $A$. Then, $$\sum\limits_{i=0}^{\infty} A^{2i}=\left( I-A^2 \right)^{-1}.\tag{$\star$}\label{star}$$

Now, let $A$ be asymmetric. It is not clear whether the power series $$F(A):=\sum_{i=0}^{\infty} A^i (A^i)^{\top}$$ admits a *compact closed-form* expression as in equation \eqref{star} (when it converges). If the matrix $A$ is normal — i.e., it commutes with its transpose — and stable, then we have that $F(A)=\left(I-AA^{\top}\right)^{-1}$. In general, that is not the case.

Framework.Let $A$ be a random matrix generated as follows. Let $\widetilde{A}\in \mathcal{G}(N,p)$ be the adjacency matrix associated with the realization of a random graph over $N$ nodes with probability $p$ of arrow drawing, i.e., the entries $\widetilde{A}_{ij}$ of $\widetilde{A}$ are i.i.d. $\mathsf{Bernoulli}(p)$ for all $i\neq j$ — this is often referred to as a binomial random graph. Set $\widetilde{A}_{ii}=0$ for all $i$. Let $A$ be obtained from $\widetilde{A}$ by:i) [off-diagonal]normalizing the entries of $\widetilde{A}$ as $A_{ij}:=\alpha_1 \widetilde{A}_{ij}/d_{\max}$, for all $i\neq j$, where $d_{\max}$ is the maximumin-flowdegree of the random graph and $0<\alpha_1<1$ is some constant;ii) [diagonal]setting the diagonal terms as $A_{ii}:=\alpha-\sum_{k\neq i} \widetilde{A}_{ik}$, where $0<\alpha_1\leq\alpha<1$. In other words, the rows of $A$ sum to $\alpha$ and its support is given by the realization of the binomial random graph on $N$ nodes.

**Question.** Is it true that $$\mathbb{P}\left(\lVert\sum_{i=0}^{\infty}A^i \left(A^i\right)^{\top} - \left(I-AA^{\top}\right)^{-1}\rVert_\text{max}>\varepsilon\right)\overset{N\rightarrow \infty}\longrightarrow 0?$$

Is there any known result in this direction? Any reference would be appreciated. The specifics of the random nature of $A$ are not important other than that of $A$ being stable and randomly generated (so that it is more natural to consider the limit $N\rightarrow\infty$).

Some numerical experiments are pointing in the direction of a *concentration* of $\sum_{i=0}^{\infty} A^i \left(A^i\right)^{\top}$ about $\left(I-AA^{\top}\right)^{-1}$ as the dimension $N$ grows large (I may attach them later, if necessary) under the referred framework.

**Context.** $F(A)$ is the covariance matrix associated with the time-series reflecting the state-evolution of certain linear stochastic dynamical systems. Having a closed form expression for the covariance in terms of $A$ is helpful for model identification purposes.