PDF of area of triangle with normally-distributed coordinates in any dimensions 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.
 A: 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.

the area of the parallelogram defined by ${\bf a},\, {\bf b}$.

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.

Same as per the previous post, it remains to work out the case in
which the parallelogram is translated from the origin.
The simulations indicates that it is  $\Gamma _{\,m - 1,\,1} (c/ \sqrt{3}) =\Gamma _{\,m - 1,\,\sqrt{3}} (c) $
but could not reach yet to prove that.
This would give in fact an expected value $E\left[ c \right] = \left( {m - 1} \right)\sqrt 3 $
for the parallelogram, which corresponds to the value you cite for the triangle,
and confirms your fitting formula.

