Distribution of $\mathbf{x}^\top \mathbb{A} \mathbf{x} + \mathbf{b}^\top \mathbf{x}$, when $\mathbf{x}\sim\mathcal{N}$? Let $\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, \Sigma)$ be a random vector following a multivariate normal density distribution. I am interested in the density function of the transformed variable
$$z = \mathbf{x}^\top \mathbb{A} \mathbf{x} + \mathbf{b}^\top \mathbf{x}$$
where $\mathbb{A}$ is a constant matrix and $\mathbf{b}$ a constant vector. Is there a name for the distribution associated with $z$?
I presume this related to the Chi-squared distribution, which is a special case of this problem. But I'm not sure my problem can be reduced to the Chi-squared distribution?
Any further ideas?
Update: Upon further thought and thanks to @whuber comment, I see the distribution of $z$ might not have a closed form. A less ambitious goal is then to compute the moments of $z$, starting with the mean, variance, ....
The mean is actually easy:
$$\langle z\rangle = \mathrm{Tr}(\mathbb{A}\Sigma) + \boldsymbol{\mu}^{\top}\mathbb{A}\boldsymbol{\mu} + \mathbf{b}^{\top}\mathbf{u}$$

$\DeclareMathOperator{\Var}{Var}$
$\DeclareMathOperator{\Cov}{Cov}$
$\DeclareMathOperator{\tr}{tr}$
The variance of $z$ can be computed as follows (suppose $A$ is symmetric, which is standard when the quadratic form is discussed. For the non-symmetric case, rewrite $z$ as $z = x'Bx + b'x$, where $B = (A + A')/2$ and replace $A$ below with the symmetric matrix $B$).
In this thread, it is shown that when $x \sim N(\mu, \Sigma)$ and $A$ is symmetric,
\begin{align}
\Var(x'Ax) = 2\tr((A\Sigma)^2) + 4\mu'A\Sigma A\mu. 
\end{align}
Therefore, to evaluate
\begin{align}
\Var(z) = \Var(x'Ax + b'x) = \Var(x'Ax) + \Var(b'x) + 2\Cov(x'Ax, b'x), \tag{1}
\end{align}
it suffices to evaluate
\begin{align}
\Cov(x'Ax, b'x) = E[x'Axb'x] - E[x'Ax]E[b'x]. 
\end{align}
Since $E[x'Ax] = \mu'A\mu + \tr(A\Sigma), E[b'x] = b'\mu$, the only difficult part left is $E[x'Axb'x]$.  To this end, write $x = \mu + y$, where $y \sim N(0, \Sigma)$, then
\begin{align}
& x'Axb'x = (\mu + y)'A(\mu + y)b'(\mu + y) \\
=& (\mu'A\mu + 2\mu'Ay + y'Ay)(b'\mu + b'y) \\
=& \mu'A\mu b'\mu + \mu'A\mu b'y + 2b'\mu\mu'Ay + 2\mu'Ayb'y + b'\mu y'Ay + y'Ayb'y
\end{align}
In my answer to this thread, it is argued that for $y \sim N(0, \Sigma)$, we have
\begin{align}
E[y'Ayb'y] = 0.
\end{align}
In addition,
\begin{align}
& E[\mu'Ayb'y] = E[y'b\mu'Ay] = \tr(b\mu'A\Sigma) = \mu'A\Sigma b,  \\
& E[b'\mu y'Ay] = b'\mu\tr(A\Sigma).
\end{align}
It thus follows that
\begin{align}
& E[x'Axb'x] = \mu'A\mu b'\mu + 2\mu'A\Sigma b + b'\mu\tr(A\Sigma), \\
& \Cov(x'Ax, b'x) = \mu'A\mu b'\mu + 2\mu'A\Sigma b + b'\mu\tr(A\Sigma) - 
(\mu'A\mu + \tr(A\Sigma))b'\mu \\ 
&\phantom{\Cov(x'Ax, b'x)} = 2\mu'A\Sigma b.
\end{align}
Substituting
\begin{align}
& \Var(x'Ax) = 2\tr((A\Sigma)^2) + 4\mu'A\Sigma A\mu, \\
& \Var(b'x) = b'\Sigma b, \\
& \Cov(x'Ax, b'x) = 2\mu'A\Sigma b
\end{align}
into $(1)$, we get
\begin{align}
\Var(z) = 2\tr((A\Sigma)^2) + 4\mu'A\Sigma A\mu + b'\Sigma b + 4\mu'A\Sigma b.
\end{align}
