Computing a determinant involving roots of unity Let $d \geq 2$ be an integer and $\xi=\exp(\frac{2\pi i}{d})$. I am trying to compute the determinant of the matrix
$$
(\xi^{ij}-1)_{1 \leq i, j \leq d-1}.
$$ Let me call it $\Delta(d)$. For small values of $d$ I get:
    $\Delta(2)=-2$
    $\Delta(3)=-3\sqrt{3}i$
    $\Delta(4)=-16i$
But I can't seem to find a general formula. How can I do this?
 A: (Too long for a comment) Replacing $\xi$ with a variable $x$, Mathematica gives the following:
\begin{align*} 
\Delta(2)&=-1+x\\
\Delta(3)&=(-1+x)^3 x (1+x)\\
\Delta(4)&=(-1+x)^6 x^4 (1+x)^2 \left(1+x+x^2\right)\\
\Delta(5)&=(-1+x)^{10} x^{10} (1+x)^4 \left(1+x^2\right) \left(1+x+x^2\right)^2\\
\Delta(6)&=(-1+x)^{15} x^{20} (1+x)^6 \left(1+x^2\right)^2 \left(1+x+x^2\right)^3 \left(1+x+x^2+x^3+x^4\right)
\end{align*}
 etc...
There is a nice pattern: for $d\ge3$, $\Delta(d)$ contains as factors exactly the cyclotomic polynomials of degree 1 to $d-1$, and $x$, raised to some powers.
Edit: I checked up to $d=12$ and, as Benjamin remarks, the exponents seems to be:


*

*triangular numbers for the factor $x-1$ (see https://en.wikipedia.org/wiki/Triangular_number)

*tetrahedral numbers for the factor $x$ (see https://en.wikipedia.org/wiki/Tetrahedral_number)


I searched in the OEIS for other exponent sequences, but it returns a lot of possibilities. For example, the exponents over $1+x$ go like this: 1,2,4,6,9,12,16,20,25,30,..., and the OEIS gives 16 results: https://oeis.org/search?q=1%2c2%2c4%2c6%2c9%2c12%2c16%2c20%2c25%2c30
Another edit: Problem solved, see at the end of T. Amdeberhan's answer.
A: Let $\Phi=(\xi^{ij})_{0\leq i,j \leq d-1}$. Then $\det\Phi=\prod_{0\leq i<j < d} (\xi^i-\xi^j)=(-1)^{\binom{d}2}\prod_{0\leq i<j < d} (\xi^j-\xi^i)$. So,
\begin{align} \prod_{0\leq i<j < d} (\xi^j-\xi^i)^2
&=\det\begin{bmatrix}1&1 &\dots &1 \\ x_0&x_1&\dots&x_{d-1} \\
\vdots&\vdots&\dots&\vdots \\ x_0^{d-1}&x_1^{d-1}&\dots&x_{d-1}^{d-1} 
\end{bmatrix}
\det\begin{bmatrix}1&x_0 &\dots &x_0^{d-1} \\ 1&x_1&\dots&x_1^{d-1} \\
\vdots&\vdots&\dots&\vdots \\ 1&x_{d-1}&\dots&x_{d-1}^{d-1} \end{bmatrix} \\
&=\det\begin{bmatrix}t_0&t_1 &\dots &t_{d-1} \\ t_1&t_2&\dots&t_d \\
\vdots&\vdots&\dots&\vdots \\ t_{d-1}&t_d&\dots&t_{2(d-1)}\end{bmatrix};
\end{align}
where $x_k=\xi^k$ and $t_k=x_0^k+x_1^k+\cdots+x_{d-1}^k$. Notice that $t_k=0$ unless $k$ is a multiple $d$, in which case it is equal to $d$. Therefore,
$$\prod_{0\leq i<j < d} (\xi^j-\xi^i)^2
=d^d\cdot\det
\begin{bmatrix}1&0&0&\dots&0&0 \\ 0&0&0&\dots&0&1 \\ 
0&0&0&\dots&1&0 \\
\vdots&\vdots&\vdots&\dots&\vdots&\vdots \\
0&1&0&\dots&0&0\end{bmatrix}=d^d(-1)^{\binom{d-1}2}.$$
We gather that
$$\det\Phi=(-1)^{\binom{d}2}\cdot d^{\frac{d}2}\cdot i^{\binom{d-1}2}
=d^{\frac{d}2}\cdot i^{\frac{(3d-2)(d-1)}2}.$$
A couple of noteworthy facts:  if $V=\prod_{0\leq i<j < d} (\xi^j-\xi^i)$ denote the Vandermonde determinant, then 
$V^2$ is the discriminant of the polynomial $f(z)=z^d-1$;
$V^2$ is a symmetric polynomial and when expanded in elementary polynomials basis then only $e_d(\xi_0,\dots,\xi_{d-1})$ survives.
By the way, in general, computing the coefficients of even-powers of the Vandermonde's determinant $V(y_1,\dots,y_n)$ symmetric polynomials basis is a difficult and big sport in mathematics and physics.
UPDATE. This is in response to EFinat-S's quest on replacing $\xi$ by an indeterminate $x$. The determinant is just
$$\Delta(d)=\det\left[x^{ij}-1\right]_{i,j}^{1,d}
=\prod_{0\leq i<j\leq d}\left(x^j-x^i\right).$$ 
POSTSCRIPT. We can generalize the determinant evaluation as follows. Start with the Vandermonde
$$\prod_{i<j}^{\pmb{0,d-1}}(y_j-y_i)=\det\begin{bmatrix}1&y_0 &\dots &y_0^{d-1} \\ 1&y_1&\dots&y_1^{d-1} \\
\vdots&\vdots&\dots&\vdots \\ 1&y_{d-1}&\dots&y_{d-1}^{d-1} \end{bmatrix},$$
then choose $y_0=1$ and subtract the first column from the rest to get
$$\det\begin{bmatrix}1&y_0-1 &\dots &y_0^{d-1}-1 \\ 1&y_1-1&\dots&y_1^{d-1}-1 \\
\vdots&\vdots&\dots&\vdots \\ 1&y_{d-1}-1&\dots&y_{d-1}^{d-1}-1 \end{bmatrix}
=\det\begin{bmatrix}1&y_1-1 &\dots &y_1^{d-1}-1 \\ 1&y_2-1&\dots&y_2^{d-1}-1 \\
\vdots&\vdots&\dots&\vdots \\ 1&y_{d-1}-1&\dots&y_{d-1}^{d-1}-1 \end{bmatrix}.$$
Consequently,
$$\det\begin{bmatrix}1&y_1-1 &\dots &y_1^{d-1}-1 \\ 1&y_2-1&\dots&y_2^{d-1}-1 \\
\vdots&\vdots&\dots&\vdots \\ 1&y_{d-1}-1&\dots&y_{d-1}^{d-1}-1 \end{bmatrix}
=\prod_{k=1}^{d-1}(y_k-1)\cdot\prod_{i<j}^{\pmb{1,d-1}}(y_j-y_i).$$ 
The latter determinant is in the form of your question, with $y_j=\xi^j$. 
A: Using the earlier responses and comments, I confirm the formula suggested by Neil Strickland:
$$\Delta(d)=d^{d/2}i^{m(d)}\qquad\text{with}\qquad m(d): = 1 + d(7-d)/2\in\mathbb{Z}.$$
Consider the $d\times d$ Vandermonde matrix
$$\Phi(d):=(\xi^{ij})_{0\leq i,j \leq d-1}.$$
Subtracting the first column from each other column, we get a matrix with first row equal to $(1,0,\dots,0)$ and lower right $(d-1)\times(d-1)$ block equal to the OP's matrix. Therefore,
$$\Delta(d)=\det\Phi(d).$$
It is straightforward to check that $\Phi(d)^\ast\cdot\Phi(d)$ equals $d$ times the identity matrix, therefore
$$ |\det\Phi(d)|^2=d^d.$$
In other words, $|\det\Phi(d)|=d^{d/2}$, and we are left with determining
$$\frac{\det\Phi(d)}{|\det\Phi(d)|}=\prod_{0\leq i<j\leq d-1}\frac{\xi^j-\xi^i}{|\xi^j-\xi^i|}.$$
Let me use the notation $e(t):=e^{2\pi it}$, familiar from analytic number theory. Then we see that
$$\xi^j-\xi^i=e\left(\frac{j}{d}\right)-e\left(\frac{i}{d}\right)
=e\left(\frac{i+j}{2d}\right)\left(e\left(\frac{j-i}{2d}\right)-e\left(\frac{i-j}{2d}\right)\right).$$
On the right hand side, $0<\frac{j-i}{2d}<\frac{1}{2}$, hence $e\left(\frac{j-i}{2d}\right)$ lies in the upper half-plane. As a result,
$$\frac{\xi^j-\xi^i}{|\xi^j-\xi^i|}=e\left(\frac{i+j}{2d}\right)i.$$
We need to calculate the product of the right hand side over the $\binom{d}{2}$ pairs $0\leq i<j\leq d-1$. By symmetry (or by brute-force calculation), the average of $i+j$ equals $d-1$, whence
$$\prod_{0\leq i<j\leq d-1}\frac{\xi^j-\xi^i}{|\xi^j-\xi^i|}=\left(e\left(\frac{d-1}{2d}\right)i\right)^{\binom{d}{2}}=e\left(\left(\frac{d-1}{2d}+\frac{1}{4}\right)\binom{d}{2}\right).$$
We calculate
$$\left(\frac{d-1}{2d}+\frac{1}{4}\right)\binom{d}{2}=\frac{(3d-2)(d-1)}{8},$$
therefore in the end
$$\Delta(d)=d^{d/2}i^{n(d)}\qquad\text{with}\qquad n(d):=(3d-2)(d-1)/2\in\mathbb{Z}.$$
This agrees with Neil Strickland's formula, upon noting that $m(d)\equiv n(d)\pmod{4}$, i.e.,
$$2+d(7-d)\equiv (3d-2)(d-1)\pmod{8}.$$
Added 1. As Alexey Ustinov remarked, $\Phi(d)$ is known as a Schur matrix. As Carlitz wrote in his 1959 paper, "this matrix is familiar in connection with Schur's derivation of the value of Gauss's sum". In fact, on page 295 of this paper, Carlitz uses the known value of Gauss's sum to find the eigenvalues of this matrix (which are all of the form $\pm\sqrt{d}$ and $\pm i\sqrt{d}$, hence one only needs to find the 4 multiplicities). This can be regarded as a refinement and an alternate proof of the above result, since the product of the eigenvalues is the determinant. 
Added 2. Carlitz referred to Landau's Vorlesungen, whose relevant part appeared in English as Landau: Elementary Number Theory (Chelsea, 1958). So I looked up this translation, and to my surprise on pages 211-212 I found essentially the same calculation as above. In fact all this is in Schur's 1921 paper, thanks to Alexey Ustinov for locating it for me (alternative link here). As Landau explains (following Schur), the determinant calculation leads to an evaluation of Gauss's sum, at least for odd $d$. The point is that one can figure out the 4 eigenvalue multiplicities from the determinant, and hence one obtains the formula for the trace of $\Phi(d)$ as well. However, this trace is nothing but Gauss's sum!
Added 3. For a more recent treatment of Schur's evaluation of the Gauss sum, see Section 6.3 in Rose: A course in number theory (2nd ed., Oxford University Press, 1994).
A: Numerical experiment (as far as $d=50$) makes it clear that $\Delta(d)=d^{d/2}i^{m_d}$ where 
$$ m_d = 1 + d(7-d)/2 \in\mathbb{Z}. $$
I haven't tried to find a proof, however.
A: Partial answer, to expand on my suggestion as requested by OP: let $\Phi=(\xi^{ij})_{0\leq i,j \leq d-1}$ be the FFT matrix; in particular, its first row and column contain all ones. Let
$$
L = \begin{bmatrix}1 \\ -1 & 1\\ -1 && 1 \\  \vdots & & & \ddots \\ -1 & &&&1\end{bmatrix}
$$
be the elementary matrix that does one step of Gaussian elimination. One can check that
$$
L\Phi = \begin{bmatrix}1 & 1 & \dots & 1\\ 0\\\vdots &  & M\\0\end{bmatrix},
$$
where $M$ is the matrix whose determinant you wish to compute. By comparing determinants, $\det M = \det \Phi$. The Vandermonde determinant formula gives
$$
\det \Phi = \prod_{0\leq i<j < d} (\xi^i-\xi^j).
$$
So the problem can be reduced to finding a closed formula for this product. I don't have an immediate solution but it looks like a manageable problem.
