Suppose that $Y_i\sim Exp(\lambda=1)$
We denote by $S$ the sum $Y_1+Y_2+Y_3$.
The joint pdf of the random variables $Y_1, Y_2, S$ is:
$$\begin{array}{l}f_{Y_1, Y_2, S}(y_1, y_2, s)=f_{Y_1, Y_2}(y_1,y_2)f_{S|Y_1=y_1, Y_2=y_2}(s)=\\ =e^{-y_1}e^{-y_2}f_{Y_3}(s-y_1-y_2)=e^{-y_1-y_2}e^{-s+y_1+y_2}=e^{-s}\end{array}$$
for $y_1\geq 0, y_2\geq 0, s\geq y_1+y_2$. Here $f_{S|Y_1=y_1, Y_2=y_2}$ stands for the conditional pdf of the rv S, under conditions $Y_1=y_1$, $Y_2=y_2$, which amounts to the pdf of $Y_3$ evaluated at $s-y_1-y_2$.
In order to get the joint distribution of the random variables $X_1=\displaystyle\frac{Y_1}{S}, X_2=\displaystyle\frac{Y_1+Y_2}{S}, S$ we consider
the transformation $\psi: (y_1, y_2, s)\mapsto (x_1, x_2, s)$, with:
$$
x_1=y_1/s,\quad x_2=(y_1+y_2)/s, \quad s=s$$
Hence $\psi^{-1}(x_1, x_2, s)=(y_1, y_2, s)$ is defined by:
$$y_1=sx_1, \quad y_2=-sx_1+sx_2, \quad s=s$$
and its Jacobian is:
$$\left|\begin{array}{rrr}s&0&x_1\\ -s&s &x_2-x_1\\0&0&1\end{array}\right|=s^2$$
Thus the joint distribution of the variables $X_1, X_2, S$ is:
$g(x_1, x_2, s)=e^{-s}s^2$, for $sx_1, sx_2-sx_1\geq 0$ and $s\geq sx_1+sx_2-sx_1$. But these conditions amount to: $x_1, x_2\geq 0$, $x_1\leq x_2\leq 1$.
To get the joint density distribution, $h_{X_1, X_2}(x_1, x_2)$, of the random variables $X_1,X_2$ we integrate with respect to $s$ the
density $g(x_1, x_2, s)$:
$$
h_{X_1, X_2}(x_1, x_2)=\int_0^\infty s^2 e^{-s}ds=2! \quad x_1, x_2\in[0,1], \quad x_1\leq x_2$$
But this is just the joint distribution of the order statistics, as you stated in your question.
If you start with $Y_k\sim Exp(\lambda)$, $k=1,2,3$, with $\lambda>0$ arbitrary you get that the joint distribution of the variables $X_1, X_2$ coincides with that of the order statistics
$U_{(1)}, U_{(2)}$, where $U_1, U_2\sim Unif[0,1]$, iff $\lambda=1$.
Best Answer
This is about the quantity $$Q:=\int_{{\mathbb R}^2\times{\mathbb R}^2}\exp\bigl(-k|x-y|^2\bigr)\>{\rm d}(x_1,x_2,y_1,y_2)\ .$$ Unfortunately $Q=\infty$. In order to show this we use the transformation $$x_1={u_1+v_1\over2},\quad y_1={v_1-u_1\over2},\quad x_2={u_2+v_2\over2},\quad y_2={v_2-u_2\over2}\ .\tag{1}$$ Then $x-y=(u_1,u_2)$, and the Jacobian of $(1)$ computes to ${1\over4}$. It follows that $$Q={1\over4}\int_{{\mathbb R}^2\times{\mathbb R}^2}\exp\bigl(-k(u_1^2+u_2^2)\bigr)\>{\rm d}(u_1,u_2,v_1,v_2)\ .$$ Here the $u$-integral is fine, but the $v$-integral is $=\infty$.