Integral of the entrywise square of the exponential of a matrix Note: I posted my question on math.stackexchange but got no answer. That is why I am asking it here.
Let $A$ be a $n\times n$ square matrix such that the real part of all eigenvalues are negative. For each $i,j$, let $\exp(At)_{ij}$ be the element $(i,j)$ of the matrix. It is well known that: 
$$ \int_0^\infty \exp(At)_{ij}dt = -(A^{-1})_{ij}$$
Is it possible to simplify a similar expression where each element is squared: 
$$ \int_0^\infty (\exp(At)_{ij})^2 dt = ??$$
I am wondering if it is possible to simplify the above expression. If it helps, I can assume that $A$ is diagonalizable.  Note that unless for one-dimensional matrices, $(\exp(At)_{ij})^2\ne\exp(2At)_{ij}$.
 A: You can do something. Here's a computation for diagonalizable $A$. Let $A=BDB^{-1}$ and let the elements of $D$ be 
$-\lambda_1,\ldots,-\lambda_d$. Then
\begin{align*}
\int_0^\infty (e^{At})_{ij}^2\,dt&=
\int_0^\infty (Be^{Dt}B^{-1})_{ij}^2\,dt\\
&=\int_0^\infty \sum_{k,k'} B_{ik}e^{-\lambda_k t}B^{-1}_{kj}
B_{ik'}e^{-\lambda_k' t}B^{-1}_{k'j}\,dt\\
&=\int_0^\infty \sum_{k,k'} B_{ik}B^{-1}_{kj}
B_{ik'}B^{-1}_{k'j}e^{-(\lambda_k+\lambda_k') t} \,dt\\
&=\sum_{k,k'} \frac{1}{\lambda_k+\lambda_k'}B_{ik}B^{-1}_{kj}
B_{ik'}B^{-1}_{k'j}.
\end{align*}
A: Given a Hurwitz matrix $\mathrm A \in \mathbb R^{n \times n}$, let
$$\Phi (t) := \exp(\mathrm A t)$$
be the state transition matrix, and let its $(i,j)$-th entry be denoted by
$$\varphi_{ij} (t) := \mathrm e_i^{\top} \Phi (t) \, \mathrm e_j$$
Hence,
$$\begin{array}{rl } \displaystyle\int_0^{\infty} \left( \varphi_{ij} (t) \right)^2 \, \mathrm d t &= \displaystyle\int_0^{\infty} \left( \mathrm e_i^{\top} \Phi (t) \, \mathrm e_j \right)^2 \, \mathrm d t\\\\ &= \displaystyle\int_0^{\infty} \mathrm e_i^{\top} \Phi (t) \, \mathrm e_j \mathrm e_j^{\top} \Phi^{\top} (t) \, \mathrm e_i \, \mathrm d t\\\\ &= \mathrm e_i^{\top} \underbrace{\left( \displaystyle\int_0^{\infty} \Phi (t) \, \mathrm e_j \mathrm e_j^{\top} \Phi^{\top} (t) \, \mathrm d t \right)}_{=: \mathrm W_c}  \mathrm e_i = \mathrm e_i^{\top} \mathrm W_c \mathrm e_i\end{array}$$
where $\mathrm W_c$ is the controllability Gramian of the pair $(\mathrm A, \mathrm e_j)$ and is the solution to the following controllability Lyapunov equation
$$\boxed{\mathrm A \mathrm W_c + \mathrm W_c \mathrm A^{\top} + \mathrm e_j \mathrm e_j^{\top} = \mathrm O_n}$$
Thus, the $n$ columns of the integral of the entrywise product
$$\int_0^{\infty} \left( \Phi (t) \circ \Phi (t) \right) \mathrm d t$$
are the diagonals of $\mathrm W_c^{(1)}, \mathrm W_c^{(2)}, \dots, \mathrm W_c^{(n)}$, where $\mathrm W_c^{(1)}, \mathrm W_c^{(2)}, \dots, \mathrm W_c^{(n)}$ are the solutions to the following $n$ controllability Lyapunov equations
$$\begin{array}{cl} \mathrm A \mathrm W_c^{(1)} + \mathrm W_c^{(1)} \mathrm A^{\top} + \mathrm e_1 \mathrm e_1^{\top} &= \mathrm O_n\\ \mathrm A \mathrm W_c^{(2)} + \mathrm W_c^{(2)} \mathrm A^{\top} + \mathrm e_2 \mathrm e_2^{\top} &= \mathrm O_n\\ \vdots  & \\ \mathrm A \mathrm W_c^{(n)} + \mathrm W_c^{(n)} \mathrm A^{\top} + \mathrm e_n \mathrm e_n^{\top} &= \mathrm O_n\end{array}$$
A: Inspired strongly by Anthony's answer, here is a formula that works for arbitrary $A$. Let $M$ be the $n^2 \times n^2$ square matrix given by $$M= A \otimes I_n + I_n \otimes A_n$$
i.e. in terms of indices
$$ M_{ij,kl}=A_{i,k} \delta_{j,l} + \delta_{i,k}A_{j,l}$$
Then because $A \otimes I_n$ and $I_n \otimes A$ commute, $$e^{Mt} = (e^{At} \otimes I_n) (I_n \otimes e^{At}) = e^{At} \otimes e^{At}$$ i.e. in indices $$(e^{M t})_{ij,kl} = (e^{A t})_{i,j}  (e^{At})_{k,l}$$  so in particular $$(e^{At})_{i,j}^2 = (e^{Mt})_{ii,jj}$$ and $$\int_t (e^{At})_{i,j}^2 = \int_t (e^{Mt})_{ii,jj} = - (M^{-1})_{ii,jj}$$
(Here I am using commas to separate the two indices of a matrix entry)
