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)$?

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}$.
