I have the following DDE system:
\begin{equation} \begin{split} \dot{x}_{1} &= -\mu x_{1}(t) + a_{11}f_{11}(\int_{-\infty}^{t} F(t-s) x_{1}(s-\tau_{1})ds) + a_{12} f_{12}(x_{2}(t-\tau_{1}))) \\ \dot{x}_{1} &= -\mu x_{2}(t) + a_{21}f_{21}(x_{1}(t-\tau_{1})) + a_{22} f_{22}(\int_{t}^{-\infty} F(t-s) x_{2} (s-\tau_{1} ds)) \\ \dot{x_{3}} &= \mu x_{1}(t-\tau_{1}) - \mu x_{3}(t)\\ \dot{x_{4}} &= \mu x_{2}(t-\tau_{1}) - \mu x_{4}(t)\\ \end{split} \end{equation}
where
\begin{equation} \begin{split} x_{3} &= \int_{-\infty}^{t} F(t-s) x_{1}(s - \tau_{1})ds\\ x_{4} &= \int_{-\infty}^{t} F(t-s) x_{2}(s - \tau_{1})ds \end{split} \end{equation}
I want to do stability analysis. The equilibrium point is $\bar{\textbf{x}}=(\bar{x_{1}},\bar{x_{2}},\bar{x_{3}},\bar{x_{4}}) = (0,0,0,0)$. We linearize the system around the point. Define $b_{ij} = a_{ij}\frac{df_{ij}}{x_{I}}$ and do the following substituions:
\begin{equation} \begin{split} x_{1}(t) = \tilde{x_{1}}(t) - \bar{x_{1}}\\ x_{2}(t) = \tilde{x_{2}}(t) - \bar{x_{2}}\\ x_{3}(t) = \tilde{x_{3}}(t) - \bar{x_{3}}\\ x_{4}(t) = \tilde{x_{4}}(t) - \bar{x_{4}}\\ \end{split} \end{equation}
Since $\bar{\textbf{x}} = \textbf{0}$ we have $x_{i}(t) = \tilde{x_{I}}(t)$ Then we would have the following system:
\begin{equation} \begin{split} \dot{x}_{1} &= -\mu x_{1}(t) + b_{11}(x_{3}(t)) + b_{12}(x_{2}(t-\tau_{1})) \\ \dot{x}_{1} &= -\mu x_{2}(t) + b_{21}(x_{1}(t-\tau_{1})) + b_{22}(x_{4}(t)) \\ \dot{x_{3}} &= \mu x_{1}(t-\tau_{1}) - \mu x_{3}(t)\\ \dot{x_{4}} &= \mu x_{2}(t-\tau_{1}) - \mu x_{4}(t)\\ \end{split} \end{equation}
Then we want to find characteristic equation and thus we calculate the following determinant:
\begin{equation} \begin{split} \begin{vmatrix} \lambda + \mu & -b_{12} e^{\lambda \tau_{1}} & -b_{11} & 0\\ -b_{21} e^{-\lambda \tau_{1}} & \lambda + \mu & 0 & -b_{22} \\ -\mu e^{-\lambda \tau_{1}} & 0 & \lambda + \mu & 0 \\ 0 & -\lambda e^{-\lambda \tau_{1}} & 0 & \lambda + \mu \\ \end{vmatrix} &= 0 \end{split} \end{equation}
My question is how do we get the determinant matrix? I know that we calculate $det(\lambda I - A)$ where A is the coefficient matrix of the system $(\dot{\textbf{x}} = A\textbf{x})$. How do we get exponentials like in the first row and second column $(-b_{12} e^{\lambda \tau})$? I mean we do $ \lambda I - A $ but we have $\lambda$ values at non-diagonal indices. How does it happen?