Determinant of $a_{i,j}=(x_i+y_j)^k$
How can I find the determinant of the matrix $A\in\mathcal{M}_n(\mathbb{R})$ with coefficients $a_{i,j}=(x_i+y_j)^k,k<n$ ?
All the $x_u,y_u$ are real numbers.
Derivating won't help, and I didn't find any good way to simplify the problem.
Result
For $k<n$,
$$\det(x_i+y_j)^{k}=\mathbb{1}_{k=n-1} (-1)^{\frac{n(n-1)}{2}} \prod_{k=1}^{n-1}k^{2k-n} \prod_{i_1 < i_2}(x_{i_2}-x_{i_1})\prod_{j_1 < j_2}(y_{j_2}-y_{j_1}) $$
Proof
Fact 1
$$\det(x_i+y_j)^{k}=\mathbb{1}_{k=n-1} Const_n \prod_{i_1 < i_2}(x_{i_2}-x_{i_1})\prod_{j_1 < j_2}(y_{j_2}-y_{j_1}) $$
where $Const_n$ is a real number that only depends on $n$ and does not depend on $x$s and $y$s.
Proof 1
If $\exists i_1, i_2: x_{i_1}=x_{i_2}$, then our determinant is $0$ because 2 rows will be the identical.
Likewise, if $\exists j_1, j_2: y_{j_1}=y_{j_2}$, then our determinant is $0$ because 2 columns will be identical.
Thus our determinant, which is clearly a polynomial, should be
$$\prod_{i_1 < i_2}(x_{i_2}-x_{i_1})\prod_{j_1 < j_2}(y_{j_2}-y_{j_1})$$ up to a factor of some polynomial in $x,y$, i.e.
$$\det(x_i+y_j)^{n-1}=SomePolynomial(x_1,\dots, x_n, y_1, \dots, y_n)\prod_{i_1 < i_2}(x_{i_2}-x_{i_1})\prod_{j_1 < j_2}(y_{j_2}-y_{j_1})$$
Note that $\prod_{i_1 < i_2}(x_{i_2}-x_{i_1})\prod_{j_1 < j_2}(y_{j_2}-y_{j_1})$ is a homogeneous polynomial of degree $n(n-1)$, which puts a lower bound on the degree of the determinant.
On the other hand, if you try and write the determinant as a sum of products over all permutations with correct sign (definition), the result will be a homogeneous polynomial of degree $nk$.
Since the total degrees of these two polynomials must agree, we must have the inequality $n(n-1)\leq nk\implies n-1 \leq k$. But $k<n$ by assumption, so either $k=n-1$ or no such polynomial exists; in the latter case we conclude that the determinant vanishes identically.
If $k=n-1$, the leading factor must have degree zero i.e. the two polynomials agree up to a constant, and this constant only depends on $n$:
$$\det(x_i+y_j)^{n-1}=Const_n\prod_{i_1 < i_2}(x_{i_2}-x_{i_1})\prod_{j_1 < j_2}(y_{j_2}-y_{j_1})$$
Fact 2
$$Const_n = (-1)^{\frac{n(n-1)}{2}} \prod_{k=1}^{n-1}k^{2k-n}$$
Proof 2
If we compute the determinant according to the definition, then expand all $(x_i+y_j)^{n-1}$, we will have a linear combination of $$\prod_{i=1}^{n} x_i^{p_i} \prod_{j=1}^{n} y_j^{q_j}$$ where $\sum_{i=1}^n p_i + \sum_{j=1}^n q_j = (n-1)n$.
The contributions of element $(x_i+y_j)^{n-1}$ are of form $$\binom{n-1}{m} x_i^{m} y_j^{n-1-m} \ \forall m \in \{0,1,\dots,n-1\}$$
Thus, for a given $\prod_{i=1}^{n} x_i^{p_i} \prod_{j=1}^{n} y_j^{q_j}$, $$\exists i^* \in \{1,2,\dots,n\} : \#\{i | p_i=p_{i^*}\}=1 \implies \exists ! j^* \in \{1,2,\dots,n\} : (p_{i^*} + q_{j^*}) = n-1$$ In other words, if among the values of powers of $x$s there is a unique one, with index $i^*$, then there will be one and only one value of $q$s equal to $n-1-p_{i^*}$; we denote its index as $j^*$. This means that one of the contributors to this permutation is part of element $(x_{i^*} + y_{j^*})^{n-1}$, $\binom{n-1}{p_{i^*}} x_{i^*}^{p_{i^*}} y_{j^*}^{n-1-p_{i^*}}$
This implies that, for the terms $\prod_{i=1}^{n} x_i^{p_i} \prod_{j=1}^{n} y_j^{q_j}$ where all $p$s are unique, there is only one permutation that contributes to it. To identify this permutation, we would pair up $p$s and $q$s (all $p$s are unique, all $q$s are unique, they all belong to $\{0,1,\dots,n-1\}$, there are $n$ of each, so pairng up exists and is unique). The indices of the pairs will be the indices of the elements of $(x_i+y_j) ^{n-1}$ that contributes to the permutation, the powers will gives use the coefficients of the contribution of each element.
We will apply this fact to one of the contributions of the main diagonal,
$$x_1^{n-1} x_2^{n-2} y_2 x_3^{n-3} y_3^2 \dots x_{n-2}^{2} y_{n-2}^{n-3} x_{n-1} y_{n-1}^{n-2} y_{n}^{n-1}$$
If we the the definition of the determinant, the contribution this term receives (it will come with "+" sign), is
$\prod_{k=0}^{n-1} \binom{n-1}{k} = \frac{\left((n-1)!\right)^n}{\left(\prod_{k=0}^{n-1} k!\right)^2} = \frac{\prod_{k=1}^{n-1}k^n}{\prod_{k=1}^{n-1}k^{2(n-k)}}=\prod_{k=1}^{n-1}k^{2k-n}$
This contribution in terms of $Const_n$ is $$(-1)^{\frac{n(n-1)}{2}}Const_n$$ (all $y$s come with "+" sign, all $x$s come with "-") sign.
Thus $$Const_n = (-1)^{\frac{n(n-1)}{2}} \prod_{k=1}^{n-1}k^{2k-n}$$
Not sure if my proof is clear, but I fear that adding more details will make it even more obscure. Questions welcome!
Here is an answer for the $(k,n)=(2,3)$ case which, as I will discuss at the end, can be readily generalized and made more elegant. Observe that the product of Vandermonde matrices
$$\begin{pmatrix}1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ 1 & x_3 & x_3^2 \end{pmatrix} \begin{pmatrix}y_1^2 & y_2^2 & y_3^2 \\ y_1 & y_2 & y_3 \\ 1 & 1 & 1 \end{pmatrix} =\begin{pmatrix} y_1^2 +x_1 y_1+x_1^2 & y_2^2 +x_1 y_2+x_1^2 & y_3^2 +x_1 y_3+x_1^2 \\ y_1^2 +x_2 y_1+x_2^2 & y_2^2 +x_2 y_2+x_2^2 & y_3^2 +x_2 y_3+x_2^2\\ y_1^2 +x_3 y_1+x_3^2 & y_2^2 +x_3 y_2+x_3^2 & y_3^2 +x_3 y_3+x_3^2 \end{pmatrix}$$
is very nearly off the desired form, with the only incongruity being that the $x_i y_j$ terms should have a coefficient of $2$. This is remedied by inserting a diagonal matrix $D=\text{diag}(1,2,1)$ between the matrices above, thereby multiplying the row of the second Vandermonde matrix by $2$. The elements of the resulting matrix are then indeed of the form $a_{ij}=(x_i+y_j)^2$.
So our matrix is the product of two Vandermonde matrices and a diagonal matrix. Recalling the form of the determinant of a Vandermonde matrix, taking the determinant thus yields the polynomial
$$-2(x_2-x_1)(x_3-x_1)(x_3-x_2)\cdot (y_2-y_1)(y_3-y_1)(y_3-y_2).$$
with a minus sign owing to the orientation of the second Vandermonde matrix.
This is a fairly nice answer, and one imagines the following generalization: For the case of $k=n-1$, the matrix is obtained by multiplying two Vandermonde matrices with a diagonal matrix of binomial coefficients $D_{ij}=\delta_{ij}\binom{n-1}{i}$. The resulting determinantal polynomial would thus be
$$-\prod_{k=1}^{n-1}\binom{n-1}{k}\prod_{1\leq j\leq k\leq n}(x_k-x_j)(y_k-y_j).$$
I imagine an elegant proof of this is possible.
Here is a refinement/generalization of the approach in my first answer. First, note that by way of the binomial theorem we may express $(x+y)^k$ for given $x,y$ in bilinear form:
\begin{align} (x+y)^k =\sum_{l=0}^k\binom{k}{l}x^l y^{k-l} &=\begin{pmatrix}1 & x & \cdots & x^{k-1} & x^k \end{pmatrix}\\ &\quad \cdot \text{diag}\left(1,k,\binom k2,\cdots,k,1 \right)\\ &\quad \cdot \begin{pmatrix}y^k & y^{k-1} & \cdots & y & 1\end{pmatrix}^T \end{align}
We can write this compactly as $(x+y)^k = \mathbf{x}^T (DR) \mathbf{y} $ with $(k+1)\times 1$ column vectors $(\mathbf{x})_{0\leq i \leq k}=x^i$, $(\mathbf{y})_{0\leq i \leq k}=y^i$, a $(k+1)\times (k+1)$ diagonal matrix $D=\text{diag}\left[\binom{k}{i}\right]_{0\leq i \leq k}$, and $R$ the $(k+1)\times (k+1)$ order-reversing permutation matrix $(R)_{0\leq i,j\leq k}=\delta_{i,k-j}$.
Applying this formula element-wise to the problem of interest, we have $a_{ij}=(x_i+y_j)^k=\mathbf{x}_i^T (DR)\mathbf{y}_j$ and so $A=X^T( DR )Y$ where we have introduced the $(k+1)\times n$ rectangular matrices $$X=\begin{pmatrix}\mathbf{x}_1 & \mathbf{x}_2 & \cdots & \mathbf{x}_n\end{pmatrix},\\ Y=\begin{pmatrix}\mathbf{y}_1 & \mathbf{y}_2 & \cdots & \mathbf{y}_n\end{pmatrix}.$$
Writing this as $X^T(DRY)$, we have reduced our problem to that of finding the determinant of a product of $(k+1)\times n$ and $n\times (k+1)$ rectangular matrices. There are three cases:
If $k+1 > n$, then the LH matrix has more rows than columns; but this is ruled out by the assumption $k<n$. (I'm inclined to post a new question for this case.)
If $k+1<n$, then the LH matrix has more columns than rows. Consequently the resulting matrix cannot have full rank (specifically, it has rank at most $n$) and the determinant vanishes.
If $k+1=n$, then all the matrices are square and the determinant of the product is the product of the determinants. For $X^T$ and $Y$, we observe that these are (up to a matrix transpose) Vandermonde matrices over the variables $\{x_{1\leq i \leq n}\}$ and $\{y_{1\leq i \leq n}\}$; the corresponding Vandermonde determinants are $\prod_{1\leq i< j\leq n}(x_j-x_i)$ and $\prod_{1\leq i< j\leq n}(y_j-y_i)$. The diagonal matrix $D$ has determinant $\prod_{i=0}^{n-1}\binom{n-1}{i}$. Since the order-reversing matrix $R$ has may be converted to the identity matrix by $\lfloor n/2 \rfloor$ row transpositions, we have an overall sign of $(-1)^{\lfloor n/2\rfloor}$. (For comparison with Yulia's answer, this may also be written as $(-1)^{n(n-1)/2}$ since both exponents have the same parity over integers.) Multiplying these factors, we have the final result for the $k=n-1$ case as
$$\displaystyle \boxed{(-1)^{\lfloor n/2\rfloor}\prod\limits_{i=0}^{n-1}\binom{n-1}{i}\prod\limits_{1\leq i< j\leq n}(x_j-x_i)(y_j-y_i)}$$