The GCD-matrix: generalizing a result of Smith?

For Pell numbers, the answer appears to be $$\det\left[\gcd(P_i,P_j)\right]_{i,j=1}^n = \prod_{k=1}^n \sum_{d|k} \mu(\frac{k}{d})\cdot P_d,$$ i.e. the product of first $n$ terms of the Moebius transform of Pell numbers. At least, this equality holds for all $n\leq 100$.

---

UPDATE

A more general statement holds:

Theorem. For any integer $n>0$ and any variables $v_1,v_2,\dots,v_n$, $$\det [v_{\gcd(i,j)}]_{i,j=1}^n = \prod_{k=1}^n \sum_{d|k} \mu(k/d)\cdot v_d.$$

The statement for Lucas sequences directly follows from this theorem, thanks to the property $\gcd(L_i,L_j)=L_{\gcd(i,j)}$.

The theorem can be proved using the following lemma:

Lemma. For any integers $0<m<n$, $$\sum_{d|n} \mu(n/d)\cdot v_{\gcd(m,d)} = 0.$$

Indeed, in matrix $V=[v_{\gcd(i,j)}]_{i,j=1}^n$, one can replace the row $V_n$ with the linear combination $\sum_{d|n} \mu(n/d)\cdot V_d$ (which does not change the determinant, and nullifies the $n$-th row except for its rightmost element), and then use induction on $n$.

Suppose $A_n$ is the set of natural numbers that divide $L_n(s,t)$ but don't divide any $L_m(s,t)$ for $m<n$. Then corollary 1 to theorem 1 in GCD closed sets and determinants of GCD matrices by Beslin and Leigh shows that $$\det\left[\gcd(L_i(s,t),L_j(s,t))\right]_{i,j=1}^n =\prod_{i=1}^n \left(\sum_{d\in A_i}\varphi(d)\right).$$ Showing that $\sum_{d\in A_n}\varphi(d)=\sum_{d\vert n}L_d(s,t)\cdot\mu\left(\frac{k}d\right)$ is a quick application of Mobius inversion and the identity $n=\sum_{d|n} \varphi(d)$.