Hello, I really think it helps to think about the entries of this matrix does $q^{d(\sigma,\tau)}$ where d denotes the distance between two permutations on the Cayley graph of the symmetric group generated by the set of all transpositions. In this light you can compare the determinant you're interested in with the one computed by Don Zagier in his article "realizability of a model in infinite statistics." In that paper, Zagier computes the analogous determinant corresponding to the Cayley graph of the symmetric group generated by adjacent transpositions (i.e. the Coxeter generators). This is a very different geometry on permutations but Zagier finds a similar factorization showing that the determinant has roots on the unit circle. Zagier's result is based on a clever factorization in the group algebra, just like the factorization into JM elements discussed above.
Another point of view on Zagier's result, relating it to hyperplane arrangements, is explained in Stanley and Hanlon's article "a q-deformation of a trivial symmetric group action." Also, if you want to see similar results in the context of the Brauer algebra, you should look up Paul Zinn-Justin's paper "Jucys-Murphy elements and Weingarten matrices." I saw that a user above posted a link to my Phd thesis - I was young and foolish then and more polished versions of those results have since appeared in a paper written jointly with Sho MatsumoMatsumoto called "Jucys-Murphy elements and unitary matrix integrals."