Wronskian of functions $\sin(nx), n=1,2,...,k$. Is it true that the Wronskian of the functions $\sin(nx), n=1,...,k$ is equal to $c(\sin(x))^p$ where $c$ is a constant and $p=1+2+...+k=k(k+1)/2$?
That is true for $k=1,2,3,4,5$. If it is true, how to find the constant $c=c(k)$?
 A: Let $W(f_1, \ldots, f_n)$ denote the Wronskian determinant of the functions $f_1, \ldots, f_n$. We can show that

$$
W (\sin(x), \sin(2x), \ldots, \sin(nx)) = 
1!2! \cdots (n-1)! (-2)^{n(n-1)/2} \sin(x)^{n(n+1)/2} \, .
$$

For example,
$$
\begin{align}
 W (\sin(x), \sin(2x)) &= -2 (\sin(x))^3 \, .\\
 W (\sin(x), \sin(2x), \sin(3x)) &= -16 (\sin(x))^{6} \, ,\\
 W (\sin(x), \ldots, \sin(4x)) &= 768 (\sin(x))^{10} \, .
\end{align}  
$$
The proof uses that
$$
 \sin(k x) =  U_{k-1}(\cos (x)) \sin(x)
$$
where $U_k$ are the Chebyshev polynomials of the second kind, and two identities for Wronskians: A “product rule”
$$
W(h f_1, \ldots, h f_n) = h^n \cdot W(f_1, \ldots, f_n)
$$
which is a consequence of the Leibniz rule for the $n$th derivative of a product (see also Why does the Wronskian satisfy $W(yy_1,\ldots,yy_n)=y^n W(y_1,\ldots,y_n)$?), and a “chain rule”
$$
 W(f_1 \circ g, \ldots, f_n \circ g)(x) = W(f_1, \ldots f_n)(g(x)) \cdot (g'(x))^{n(n-1)/2} \, .
$$
which is a consequence of Faà di Bruno's formula for the $n$th derivative of a composite function (compare also About a chain rule for Wronskians).
Now we can argue as follows:
$$
\begin{align}
&W (\sin(x), \sin(2x), \ldots, \sin(nx)) \\
&\quad = W(U_0(\cos(x))\sin(x), U_1(\cos (x)) \sin(x), \ldots,  U_{n-1}(\cos (x)) \sin(x)) \\
&\quad = (\sin(x))^n W(U_0(\cos (x)), U_1(\cos (x)), \ldots,  U_{n-1}(\cos (x))) \, \\
&\quad = (\sin(x))^n W(U_0(t), U_1(t), \ldots,  U_{n-1}(t)) |_{t=\cos(x)} (-\sin(x))^{n(n-1)/2} \, .
\end{align}
$$
Each $U_k$ is a polynomial of degree $k$ with the leading coefficient $2^k$, so that $W(U_0, U_1, \ldots, U_{n-1})$ is the determinant of a triangular matrix with the entries $U_k^{(k)}(t) = k!2^k$, $k=0, \ldots, n-1$ on the diagonal. It follows that
$$
 W (\sin(x), \sin(2x), \ldots, \sin(nx)) = (\sin(x))^n \cdot (-\sin(x))^{n(n-1)/2}
\cdot \prod_{k=0}^{n-1} k!2^k
$$
and that is the claimed formula.
A: Consider the Wronskian of $e^{inx}$.
\begin{align}W(e^{inx})&=\det\begin{pmatrix}e^{ix}&e^{2ix}&\cdots&e^{inx}\\
ie^{ix}&2ie^{2ix}&\cdots&ine^{inx}\\
\vdots\\
i^{n-1}e^{ix}&(2i)^{n-1}e^{2ix}&\cdots&(in)^{n-1}e^{inx}\end{pmatrix}\\
&=e^{ix(1+2+\cdots+n)}(i^{1+2+\cdots+n-1})\det\begin{pmatrix}1&1&\cdots&1\\1&2&\cdots&n\\\vdots\\1&2^{n-1}&\cdots&n^{n-1} \end{pmatrix}\\
&=e^{ixn(n+1)/2}i^{n(n-1)/2} d(n)\end{align}
The constants $c(n)$ turn out to be $2^{n(n-1)/2}d(n)$, so there is definitely a connection.
Since $\sin kx=(e^{ikx}-e^{-ikx})/2i$, we can use the linearity properties of the determinant as follows:
\begin{align}
W(sin(nx))&=\frac{1}{(2i)^{n}}\det\begin{pmatrix}e^{ix}-e^{-ix}&e^{2ix}-e^{-2ix}&\cdots\\
\vdots\\
e^{ix}\pm e^{-ix}&2^{n-1}(e^{2ix}\pm e^{-2ix})&\cdots\end{pmatrix}\\
&=\frac{1}{(2i)^{n}}\sum_{\sigma\in 2^n} W(\sigma)\\
&=2^{-n}i^{n(n-1)/2-n}\sum_\sigma d(\sigma)e^{ixf(\sigma)}
\end{align} where $W(\sigma)$ is the Wronskian of $\pm e^{\pm ix},\pm e^{\pm2ix},\ldots,\pm e^{\pm inx}$, with $\sigma=(\pm1,\ldots,\pm1)$ is a choice of signs for the exponentials and $d(\sigma)$ is the determinant of the corresponding matrix. For example $$W(1,-1,1)=\det\begin{pmatrix}e^{ix}&-e^{-2ix}&e^{3ix}\\ie^{ix}&2ie^{-2ix}&3ie^{3ix}\\i^2e^{ix}&-(2i)^2e^{-2ix}&(3i)^2e^{3ix}\end{pmatrix}=e^{ix(1-2+3)}i^{1+2}\det\begin{pmatrix}1&-1&1\\1&2&3\\1&-2^2&3^2\end{pmatrix}$$
Each choice $\sigma$ has a complementary choice $\sigma'=-\sigma$, so that $f(\sigma')=-f(\sigma)$. Moreover $f(\sigma)$ are all even or all odd, because a single change in sign leads to a difference of $2$ for $f$.
Also, $d(\sigma')=\pm d(\sigma)$ since the corresponding matrices have rows that are $\pm1$ of each other. So $d(\sigma')=d(\sigma)$ when $n=4k$ or $4k-1$; $d(\sigma')=-d(\sigma)$ when $n=4k+1$ or $n=4k+2$.
Thus $$W(\sin nx)=2^{-n}i^{(n-1)(n-2)/2-n}\sum_{\sigma\in2^{n-1}} d(\sigma)(e^{ixf(\sigma)}\pm e^{-ixf(\sigma)})$$
So... one more step...

As it turns out, $d(\sigma)=1.2!\cdots n!\binom{n(n+1)/2}{2k}$ and $\sum_\sigma d(\sigma)e^{ixf(\sigma)}=\prod_{k=1}^nk!(e^{ix}-e^{-ix})^{n(n+1)/2}$.
