$\newcommand{\+}{^{\dagger}}
\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
\newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
\newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
\newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
\newcommand{\dd}{{\rm d}}
\newcommand{\down}{\downarrow}
\newcommand{\ds}[1]{\displaystyle{#1}}
\newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
\newcommand{\fermi}{\,{\rm f}}
\newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
\newcommand{\half}{{1 \over 2}}
\newcommand{\ic}{{\rm i}}
\newcommand{\iff}{\Longleftrightarrow}
\newcommand{\imp}{\Longrightarrow}
\newcommand{\isdiv}{\,\left.\right\vert\,}
\newcommand{\ket}[1]{\left\vert #1\right\rangle}
\newcommand{\ol}[1]{\overline{#1}}
\newcommand{\pars}[1]{\left(\, #1 \,\right)}
\newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
\newcommand{\pp}{{\cal P}}
\newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
\newcommand{\sech}{\,{\rm sech}}
\newcommand{\sgn}{\,{\rm sgn}}
\newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
\newcommand{\ul}[1]{\underline{#1}}
\newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
\newcommand{\wt}[1]{\widetilde{#1}}$
\begin{align}
\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\rm P}\pars{\xi}&=
\int_{-\infty}^{\infty}\dd x\,{1 \over \root{2\pi}\sigma}\,
\exp\pars{-\,{x^{2} \over 2\sigma^{2}}}
\\ & \
\int_{-\infty}^{\infty}\dd y\,{1 \over \root{2\pi}\sigma}\,
\exp\pars{-\,{y^{2} \over 2\sigma^{2}}}\
\overbrace{\delta\pars{\xi - xy}}^{\ds{\delta\pars{y - \xi/x} \over \verts{x}}}
\\[3mm]&={1 \over 2\pi\sigma^{2}}\int_{-\infty}^{\infty}
\exp\pars{-\,{1 \over 2\sigma^{2}}\bracks{x^{2} + {\xi^{2} \over x^{2}}}}\,
{\dd x \over \verts{x}}
\\[3mm]&={1 \over \pi\sigma^{2}}\int_{0}^{\infty}
\exp\pars{-\,{1 \over 2\sigma^{2}}\bracks{x^{2} + {\xi^{2} \over x^{2}}}}\,
{\dd x \over x}\tag{1}
\end{align}
$\ds{\delta\pars{x}}$ is the
Dirac Delta Function.
With $\ds{x \equiv A\expo{\theta/2}\,,\quad A > 0\,,\quad\theta \in {\mathbb R}}$:
\begin{align}
{\rm P}\pars{\xi}&={1 \over \pi\sigma^{2}}
\int_{-\infty}^{\infty}
\exp\pars{-\,{1 \over 2\sigma^{2}}
\bracks{A^{2}\expo{\theta} + {\xi^{2} \over A^{2}}\expo{-\theta}}}\,
\pars{A\expo{\theta/2}\,\dd\theta/2 \over A\expo{\theta/2}}
\end{align}
We can choose $\ds{A}$ such that
$\ds{A^{2} = {\xi^{2} \over A^{2}}\quad\imp\quad A = \verts{\xi}^{1/2}}$:
\begin{align}
{\rm P}\pars{\xi}&={1 \over 2\pi\sigma^{2}}
\int_{-\infty}^{\infty}
\exp\pars{-\,{\verts{\xi} \over \sigma^{2}}\cosh\pars{\theta}}\,\dd\theta
={1 \over \pi\sigma^{2}}
\int_{0}^{\infty}
\exp\pars{-\,{\verts{\xi} \over \sigma^{2}}\cosh\pars{\theta}}\,\dd\theta
\end{align}
$$\color{#00f}{\large%
{\rm P}\pars{\xi} = {1 \over \pi\sigma^{2}}\,
{\rm K}_{0}\pars{\verts{\xi} \over \phantom{2}\sigma^{2}}}
$$
where $\ds{{\rm K}_{\nu}\pars{z}}$ is a
Second Kind Bessel Function.
$XX^T$ is Wishart $W_n(\Sigma,1),$ so by a property of the Wishart distribution $(a^TX)^2=a^TXX^Ta$ is $W_1(a^T\Sigma a,1)$. Comparing the pdf of this distribution to that of a chi-square, you find that they are the same up to a scaling.
Best Answer
Lets assume you have $X=(X_1, \dots, X_n)$ a random vector with multinormal distribution with expectation vector $\mu$ and covariance matrix $\Sigma$. We are interested in the quadratic form $Q(X)= X^T A X = \sum \sum a_{ij} X_i X_j$. Define $Y = \Sigma^{-1/2} X$ where we are assuming $\Sigma$ is invertible. Write also $Z=(Y-\Sigma^{-1/2} \mu)$, which will have expectation zero and variance matrix the identity.
Now $$ Q(X) = X^T A X= (Z+\Sigma^{-1/2} \mu)^T \Sigma^{1/2}A\Sigma^{1/2} (Z+\Sigma^{-1/2} \mu). $$ Use the spectral theorem now and write $\Sigma^{1/2}A \Sigma^{1/2} = P^T \Lambda P$ where $P$ is an orthogonal matrix (so that$P P^T=P^T P=I$) and $\Lambda$ is diagonal with positive diagonal elements $\lambda_1, \dotsc, \lambda_n$. Write $U = P Z$ so that $U$ is multivariate normal with identity covariance matrix and expectation zero.
We can compute $$ Q(X) = (Z+\Sigma^{-1/2} \mu)^T \Sigma^{1/2}A\Sigma^{1/2} (Z+\Sigma^{-1/2} \mu) \\ = (Z+\Sigma^{-1/2} \mu)^T P^T \Lambda P (Z+\Sigma^{-1/2} \mu) \\ = (PZ+P\Sigma^{-1/2} \mu)^T \Lambda (PZ+P\Sigma^{-1/2} \mu) \\ = (U+b)^T \Lambda (U+b) $$ where now $b = P \Sigma^{-1/2} \mu $. (There was a small typo in above defs of $U$ and $b$, now corrected.) So: $$ Q(X) = X^T A X = \sum_{j=1}^n \lambda_j (U_j+b_j)^2 $$ In your case, $\mu=0$ so $b=0$ so your quadratic form is a linear combination of independent chi-square variables, each with one degree of freedom. In the general case, we will get a linear combination of independent non-central chisquare variables.
If you want to work numerically with that distribution, there is an CRAN package (that is, package for R) implementing it, called
CompQuadForm.If you want (much) more detail, there is a book dedicated to the topic, Mathai & Provost: "Quadratic forms in random variables".