Question
What is the probability distribution function (PDF) of the absolute area of a triangle with normally-distributed coordinates in $\mathbb{R}^m$ $(m \in \mathbb{N}, m\ge2)$ ? A conjecture is given that can be proved or might help to find the correct solution.
The triangle vertices in $\mathbb{R}^m$ are $$ \mathbf{\mathrm{X}_1} =(x_1^1,\ldots, x_1^m),\;\; \mathbf{\mathrm{X}_2}=(x_2^1,\ldots,x_2^m),\;\; \mathbf{\mathrm{X}_3}=(x_3^1,\ldots,x_3^m)$$ where $x_i^j$ are independent standard normal distributed variables $$x_i^j\sim\mathcal{N}(0,1)$$ The non-oriented area $A$ of a random triangle instance in $\mathbb{R}^m$ is $$A=\|\mathbf{\mathrm{X}_1}-\mathbf{\mathrm{X}_3}\|\, \|\mathbf{\mathrm{X}_2}-\mathbf{\mathrm{X}_3}\| \frac{\sin\alpha}{2}\tag{1}$$ where $\|\cdot\|$ is the Euclidean norm and $\alpha$ is the angle at $\mathbf{\mathrm{X}_3}$. The expectation value of $A$ is $$\mathbb{E}[A]=\frac{\sqrt{3}}{2}\left(m-1\right)\tag{2}$$ A proof of eq.(2) can be found in this Cross Validated post.
Conjecture for PDF
The probability distribution of empirical data of $A$ for any tested $m$ can be fitted quite well with the function $$f(A)=\frac{k^{m-1}} {A \mathrm{e}^k (m-2)! }\, \, \text{ with} \, \ k=\frac{2A}{\sqrt{3}} \tag{3}$$
that fulfills also the conditions of a PDF $$\int_0^\infty f(A) \mathrm{d}A=1\ \,\, \text{and}\ \int_0^\infty A f(A)\mathrm{d}A=\mathbb{E}[A]$$
Experimental data is indistinguishable from eq.(3) but a proof is missing.
Related question
In a similar case the PDF for the volume of a tetrahedron in $\mathbb{R}^3$ was tried to solve in this Math SE post but no full solution was given so far.
Concerning a triangle with a vertex at the origin , thus a parallelogram defined by two vectors, or otherwise the modulus of their cross-product, an answer for the case 3D has been already provided in my answer to the related post on the PDF of the volume of a tetrahedron, which is $$ \eqalign{ & \Pr \left( {\left| {{\bf a} \times {\bf b}} \right| \in \left[ {c,c + dc} \right)} \right) = p_c (c)dc\quad \left| \matrix{ \,{\bf a},{\bf b} \sim \left( {{\cal N}_1 \,} \right)^{\,3} \hfill \cr \;0 \le c \hfill \cr} \right. = \cr & = {2 \over {\sqrt {2\pi } }}cdc \int_{v = 0}^\infty {e^{\, - \,{1 \over 2}\left( {v^{\,2} + c^{\,2} /v^{\,2} } \right)} \,dv} = \cr & = c\,e^{\, - \,\,c} dc \quad \Rightarrow \quad p_c (c) = \Gamma _{\,2,\,1} (c) \cr} $$
The derivation employed there can be extended to $m$-D as follows.
The density in the $m$-space reads as $$ \eqalign{ & {{p_p} _m} ({\bf V})dV = \,\left( {{\cal N}_{\sigma ^{\,2} } \,(x_{\,1} ){\cal N}_{\sigma ^{\,2} } \,(x_{\,1} ) \cdots {\cal N}_{\sigma ^{\,2} } \,(x_{\,m} )} \right)\,dx_{\,1} \,dx_{\,2} \, \cdots dx_{\,m} = \cr & = \left( {{1 \over {\sigma \sqrt {2\pi } }}} \right)^{\,m} e^{\, - \,{1 \over 2} \left( {{r \over \sigma }} \right)^{\,2} } \,dA_{\,m - 1} (r)\,dr = \cr & = \left( {{1 \over {\sqrt {2\pi } }}} \right)^{\,m} e^{\, - \,{1 \over 2}\left( {{r \over \sigma }} \right)^{\,2} } \,dA_{\,m - 1} (r/\sigma )\,d(r/\sigma ) = \cr & = \left( {{1 \over {\sqrt {2\pi } }}} \right)^{\,m} e^{\, - \,{1 \over 2}\rho ^{\,2} } \,dA_{\,m - 1} (\rho )\,d\rho \cr} $$ with the surface of a $m-1$-sphere (in a $m$- space) being $$ A_{\,m - 1} (\rho ) = {{2\pi ^{\,{m \over 2}} } \over {\Gamma \left( {{m \over 2}} \right)}} \rho ^{\,m - 1} $$ and $dA_{m-1}(r)$indicating an elementary surface on the sphere of constant radius $r$.
So, taking a unitary variance and using $r$ for simplicity, upon integrating over the surface at constant $r$ we get $$ \eqalign{ &{{ p_r} _m} (r)dr = \left( {{1 \over {\sqrt {2\pi } }}} \right)^{\,m} {{2\pi ^{\,{m \over 2}} } \over {\Gamma \left( {{m \over 2}} \right)}} r^{\,m - 1} e^{\, - \,{1 \over 2}r^{\,2} } \,\,dr = \cr & = {1 \over {2^{{m \over 2} - 1} \Gamma \left( {{m \over 2}} \right)}} r^{\,m - 1} e^{\, - \,{1 \over 2}r^{\,2} } \,\,dr = \chi _m \,(r)\,dr \cr} $$ as it shall be.
Let's pass to the "cross-product" of two vectors ${\bf a}, \, {\bf b}$ having such a spherical distribution.
In the $m$-D space we generalize the cross product by the wedge product whose modulus can be written as $$ \eqalign{ & \left| {\,{\bf a} \wedge {\bf b}\,} \right| = \left| {\,{\bf a}\,} \right|\left| {\,{\bf b}\,} \right|\sin \left( {\angle {\bf a},{\bf b}} \right) = \left| {\,{\bf a}\,} \right|\left| {\,{\bf b}_{\, \bot {\bf a}} \,} \right| \cr & = a\sqrt {b^{\,2} - \left( {{\bf b} \cdot {{\bf a} \over a}} \right)^{\,2} } = \sqrt {a^{\,2} b^{\,2} - \left( {{\bf b} \cdot {\bf a}} \right)^{\,2} } \cr} $$ with a clear meaning of the symbols, i.e.
We fix the first to have modulus $\left[ {a,a + da} \right)$.
Then all the vectors which have a contant wedge-product modulus $\left[{c, c+dc}\right)$ with that will be those which end into a cylindrical circular shell around $\bf a$, with radius $\left[ {c/a, c/a+dc/a}\right)$.
In the axial direction they are ${\mathcal N}_1$ distributed, and integrating over that their probability sums to
$$ \chi _{m - 1} \,(c/a)\,dc/a $$ Therefore we obtain $$ \bbox[lightyellow] { \eqalign{ & \Pr \left( {\left| {{\bf a} \wedge {\bf b}} \right| \in \left[ {c,c + dc} \right)} \right) = {{p_c} _m} (c)dc\quad \left| \matrix{ \,{\bf a},{\bf b} \sim \left( {{\cal N}_1 \,} \right)^{\,m} \hfill \cr \;0 \le c \hfill \cr} \right. = \cr & = \int\limits_a {\chi _m \,(a)\,da\,\chi _{m - 1} \,(c/a)\,dc/a\,} = \cr & = \int_{a = 0}^\infty {{{a^{\,m - 1} } \over {2^{\,m/2 - 1} \Gamma \left( {{m \over 2}} \right)}}e^{\, - \,a^{\,2} /2} \,da {{\left( {c/a} \right)^{\,m - 2} } \over {2^{\,m/2 - 3/2} \Gamma \left( {{{m - 1} \over 2}} \right)}} e^{\, - \,\left( {c/a} \right)^{\,2} /2} \,dc/a\,} = \cr & = dc{{c^{\,m - 2} } \over {2^{\,m - 5/2} \Gamma \left( {{m \over 2}} \right) \Gamma \left( {{{m - 1} \over 2}} \right)}} \int_{a = 0}^\infty {\,e^{\, - \,{1 \over 2}\left( {a^{\,2} + {{c^{\,2} } \over {a^{\,2} }}} \right)} da\,} = \cr & = dc{{c^{\,m - 2} } \over {2^{\,m - 5/2} \Gamma \left( {{m \over 2}} \right) \Gamma \left( {{{m - 1} \over 2}} \right)}}\sqrt {{\pi \over 2}} \,e^{\, - \,\,\sqrt {c^{\,2} } } = \cr & = dc{{c^{\,m - 2} } \over {{{2^{\,2\left( {{{m - 1} \over 2}} \right) - 1} } \over {\sqrt \pi }} \Gamma \left( {{{m - 1} \over 2} + {1 \over 2}} \right)\Gamma \left( {{{m - 1} \over 2}} \right)}}\, e^{\, - \,\,\sqrt {c^{\,2} } } = \cr & = dc{{c^{\,m - 2} } \over {\;\Gamma \left( {m - 1} \right)}}\,e^{\, - \,\,c} = \Gamma _{\,m - 1,\,1} (c)\,dc \cr} }$$ which, as far as the simulation on my computer can go, checks.