Thanks! I actually forgot to specify a very important detail in the question, i.e. I am talking about probabilities, so everything is non negative. Anyway, the question can be generalised on how to interpolate between two different continuous linear operators, in general

Yeah @Christian Remling I figured that out don't worry, and thanks again! I just needed the input to look at it as an ODE!

Fun fact is that I am actually dealing with an Euler method, I don't know why I didn't think to simply check what ODE I was solving.... Thanks!

I just want to extend to $[0,T] \times \mathbb{R}^d$, or whatever the initial domain in space was (consider a torus if you like it more). I didn't write it in the question, but the interval $[0,T]$ is for the time variable and $\mathbb{R}^d$ for the space variables. (Edited the question,hope it's clearer).

A question: the formula you provide holds in full generality, right? No hypothesis on triangulability of the matrix or something like that?

I think this might settle it. $\sigma_1$ is the trace, while for the other $\sigma_i$s, I can bound the eigenvalues in terms of the $L^\infty$ norm of the matrix (Gershgoring theorem if I am not mistaken). If I make the entries of $A$ scale like, say, $\varepsilon ^{-1/2}$, the higher order terms in the sum scale as $\varepsilon$ or as $o (\varepsilon)$.