Determinant of a gcd matrix

After testing the first cases, one is led to conjecture the formula :
$$\operatorname{det}\left(A_n\right)=\prod_{i=1}^n \varphi(i)$$
where $\varphi$ is Euler's totient function.
One first trick to make this function appear in the coefficients of the matrix is to remember that, for any integer $m\geq 1$, we have that :
$$m=\sum_{d | m} \varphi(d)$$
Applying this to the gcd, we obtain :
$$\operatorname{gcd}(i, j)=\sum_{d \mid \operatorname{gcd}(i, j)} \varphi(d)=\sum_{d \mid i \text { and }d|j} \varphi(d)$$
This can be rewritten as a sum over all $1\leq d\leq n$ by introducing the symbols $\delta_{d|i}$, which are 1 if $d|i$ and 0 if not :
$$\operatorname{gcd}(i, j)=\sum_{d=1}^{n} \delta_{d \mid i} \varphi(d) \delta_{d \mid j}$$
But this is reminiscent of the formula for the coefficient of a matrix of the form $L D U$ where $D$ is diagonal. More precisely here, we have :
$$\operatorname{gcd}(i, j)=\left(L \operatorname{diag}(\varphi(1), \ldots, \varphi(n)) L^{\top}\right)_{ij}$$
where
$L$ is the upper triangular matrix with coefficients $L_{ij}=\delta_{j|i}$. So the determinant of our matrix is given by :
$$\operatorname{det}\left(A_n\right)=\prod_{i=1}^n \varphi(i) \cdot \operatorname{det}(L)^2$$
But it is straightforward to see that $L$ is a triangular matrix with 1s on the diagonal, so its determinant is 1. We have arrived at the wanted result.

NB: The decomposition of our matrix $A_{n}$ written above is called the $LU$ decomposition, where $L$ stands for "lower triangular" and $U$ for "upper triangular". It is really useful to compute some determinants and its unicity can lead to some beautiful formulas.