We will try to show that the numbers $1/4$ and $3/4$ for $1$ point, and the following ones fit in the following scheme, which starts with the
head value $J(0,1)=1$:
$$
1\\
\frac 14\ ,\ \frac 34\\
\frac 2{18}\ ,\ \frac 5{18}\ ,\ \frac {11}{18}\\
\frac 3{48}\ ,\ \frac 7{48}\ ,\ \frac {13}{48}\ ,\ \frac {25}{48}\\
$$
and so on, but the good way to write this is:
$J(n,1)=1/(n+1)^2$ (which will be shown in detail below), and
$$
\boxed{\qquad
(*)\qquad\qquad
J(n,n+1-k) =
\frac 1{n+1}
\left(\frac 1{n+1}+\dots+\frac 1k\right)\ ,
\qquad}
$$
and explicitly we get the following values for the first possible indices,
arranged in a triangle,
$$
\frac 11\left(\frac11\right)
\\
\frac 12\left(\frac 12\right)\ ,\
\frac 12\left(\frac 12+\frac 11\right)
\\
\frac 13\left(\frac 13\right)\ ,\
\frac 13\left(\frac 13+\frac 12\right)\ ,\
\frac 13\left(\frac 13+\frac 12+\frac 11\right)
\\
\frac 14\left(\frac 14\right)\ ,\
\frac 14\left(\frac 14+\frac 13\right)\ ,\
\frac 14\left(\frac 14+\frac 13+\frac 12\right)\ ,\
\frac 14\left(\frac 14+\frac 13+\frac 12+\frac 11\right)
\\
\frac 15\left(\frac 15\right)\ ,\
\frac 15\left(\frac 15+\frac 14\right)\ ,\
\frac 15\left(\frac 15+\frac 14+\frac 13\right)\ ,\
\frac 15\left(\frac 15+\frac 14+\frac 13+\frac 12\right)\ ,\
\frac 15\left(\frac 15+\frac 14+\frac 13+\frac 12+\frac 11\right)
\\
$$
and so on.
The harmonic numbers are thus the hidden tenor.
The reader in hurry can skip some paragraphs and go to the compact solution int he sequel.
The following long and possibly boring exposition gives the sincere solution, showing how it was found.
One observation maybe at this point.
The above numbers are building the triangle
(1) divided by 1
(1, 3) divided by 4
(2, 5, 11) divided by 18
(3, 7, 13, 25) divided by 48
(12, 27, 47, 77, 137) divided by 300
(10, 22, 37, 57, 87, 147) divided by 360
(60, 130, 214, 319, 459, 669, 1089) divided by 2940
(105, 225, 365, 533, 743, 1023, 1443, 2283) divided by 6720
(280, 595, 955, 1375, 1879, 2509, 3349, 4609, 7129) divided by 22680
(252, 532, 847, 1207, 1627, 2131, 2761, 3601, 4861, 7381) divided by 25200
(2520, 5292, 8372, 11837, 15797, 20417, 25961, 32891, 42131, 55991, 83711) divided by 304920
and any row can be obtained from the previous one knowing only its first entry, which is $1/(n+1)^2$. We subtract this first entry from the remaining ones,
and the new data is proportional to the previous row. (Using the harmonic numbers representation this becomes obvious.) For instance, if we subtract $2$, the first entry in $(2,5,11)$ from each other entry, we get $(0,3,9)$, which is proportional to the previous row, $(1,3)$ (after forgetting the first zero entry). Also, starting with $(3,7,13,25)$ we get $(0,4,10,25)$ and $(4,10,22)\sim(2,5,11)$.
First, let us get the theoretical framework.
There is indeed a lot of literature on spacings, here i will use [P] for R. Pyke, Spacings, 1965 . The marked passages below are cited from [P].
I was somehow surprised to see in [P], §4, Construction of Spacings, some models for the Spacings. One of them is as follows, and i will use $n$ (instead of $N$ as in the OP) for the number of points to cut the segment $[0,1]$ into pieces.
The following is extracted from [P], §2.1 and §4.1.
§2.1.
Let $X=(X_1,X_2,\dots,X_n)$ be a tuple/family of independent r.v.'s on $(0,1)$. The density function $f_X$ of $X$ in $x=(x_1,x_2,\dots,x_n)$ is then
$$f_X(x) =
\begin{cases}
1 &\text{ if }x\in[0,1]^n\ ,\\
0 &\text{else .}
\end{cases}
$$
Let $U=(U_1,U_2,\dots,U_n)$ be the order statistics of $X$,
obtained by arranging the components of $X$ in increasing order. The density function $f_U$ of $U$ in $u=(u_1,u_2,\dots,u_n)$ is then
$$f_U(u) =
\begin{cases}
n! &\text{ if }0\le u_1\le u_2\le\dots\le u_n\le 1\ ,\\
0 &\text{else .}
\end{cases}
$$
Set $U_0=0$, $U_{n+1}=1$,
and let $D=(D_1,D_2,\dots,D_{n+1})=(U_1-U_0,U_2-U_1,\dots,U_{n+1}-U_n)$ be the spacings statistics of $U$.
Then $D$ has a singular distribution, since
$D_1+D_2+\dots+D_{n+1}=1$, but when restricted to this hyperplane has (w.r.t. the $n$--dimensional Lebesgue mass) the density in
$d=(d_1,d_2,\dots,d_{n+1})$ with $\sum d:=d_1+d_2+\dots+d_{n+1}=1$
$$f_D(d) =
\begin{cases}
n! &\text{ if }0\le d_1, d_2,\dots,d_{n+1}\text{ and }\sum d=1\ ,\\
0 &\text{else .}
\end{cases}
$$
Therefore most probability statements about $D$ may, theoretically,
be obtained by computing the volume of a subset of $\sum d=1$,
or equivalently of the simplex $0\le d_1,\dots,d_n$ and $d_1+\dots+d_n\le 1$. (Only $n$ variables.)
and also
§4.1
Let $Y_1,Y_2,\dots,Y_{n+1}$ be independent exponential random variables (r.v.) with mean $1$. Set
$$
S = Y_1+Y_2+\dots+Y_{n+1}
$$
for their sum and
$$
D_j = \frac {Y_j}S\ ,\qquad 1\le j\le n+1\ .
$$
Then $(D_1,D_2,\dots,D_{n+1})$ is distributed as the set of $(n+1)$ spacings determined by $n$ independent uniform r.v. .
This means for us that we have to compute for a fixed $n$ the integrals (in the variable $x$ instead of the above $d$, so that i can write without a drunken syndrom that $dx$ as usual)
$$
\begin{aligned}
J(n,k)
&:=\int_{\substack{x\in[0,1]^n}\\\sum x\le 1}(\text{$k$. value after ordering $(x_1,x_2,\dots, x_n,1-\sum x)$})\; dx\\
&=\int_{\substack{y\in[0,1]^{n+1}}\\\sum y= 1}(\text{$k$. value after ordering $(y_1,y_2,\dots, y_n,y_{n+1})$})\; dy\\
&=(n+1)!
\int_{\substack{0\le y_1\le y_2\le \dots\le y_{n+1}\le 1\\
y\in[0,1]^{n+1}\\
\sum y= 1}}\qquad y_k\; dy
\ .
\end{aligned}
$$
Here, $x=(x_1,x_2,\dots,x_n)$ has "only" $n$ components, and $y$ has $(n+1)$ components. We will use a "projected version" of it.
It is always a good idea to let the computer repeat some Monte-Carlo simulations.
The following code:
def streichholzBrechStatistik_np(n, size=10**7):
v = np.random.random_sample(n*size)
zeeros = np.zeros(size, dtype='float').reshape( (size, 1) )
ooones = np.ones (size, dtype='float').reshape( (size, 1) )
sample = np.column_stack( (zeeros, v.reshape(size, n), ooones) )
sample . sort()
parts = sample[ : , 1:] - sample[ : , :-1 ]
parts . sort()
return sum(parts) / size
for n in (1, 2, 3, 4, 5, 6):
print(n, streichholzBrechStatistik_np(n))
It gave me this time:
1 [0.25008897 0.74991103]
2 [0.11113685 0.2777319 0.61113125]
3 [0.06250226 0.14584965 0.27083017 0.52081792]
4 [0.04000254 0.08999392 0.15666303 0.25666772 0.4566728 ]
5 [0.02778191 0.06111756 0.10280426 0.15832957 0.24164264 0.40832405]
6 [0.0204181 0.04421353 0.07277986 0.10849184 0.15611271 0.22754258 0.37044138]
It remains to guess the pattern.
In the following we will compute detailed the values $J(n,1)$.
J(n,1) in the case n = 1 :
$$
J(1,1)
=
2!\int_{\substack{0\le y_1\le y_2\le 1\\y_1+y_2=1}} y_1\; dy
=
2!\int_{y_1\in[0,1/2]}y_1\; dy_1
=2!\left[\ \frac 12y_1^2\ \right]_0^{1/2}
=\frac 14\ ,
$$
and $J(1,2)=1-J(1,1)=\frac 34$.
J(n,1) in the case n = 2 :
It is maybe good to compare the points below with those in a barycentric subdivision of the $2$--simplex
as in the raw ascii picture
(1:0:0)
*
/|\
/ | \
/ | \
/ | \
/ * (1:1:1)
/ | \
/ | \
*-------*-------*
(0:1:0) (0:1:1) (0:0:1)
and note that the part with $0\le y_1\le y_2\le y_3\le 1$ corresponds to
one of the six triangles delimited by the barycentric points (vertices,
mid points of sides, and centroid). If the solution is given strictly in terms of the triangle, than it may be that we have a case separation, which may be hard to generalize. But let us give this computation, since it is suggerated in the OP.
The following integral is built w.r.t. the $2$--simplex part:
$0\le y_1\le y_2\le y_3\le 1$ and $y_1+y_2+y_3=1$. We may want to have only two
inequalites in the variables $y_1,y_2$. So we have to get rid of $y_3$.
I prefer to project onto the $Oy_1y_2$ plane. Then we have a "polytope", and its vertices are obtained by solving a system of two equations, obtained by forgetting each one among the three boundary conditions: $0=y_1$, and $y_1=y_2$, and $y_2=1-y_1-y_2$. The points we obtain are $\left(frac 12,\frac 12\right) $,
and $(0,1)$, and $(0,0)$. In a similar way, if we insist to fix $y_1=a$, then this value $a$ gets "bounded" from $1=y_1+y_2+y_3$, so $y_2$ varies between
$a$ and $\frac 12(1-3a)$. (The last $\frac 12$ corresponds to $y_2=y_3$ before we project onto $Oy_1y_2$.)
(The three points correspond to $(0:0:1), (0:1:1)$, and $(1:1:1)$.)
Then we have
$$
\begin{aligned}
J(2,1)
&=
\frac
{\displaystyle\int_0^{1/3}y_1\; dy_1\int_{y_1}^{(1-y_1)/2}dy_2}
{\displaystyle\int_0^{1/3}1\; dy_1\int_{y_1}^{(1-y_1)/2}dy_2}
\\
&=
\frac
{\displaystyle\int_0^{1/3}a\; da\cdot\frac 12(1-3a)}
{\displaystyle\int_0^{1/3}1\; da\cdot\frac 12(1-3a)}
\\
&=\frac 19\ .
\end{aligned}
$$
J(n,1) in the case n = 3 :
We calculate only $J(3,1)$. We fix some $y_1=a$, and of course $a\in[0,\ 1/4]$.
Which is the section of $y_1=a$ in the already projected polytope:
$$
0\le y_1\le y_2\le y_3\le 1-y_1-y_2-y_3\ ?
$$
For the fixed $y_1=a$, the three boundary conditions in
$y_2,y_3$ are $a=y_2$, and $y_2=y_3$, and $y_3=1-a-y_2-y_3$. The points we obtain are
$$
\begin{aligned}
\text{For }& y_2=y_3\ ,
\ y_3=1-a-y_2-y_3 &&\to \left(\frac 13(1-a),\ \frac 13(1-a)\right)\ ,
\\
\text{For }& a=y_2\ ,
\ y_3=1-a-y_2-y_3 &&\to \left(a,\ \frac 12(1-2a)\right)\ ,
\\
\text{For }& a=y_2\ ,
\ y_2=y_3 &&\to \left(a,\ a\right)\ .
\end{aligned}
$$
The corresponding $2$-volume is related to the determinant
$$
\begin{vmatrix}
\frac 13(1-a) & \frac 13(1-a) & 1\\
a & \frac 12(1-2a) & 1\\
a & a & 1
\end{vmatrix}
=
\begin{vmatrix}
\frac 13(1-a)-a & \frac 13(1-a)-a & 0\\
0 & \frac 12(1-2a)-a & 0 \\
a & a & 1
\end{vmatrix}
=
\begin{vmatrix}
\frac 13(1-4a)-a & * & 0\\
0 & \frac 12(1-4a) & 0 \\
a & a & 1
\end{vmatrix}
=
\begin{vmatrix}
\frac 13(1-4a)-a & * \\
0 & \frac 12(1-4a)
\end{vmatrix}
=
\frac 1{3!}(1-4a)^2
\ .
$$
A simpler way to argue would be, that $a$ is extracted from ($y_1$ and)
$y_2,y_3,y_4$, so the scaled $2$--simplex has the length scale proportional to
$(1-4a)$, so the $2$-volume is proportional to $(1-4a)^2$.
Then we have
$$
J(3,1)
=
\frac
{\displaystyle\int_0^{1/4}a\; da\cdot\frac 1{3!}(1-4a)^2}
{\displaystyle\int_0^{1/4}1\; da\cdot\frac 1{3!}(1-4a)^2}
=\frac 1{16}
\ .
$$
J(n,1) in the general case:
The short calculus above will be done for a general $n$. Again we set $y=a$.
Then the corresponding $(n-1)$-simplex has a volume $\displaystyle\sim (1-(n+1)a)^{n-1}$,
so we can immediately compute:
$$
\begin{aligned}
J(n,1)
&=
\frac
{\displaystyle\int_0^{1/(n+1)}a\; da\cdot(\dots)(1-(n+1)a)^{n-1}}
{\displaystyle\int_0^{1/(n+1)}1\; da\cdot(\dots)(1-(n+1)a)^{n-1}}
\\
&\qquad\text{Substitution: }t=(n+1)a
\\
&=
\frac
{\displaystyle\int_0^1\frac 1{n+1}t\; dt\cdot (1-t)^{n-1}}
{\displaystyle\int_0^1 1\; dt\cdot (1-t)^{n-1}}
=
\frac
{\displaystyle\int_0^1\frac 1{n+1}(1-u)\; du\cdot u^{n-1}}
{\displaystyle\int_0^1 1\; du\cdot u^{n-1}}
\\
&=
\frac
{\displaystyle\frac 1{n+1}\left[\frac 1n-\frac 1{n+1}\right]}
{\displaystyle\frac 1n}
=\frac 1{(n+1)^2}
\ .
\end{aligned}
$$
Compact solution:
From these computations we can conclude with the observation that the proportion
$$
J(n+1,2)-J(n+1,1) \ :\
J(n+1,3)-J(n+1,1) \ :\
J(n+1,4)-J(n+1,1) \ :\
\dots
$$
is the same one as
$$
J(n,1) \ :\
J(n,2) \ :\
J(n,3) \ :\
\dots
$$
(using marginal / conditional densities), which is more or less the same
as computing inductively the integrals (with the same technical ingredients as it was the case for $J(n.1)$), so let us compare for some $k>1$ (and giving the sense with the convention $y_k=0$, $J(n,0)=0$ in case we have to deal with some $k-2$, or considering this case separately)
$$
\begin{aligned}
&J(n+1,k)-J(n+1,k-1)
\\
&\qquad
=
\frac
{\displaystyle
{\color{blue}\int}\int\dots\int_D
{\color{blue}{dy_1}}\; dy_2\;dy_3\; \dots\; dy_k\;
\cdot\; {\color{blue}{(y_{k}-y_{k-1})}}\;\cdot(1{\color{blue}-y_1}-y_2-\dots-y_k)^{(n+2)\ -k\ -1}
}
{\displaystyle
{\color{blue}\int}\int\dots\int_D
{\color{blue}{dy_1}}\; dy_2\;dy_3\; \dots\; dy_k\;
\cdot\; 1\;\cdot(1{\color{blue}{-y_1}}-y_2-\dots-y_k)^{(n+2)\ -k\ -1}
}
\\[2mm]
&\qquad\qquad\text{ and}
\\[2mm]
&J(n,k-1)-J(n,k-2)
\\
&\qquad
=\frac
{\displaystyle
\int\dots\int_E
dz_2\;dz_3\; \dots\; dz_k\;
\cdot\; \color{red}{z_k-z_{k+1})}\;\cdot(1-z_2-\dots-z_k)^{(n+1)\ -(k-1)\ -1}
}
{\displaystyle
\int\dots\int_E
dz_2\;dz_3\; \dots\; dz_k\;
\cdot\; 1\;\cdot(1-z_2-\dots-z_k)^{(n+1)\ -(k-1)\ -1}
}
\\[2mm]
&\qquad\qquad\text{ taken on}
\\[2mm]
D&=D(n+1,k)=
\Big\{\ y\ :
{\color{blue}{0\le y_1}}
\le y_2\le y_3\le\dots\le y_k
\le
\frac
{1{\color{blue}{-y_1}}-y_2-y_3-\dots-y_k}
{\text{number of remained variables}} \
\Big\}\ ,
\\
E&=D(n,k-1)=
\Big\{\ z\ :
0
\le z_2\le z_3\le\dots\le z_k
\le
\frac
{1-z_2-z_3-\dots-z_k}
{\text{number of remained variables}}
\ \Big\}
\ .
\end{aligned}
$$
(Note that $z_1$ is missing in the last difference, so that the function to be integrated is involving the same indices,
$k$, and $k-1$.)
We use now the following substitution which passes from $D$ to $E$,
from the variable $y\in D$ to the variable $z\in E$.
Consider $y_1$ from $y\in D$ as a parameter and separate the integral on $D$ (Fubini)
as an integral in $y_1$ for $y_1\in[0,\ 1/(n+2)$, and an integral in $y'=(y_2,y_3,\dots,y_k)$. This parameter, $y'$, takes values in a range, that is constrained to $y_1\le y_2\le \dots\le y_k$, and to the fact that it can be extended to some $(y_1,y',y'')=(y_1,y_2,\dots,y_k,y_{k+1},\dots, y_{n+2})$ using "remained variables" with sum of the components equal to one. The variables in $y''$ are each at least $y_k$, which gives the upper variation boundary for $y_k$.
The number of the "remained variables" is the same for all integrals above.
We pass from $y'$ to $z$ in the $y$-integrals via the linear (diagonal) substitution
$$
\begin{aligned}
y_2 &= (1-y_1)z_2\ ,\\
y_3 &= (1-y_1)z_2\ ,\\
&\qquad\qquad\text{ and so on up to }\\
y_k &= (1-y_1)z_2\ .
\end{aligned}
$$
The $y$-integrals can now be separated,
as an integral in $y_1\in[0,\ 1/(n+2)]$, and one in the new variable $z\in E$.
The integral in the numerator introduces a proportionality factor,
$$
\int_0^{1/(n+2)}
dy_1\;{\color{blue}{(1-y_1)}}\;\cdot(1{\color{blue}-y_1})^{(n+2)\ -k\ -1}
\ .
$$
Inductively on $k$ we obtain the claimed values in $(*)$.
