Computing a determinant involving roots of unity

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).

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.

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.