EDIT. I think that your new conjecture is true.
$\textbf{Proposition 1}$. Let $G\in GL_n(\mathbb{C})$ and $M$ be a Hermitian $\geq 0$ $n\times n$ matrix with $rank(M)=p$. Let $spectrum(G)=(\lambda_i)$ where $|\lambda_1|\geq\cdots\geq |\lambda_n|$.
Then $\lim_{k\to \infty} \mathrm{tr}(G^kM{G^*}^k)^{1/2k}=|\lambda_i|$ for some $i$ s.t. $1\leq i\leq n-p+1$.
$\textbf{Proof}$. I don't write some details. $(e_i)_i$ denotes the canonical basis of $\mathbb{C}^n$.
An easy lemma
$\textbf{Lemma}$. i) If $M\leq N$, then $\mathrm{tr}(G^kM{G^*}^k)^{1/2k}\leq \mathrm{tr}(G^kN{G^*}^k)^{1/2k}$.
ii) If $M$ is changed with $\alpha M$ ($\alpha>0$), then the studied limit is the same.
We may assume that $M=diag(m_1,\cdots,m_p,0_{n-p})$; according to the lemma, we may assume that the $m_i$ are $1$ ($diag(1/2,1/2,0,0)\leq M=diag(1/2,2,0,0)\leq diag(2,2,0,0)$
and we multiply by $2$). In particular, $M^2=M$.
$\mathrm{tr}(G^kM{G^*}^k)=\mathrm{tr}({G^*}^kG^kM^2)=tr(M{G^*}^kG^kM)$ and
$E_k(u)=u^*M{G^*}^kG^kMu=(Mu)^*({G^*}^kG^k)(Mu)$. Note that $Mu$ goes through $span(e_1,\cdots,e_p)$.
Let $(v_i)_i$ be an orthonormal basis of eigenvectors of ${G^*}^kG^k$ associated to the singular values ${\sigma_1}^2(G^k)\geq \cdots\geq{\sigma_n}^2(G^k)$.
$Z=span(e_1,\cdots,e_p)\cap span(v_1,\cdots,v_{n-p+1})$ is a vector-space of dimension $\geq 1$.
If $Mu\in Z\setminus\{0\}$, then $Mu=\sum_{i\leq n-p+1}a_iv_i$ and
$${\sigma_{n-p+1}}^2(G^k)\leq \dfrac{E_k(u)}{||Mu||^2}=\dfrac{\sum_{i\leq n-p+1}\sigma_i^2(G^k)|a_i|^2}{\sum_{i\leq n-p+1}|a_i|^2}\leq {\sigma_1}^2(G^k).$$
Let $f(u)=\min\{i;a_i\not= 0\}$ and $j=min_u f(u)$.
Then $j$ is constant for $k$ great enough (to check) and $\dfrac{E_k(u)}{||Mu||^2}$ behaves essentially like ${\sigma_j}^2(G^k)$ (up to a factor) -to check-.
Note that, if the $(b_i)$'s are $>0$, then $\lim_{k\to \infty} (b_1^{2k}+\cdots+b_n^{2k})^{1/2k}=\sup(b_1,\cdots,b_n)$.
Finally, roughly speaking, we seek
$\lim_{k\to \infty} ({\sigma_j}^2(G^k))^{1/2k}=\lim_{k\to \infty} (\sigma_j(G^k))^{1/k}$.
According to a result by Yamamoto, the above limit is $|\lambda_j|$. $\square$
$\textbf{Proposition 2}.$ For a generic matrix $G$,
$\lim_{k\rightarrow +\infty}(tr(G^kMG^{*k}))^{1/2k}=\rho(G)$.
$\textbf{Proof}$. Consider a generic $G$ (for example, randomly choose it).
Up to an orthonormal change of basis, we may assume that $M=diag(m_1,\cdots,m_p,0_{n-p})$ where $m_i>0$ and $p\geq 1$.
Let $u=(u_i)$ be a unitary vector s.t. $G^*u=\overline{\lambda} u$ where $|\lambda|=\rho(G)$. Generically, $u$ is unique (up to a factor) and $G$ is diagonalizable.
Then $E(u)=u^*G^kMG^{*k}u=u^*\lambda^k(\overline{\lambda^k} Mu)=\rho(G)^{2k}u^*Mu$ and
$E(u)=\rho(G)^{2k}\sum_i m_i|u_i|^2$.
Generically, at least one of the $u_i$ is non-zero (in fact all the $u_i$) and
$E(u)\sim a\rho(G)^{2k}$ when $k\rightarrow +\infty$ where $a=\sum_i m_i|u_i|^2$.
If $v$ is another unitary vector, then
$E(v)=O(||G^{*k}||||G^k||)=O(\rho(G)^{2k})$.
Then, for $k$ great enough, $tr(G^kMG^{*k})\geq \dfrac{a}{2}\rho(G)^{2k}$ and $tr(G^kMG^{*k})=O(\rho(G)^{2k})$ and we are done. $\square$
Remark. if you choose (for example) $G$ and $M$ as diagonal matrices, then, of course, we may find another limit than the one above.
Best Answer
So startting from this post that deals with the homogeneous case, we take variation of constants to handle with $D$: (the verification of the formula will also be part of the calculation).
$$ \sigma=e^{tA}C(t)(e^{tA})^T\\ \Rightarrow \sigma'=Ae^{tA}C(t)(e^{tA})^T+e^{tA}C'(t)(e^{tA})^T+e^{tA}C(t)(e^{tA}A)^T $$
here we have used basic properties of the matrix exponential and the matrix product rule, see e.g. formulas (33), (39) and (476) in the matrix cookbook. Using the ODE we see that
$$ \sigma'=\underbrace{Ae^{tA}C(t)(e^{tA})^T}_{=A\sigma}+e^{tA}C'(t)(e^{tA})^T+\underbrace{e^{tA}C(t)(e^{tA}A)^T}_{=\sigma A^T}\stackrel{!}{=}A\sigma+\sigma A^T +D\\ \Leftrightarrow e^{tA}C'(t)(e^{tA})^T=D\\\ \Leftrightarrow C'(t)= e^{-tA}D(e^{-tA})^T\\ \Leftrightarrow C=\int_0^t e^{-sA}D(e^{-sA})^Tds+C_0 $$
This gives $$ \sigma=e^{tA}\left(\int_0^t e^{-sA}D(e^{-sA})^Tds\right)(e^{tA})^T+e^{tA}\sigma_0(e^{tA})^T $$
as the full solution of the matrix ODE.