The identity is clearly invariant to scaling $a$ by a constant. (One has to check the constant appears in the determinant to power $n(n+1)$). So wlog we may assume $a(1)=1$.

We can expand $\log a$ as a formal power series in $\log q$, $$ a(q)= e^{ \sum_{i} a_i (\log q)^i}$$ $$a (q^j) = e ^{ \sum_i a_i j^i (\log q)^j}$$ $$\prod_{j=1}^n a(q^j) = e^{ \sum_i a_i (\sum_{j=0}^n j^i) (\log q)^j }$$

the key thing being that $(\sum_{j=0}^n j^i)$ is a polynomial in $n$ of degree $\leq n+1$. Hence if we expand this product $f(n,q)$ out as a power series in $\log q$, the coefficient of the $i$th power of $\log q$ will be a polynomial in $n$ of degree $\leq 2i$, whose leading term is a power of $a_1$.

Hence the coefficient of $\log q^i$ in the power series of $f(x+y,q)$ is a sum of monomials in $x$ times monomials in $y$ of total degree $\leq i$.

When we take a determinant of a sum of rank one matrices, in this case a sum indexed by powers of $\log q$ and then pairs of monomials, it's the same as the sum of the determinants of all sums of $n+1$ matrices from the set. These determinants are nonzero only if the degreees of all the monomials in $x$ and monomials in $y$ appearing are distinct, which means the total degree in each variable must be at least ${n+1 \choose 2}$, for a total degree of at least $2 {n+1 \choose 2}$, and thus total degree in $\log q$ of at least ${n+1 \choose 2}$. This is only sharp if the monomials we study are the highest degree terms, which we saw come only from $a_1$.

So the answer depends only on $a_1$, and the calculations you include in your answer finish the job.