Possible Close-form solution for 2 dimensional $\frac{d\Sigma}{dt} = \mathbf{A}\Sigma + \Sigma \mathbf{A}^T + \mathbf{B}\mathbf{B}^T$

When $\mathbf{A},\mathbf{B}$ are 2x2 matrix and $\Sigma(t)$ are PSD, can we expect a close form solution for the following ODE, $$ \frac{d\Sigma}{dt} = \mathbf{A}\Sigma + \Sigma \mathbf{A}^T + \mathbf{B}\mathbf{B}^T. $$

The problem is from solving the covariance of a time-invariant SDE $$ x = \mathbf{A}x + \mathbf{B}dw. $$ When $x$ is in two dimensions, I guess there could be a close-form solution for $\Sigma$. But I am not sure how to solve it.

$ \def\L{\left}\def\R{\right} \def\LR#1{\L(#1\R)} \def\BR#1{\Big(#1\Big)} \def\vecc#1{\operatorname{vec}\LR{#1}} \def\Mat#1{\operatorname{Mat}\LR{#1}} \def\Reshape#1{\operatorname{Mat}\!\BR{#1}} \def\diag#1{\operatorname{diag}\LR{#1}} \def\Diag#1{\operatorname{Diag}\LR{#1}} \def\trace#1{\operatorname{Tr}\LR{#1}} \def\qiq{\quad\iff\quad} \def\grad#1#2{\frac{\p #1}{\p #2}} $Vectorize the equation $$\eqalign{ s &= \vecc{\Sigma} \quad\qiq \Sigma = \Mat{s} \\ b &= \vecc{BB^T} \\ \dot s &= \BR{I\otimes A+A\otimes I}s \;+\; b \\ \dot s &= Cs \;+\; b \\ }$$ Substitute the dependent variable and solve the resulting ODE by inspection $$\eqalign{ x &= s + C^{-1}b \\ \dot x &= \dot s \;=\; Cx \qquad\qquad\qquad\qquad\qquad\quad \\ x(t) &= e^{Ct}\,x(0) \\ }$$ Recover the original variable and de-vectorize to matrix form $$\eqalign{ s(t) &= e^{Ct}s(0) + \LR{e^{Ct}-I}C^{-1}b \qquad\qquad \\ \Sigma(t) &= \Reshape{s(t)} \\ }$$

Update

Since $C$ is the Kronecker sum of $A$ with itself, its exponential can be written in a factored form $$C=\LR{A\oplus A} \qiq e^{Ct} = \LR{e^{At}\otimes e^{At}}$$ which allows the solution to be de-vectorized algebraically $$\eqalign{ V &= \Mat{C^{-1}b} \qiq \vecc{V} = C^{-1}b \\ \Sigma(t) &= e^{At}\,\Sigma(0)\,e^{At} + e^{At}V e^{At} - V \\ }$$ Coincidentally, this provides a solution of the matrix integral in the other answer.

Notice that \begin{align*}\dfrac{d}{dt}\left[e^{-t\mathbf{A}}\mathbf{\Sigma}e^{-t\mathbf{A}^T}\right] &= -e^{-t\mathbf{A}}\mathbf{A}\mathbf{\Sigma}e^{-t\mathbf{A}^T} + e^{-t\mathbf{A}}\dfrac{d\mathbf{\Sigma}}{dt}e^{-t\mathbf{A}^T} - e^{-t\mathbf{A}}\mathbf{A}^T\mathbf{\Sigma}e^{-t\mathbf{A}^T} \\ &= e^{-t\mathbf{A}}\left(-\mathbf{A}\mathbf{\Sigma}+\dfrac{d\mathbf{\Sigma}}{dt} -\mathbf{\Sigma}\mathbf{A}^T \right)e^{-t\mathbf{A}^T} \\ &= e^{-t\mathbf{A}}\mathbf{B}\mathbf{B}^Te^{-t\mathbf{A}^T}.\end{align*}

Hence, $$e^{-t\mathbf{A}}\mathbf{\Sigma}(t)e^{-t\mathbf{A}^T} - \mathbf{\Sigma}(0) = \int_{0}^{t}e^{-s\mathbf{A}}\mathbf{B}\mathbf{B}^Te^{-s\mathbf{A}^T}\,ds,$$ i.e., $$\mathbf{\Sigma}(t) = e^{t\mathbf{A}}\mathbf{\Sigma}(0)e^{t\mathbf{A}^T} + \int_{0}^{t}e^{(t-s)\mathbf{A}}\mathbf{B}\mathbf{B}^Te^{(t-s)\mathbf{A}^T}\,ds.$$