I am trying to compute the following limit
$$\lim_{\epsilon \to 0} \frac{1}{\epsilon} \mbox{tr} \left( S(e^{-\epsilon C}-e^{-\epsilon B})^T(Q_\epsilon)^{-1}(e^{-\epsilon C}-e^{-\epsilon B}) \right)$$
where
$Q_\epsilon,B,C$ are real $n \times n$ matrices.
The spectra of $B$ and $C$ have strictly positive real parts.
$C= S B^T S^{-1}$ for some symmetric positive definite matrix defined as $$S = \int_0^\infty e^{-t B} \sigma \sigma^T e^{-t B^T}\:dt$$
Let $$Q_\epsilon = \int_0^\epsilon e^{-t C} \sigma \sigma^T e^{-t C^T}\:dt$$ where $\sigma$ is an $n \times n$ singular matrix. However, $Q_\epsilon$ has the property of being invertible for $\epsilon \neq 0$. This is equivalent to Kalman's rank condition -- i.e., the block matrix $[\sigma, B \sigma,\ldots B^{n-1}\sigma]$ has rank $n$ -- which we assume as well.
The following relations holds: $BS+SB^T = \sigma \sigma^T$
My conjecture from numerical experiments is that the limit equals $+\infty$, but I haven't been able to prove it.
My attempt so far:
Let $D = \sigma \sigma^T$.
We can show that necessarily, $BS$, $BD$, $CS$, $CD$ are not symmetric.
Since the limit $\lim_{\epsilon \to 0} \frac{1}{\epsilon}(e^{-\epsilon C}-e^{-\epsilon B})^T$ is straightforward, the harder bit is to compute $\lim_{\epsilon \to 0}(Q_\epsilon)^{-1}(e^{-\epsilon C}-e^{-\epsilon B})$
Also, \begin{align*} Q_\epsilon &= \int_0^\epsilon \sum_{n=0}^\infty \left(\frac{(-t)^n}{n!}C^n\right)D \sum_{n=0}^\infty\left( \frac{(-t)^n}{n!}(C^T)^n\right)\:dt \\ &= \int_0^\epsilon \sum_{n=0}^\infty \sum_{k=0}^n \frac{(-t)^{n-k}}{(n-k)!}C^{n-k}D\frac{(-t)^{k}}{k!}(C^T)^{k}\:dt \\ &=\sum_{n=0}^\infty \sum_{k=0}^n \frac{(-1)^n}{k!(n-k)!}\int_0^\epsilon t^nC^{n-k}D(C^T)^{k}\:dt \\ &=\sum_{n=0}^\infty \sum_{k=0}^n \frac{(-1)^n\epsilon^{n+1}}{k!(n-k)!(n+1)}C^{n-k}D(C^T)^{k} \\ \end{align*}
Similarly, \begin{align*} e^{-\epsilon C}-e^{-\epsilon B} &= \sum_{n=1}^\infty \frac{\epsilon^{n}}{n!}((-C)^n-(-B)^n)\\ &= \sum_{n=1}^\infty \frac{(-\epsilon)^{n}}{n!}(C^n-B^n) \end{align*}
So we can reformulate the problem as a quotient of matrix power series:
\begin{align*} \lim_{\epsilon \to 0}\left(\sum_{n=0}^\infty \sum_{k=0}^n \frac{(-1)^n\epsilon^{n}}{k!(n-k)!(n+1)}C^{n-k}D(C^T)^{k}\right)^{-1}\left(\sum_{n=1}^\infty \frac{(-\epsilon)^{n-1}}{n!}(C^n-B^n) \right) \end{align*}
This seems quite difficult since I can't just let $\epsilon \to 0$ as $D$ is singular by assumption.
Define $$\eqalign{ F_\epsilon &= \left(e^{-\epsilon C}-e^{-\epsilon B}\right) \;&\implies\; &\dot F_\epsilon=\frac{dF_\epsilon}{d\epsilon} = \left(Be^{-\epsilon B} - Ce^{-\epsilon C}\right) \\ P_\epsilon &= F_\epsilon SF_\epsilon^T \;&\implies\;&\dot P_\epsilon = \dot F_\epsilon SF_\epsilon^T + F_\epsilon S\dot F_\epsilon^T \\ Q_\epsilon &= \int_0^\epsilon q(t)\,dt \;&\implies\;&\dot Q_\epsilon = q(\epsilon) = e^{-\epsilon C}\sigma\sigma^Te^{-\epsilon C^T} \\ \phi(\epsilon) &= {\rm Tr}(P_\epsilon\,Q_\epsilon^{-1}) &=\quad&P_\epsilon:Q_\epsilon^{-1} \\ }$$ where a colon is used in the final expression as a convenient product notation for the trace.
First calculate the derivative $$\eqalign{ \dot\phi &= \frac{d\phi}{d\epsilon} \\ &= \dot P_\epsilon:Q_\epsilon^{-1} + P_\epsilon:\dot Q_\epsilon^{-1} \\ &= \dot P_\epsilon:Q_\epsilon^{-1} - P_\epsilon:(Q_\epsilon^{-1}\dot Q_\epsilon\,Q_\epsilon^{-1})\\ }$$ L'Hopital's Rule tells us that $$\eqalign{ \lim_{\epsilon\to 0}\left(\frac{\phi}{\epsilon}\right) &= \lim_{\epsilon\to 0}\,\dot\phi \\ }$$ But the problem, as you noted in your question, is that $Q^{-1}_\epsilon$ is singular at $\epsilon=0$.
The other terms in the derivative expression are nicely behaved and are straightforward to calculate. However, the effort is wasted because, when the limit is taken, that singular matrix pops up in every term.