Here my answer for the case where $K$ divides $N$ :

I consider the intervall $[-1,1]$ instead of $[0,1]$ .

Let $A$ be the $N\times K$ matrix with elements $x_i^{j-1}, 1 \leq i \leq N, 1 \leq j \leq K$

By the Cauchy-Binet theorem the function we want to maximize equals $det(A^T A)$ .

Next I use the result of Fejer 1932 (see http://www.math.technion.ac.il/hat/fpapers/fejerpisa.pdf) :

Let $y_i , 1 \leq i \leq K$ be the zeros of the polynomial $(1-x^2)P'_{N-1}(x)$ where $P_k$ is the $k$-th Legendre polynomial, and let $l_i(x)$ be the fundamental Lagrange interpolating polynomials associated to these points.

Then it holds :

$$ \sum_{i=1}^K l_i(x)^2 \leq 1 $$ for $-1 \leq x \leq 1$ .

Now I follow the paper of Bos, Taylor and Wingate cited in a comment :

Since $$x_i^{j-1} = \sum_{k=1}^K y_k^{j-1} l_k(x_i)$$, we can write $A = B C$, where $B$ has the matrix elements $l_k(x_i), 1 \leq i \leq N, 1 \leq k \leq K$ and $C$ has the matrix elements $y_k^{j-1} , 1 \leq k \leq K, 1 \leq j \leq K$ .

Therefore $$det(A^T A) = det(B^T B) \prod_{1 \leq i < j \leq K} (y_i - y_j)^2$$ .

Now I use Hadamard's inequality, the fact that the geometric mean is less or equal the arithmetic mean and Fejer's inequality and obtain :

$$det(B^T B) \leq \prod_{1 \leq k \leq K} \sum_{i=1}^N l_k(x_i)^2 \leq (\dfrac{1}{K} \sum_{k=1}^K \sum_{i=1}^N l_k(x_i)^2)^K \leq (\dfrac{N}{K})^K$$

Here equality holds iff the square of the euclidean norm of each column vector of $B$ equals $N/K$ and iff they are pairwise orthogonal and this is the case iff $\lbrace x_i : 1 \leq i \leq N\rbrace = \lbrace y_j : 1 \leq j \leq K\rbrace$ and $\vert \lbrace i : x_i = y_j\rbrace\vert = N/K$ for each $j$ (note that equality in Fejer's inequality holds iff $x = y_j$ for a $j$).

Here my answer for the case where $K$ divides $N$ :

I consider the intervall $[-1,1]$ instead of $[0,1]$ .

Let $A$ be the $N\times K$ matrix with elements $x_i^{j-1}, 1 \leq i \leq N, 1 \leq j \leq K$

By the Cauchy-Binet theorem the function we want to maximize equals $det(A^T A)$ .

Next I use the result of Fejer 1932 (see http://www.math.technion.ac.il/hat/fpapers/fejerpisa.pdf) :

Let $y_i , 1 \leq i \leq K$ be the zeros of the polynomial $(1-x^2)P'_{K-1}(x)$ where $P_k$ is the $k$-th Legendre polynomial, and let $l_i(x)$ be the fundamental Lagrange interpolating polynomials associated to these points.

Then it holds :

$$ \sum_{i=1}^K l_i(x)^2 \leq 1 $$ for $-1 \leq x \leq 1$ .

Now I follow the paper of Bos, Taylor and Wingate cited in a comment :

Since $$x_i^{j-1} = \sum_{k=1}^K y_k^{j-1} l_k(x_i)$$, we can write $A = B C$, where $B$ has the matrix elements $l_k(x_i), 1 \leq i \leq N, 1 \leq k \leq K$ and $C$ has the matrix elements $y_k^{j-1} , 1 \leq k \leq K, 1 \leq j \leq K$ .

Therefore $$det(A^T A) = det(B^T B) \prod_{1 \leq i < j \leq K} (y_i - y_j)^2$$ .

Now I use Hadamard's inequality, the fact that the geometric mean is less or equal the arithmetic mean and Fejer's inequality and obtain :

$$det(B^T B) \leq \prod_{1 \leq k \leq K} \sum_{i=1}^N l_k(x_i)^2 \leq (\dfrac{1}{K} \sum_{k=1}^K \sum_{i=1}^N l_k(x_i)^2)^K \leq (\dfrac{N}{K})^K$$

Here equality holds iff the square of the euclidean norm of each column vector of $B$ equals $N/K$ and iff they are pairwise orthogonal and this is the case iff $\lbrace x_i : 1 \leq i \leq N\rbrace = \lbrace y_j : 1 \leq j \leq K\rbrace$ and $\vert \lbrace i : x_i = y_j\rbrace\vert = N/K$ for each $j$ (note that equality in Fejer's inequality holds iff $x = y_j$ for a $j$).