Higher order Lyapunov equation

Question

Let $A$ be a (finite) Hurwitz matrix.

In this related question of mine, (see also https://en.wikipedia.org/wiki/Lyapunov_equation) it is shown that $$ \int_0^\infty \sum_{j,k} (e^{At})_{ij} Q_{jk} (e^{At})_{\ell,k}dt = W_{i,\ell}, $$ where the matrix $W$ is the solution of some matrix equation $AW+WA^T+Q = 0$. Numerically, the solution of such an equation can be computed in $O(d^3)$, where $d$ is the dimension of my system.

I am wondering if this result can be easily extended to multi-linear maps? In particular, is it possible to compute (numerically) the quantity $$ \int_0^\infty \sum_{j_1,j_2,j_3}Q_{j_1,j_2,j_3}(e^{At})_{i_1,j_1} (e^{At})_{i_2,j_2}(e^{At})_{i_3,j_3}dt $$ by an efficient procedure?

Also, if we go to higher order, like $$ \int_0^\infty \sum_{j_1,j_2,\dots,j_k}Q_{j_1,\dots,j_k}(e^{At})_{i_1,j_1} \dots (e^{At})_{i_k,j_k}dt, $$ can we compute the solution of such an equation with a (relatively efficient) numerical procedure? If yes, how does its complexity grow with $k$?

[EDIT] : after my answer below, there still remains a question whether a more efficient algorithm exists (like the Lyapunov equation for $k=2$).

Answer 1

Inspired by Will's answer of Integral of the entrywise square of the exponential of a matrix I realized that the above equation can also be computed by using the Kronecker product.

Let $$M= A \oplus A \oplus A = A \otimes I_n \otimes I_n + I_n \otimes A \otimes I_n + I_n \otimes I_n \otimes A.$$ In terms of indices, this means: $$M_{i_1i_2i_3,j_1j_2j_3} = A_{i_1j_1}\delta_{i_2j_2}\delta_{i_3j_3} + \delta_{i_1j_1}A_{i_2j_2}\delta_{i_3j_3}+ \delta_{i_1j_1}\delta_{i_2j_2}A_{i_3j_3}$$ We have: $$\exp(Mt) = \exp(At)\otimes \exp(At) \otimes \exp(At).$$ In terms of indices this means that: $$\exp(Mt)_{i_1i_2i_3,j_1j_2j_3} = \exp(At)_{i_1,j_1}\exp(At)_{i_2,j_2}\exp(At)_{i_3,j_3}$$ In particular, $$ \int_0^\infty \exp(Mt)dt = -M^{-1},$$ which gives the solution to the above integral in $O(n^9)$.

This can be easily generalized to higher order, with a complexity in $O(n^{3k})$. For $k=2$, it is less efficient than solving the Lyapunov equation (which can be done in $O(n^3)$).