Identity with nested sum taken over divisors of $\gcd$'s

For computational reasons, I want to show that the following holds true:

Let $n_1,n_2,N\in \mathbb{N}$. One has $$\Large \sum_{a\mid \gcd(n_1,N)}\sum_{b\mid \gcd(n_2,\frac{n_1N}{a^2})} ab =\sum_{g\mid \gcd(n_1,n_2)}\sum_{d\mid \gcd(\frac{n_1n_2}{g^2},N)}gd. $$

It is true numerically, but I do not know how to prove it mathematically

Any comment ? Thanks.

PS: Note that all sums in the above equation runs over positive divisors only.

barto's comment to the OP's question suggest exchanging the names of $n_1$ and $N$ : the question becomes $$ \sum_{a|(n_1, N)}\sum_{b|(n_2, n_1N/a^2)}ab \overset{?}{=} \sum_{g|(n_2, N)}\sum_{d|(n_1, n_2N/g^2)}gd .$$ In other words, the left-hand side must be unchanged if we swap $n_1$ and $n_2$. In terms of the generating series, the identity you want to show is therefore equivalent to proving that for each $N\geq 1$, the function $$ (1) \qquad F(s, w) = \sum_{n_1, n_2\geq 1}n_1^{-s} n_2^{-w} \sum_{a|(n_1, N)}\sum_{b|(n_2, n_1N/a^2)} ab $$ satisfies $F(s, w) = F(w, s)$. That's by unicity of the coefficients of a Dirichlet series, granted that it has finite abscissa of absolute convergence. To check this, note that the coefficients of $F(s, w)$ are positive, and $$ F(s, w) \leq \sum_{n_1, n_2\geq 1}n_1^{-s}n_2^{-w} \sum_{a|n_1}a\sum_{b|n_2}b \leq \sum_{n_1, n_2\geq 1}d(n_1)n_1^{1-s}d(n_2)n_2^{1-w} = \zeta(s-1)^2\zeta(w-1)^2 $$ where $d(n)$ stands for the divisor counting function. If $s$ and $w$ are both $>2$, the left hand side is finite so this implies absolute convergence of $F(s, w)$, and you are guaranteed unicity of coefficients.

The standard procedure to evaluate the LHS of (1) to switch summation, and use the fact that $(a, b)|c$ is equivalent to $a|c$ and $b|c$.

Push the summations over $n_1$ and $n_2$ to the end. You have $$ F(s, w) = \sum_{\substack{a, b\geq 1\\ a|N}}\sum_{\substack{n_1, n_2\geq 1\\ a|n_1, b|n_2 \\ a^2b | n_1N}} abn_1^{-s}n_2^{-w}. $$ Rewrite $n_1 = an_1'$ and $n_2=bn_2'$, then $$ F(s, w) = \sum_{\substack{a, b\geq 1\\ a|N}}\sum_{\substack{n_1', n_2'\geq 1\\ ab | n_1'N}} a^{1-s}b^{1-w}(n_1')^{-s}(n_2')^{-w}. $$ Note that $$ ab|n_1'N \Leftrightarrow \frac{ab}{(ab, N)}\Big|n_1'\frac{N}{(ab, N)} $$ so that actually $ab/(ab, N) | n_1'$. Rewriting $n_1' = (ab/(ab, N))n_1''$, you get $$ F(s, w) = \sum_{\substack{a, b\geq 1\\ a|N}}a^{1-2s}b^{1-w-s}(ab, N)^s \sum_{\substack{n_1'', n_2'\geq 1}} (n_1'')^{-s}(n_2')^{-w} = H(s,w)\zeta(s)\zeta(w) $$ where $$ H(s, w) = \sum_{\substack{a, b\geq 1\\ a|N}}a^{1-2s}b^{1-w-s}(ab, N)^s $$ and we want to show that $H(s, w)=H(w, s)$. Since $a|N$, then $(ab, N)=a(b, N/a)$ and so $$ H(s, w) = \sum_{\substack{a, b\geq 1\\ a|N}}a^{1-s}b^{1-w-s}(b, N/a)^s $$ Let $u=(b, N/a)$, and write $b=ub'$, then $$ H(s, w) = \sum_{\substack{a, b', u\geq 1\\au|N \\(b', N/(au))=1}}u^{1-w}a^{1-s}(b')^{1-w-s} .$$ Upon writing $d=au$, this is rewritten as $$ H(s, w) = \sum_{\substack{d|N, b\geq 1 \\ (b, N/d)=1}} \Big(\sum_{d=au}a^{1-s}u^{1-w}\Big)b^{1-w-s} $$ whose symmetry is obvious, because $a$ and $u$ play the same role.

You can probably reverse-engineer the changes of variables to deduce a magic 1-1 correspondance between $(a, b)$ and $(d, g)$ in your original sum and obfuscate the rest.

$\newcommand{\lcm}[0]{\mathrm{lcm}}$Allow me to rewrite this as $$\large \sum_{a \mid \gcd(x, z)} \ \sum_{b \mid \gcd(y,\frac{x z}{a^2})} a b = \sum_{c \mid \gcd(x, y)} \ \sum_{d\ \mid \gcd(z, \frac{x y}{c^2})} c d. $$

I will provide a proof (barring mistakes) in the special case when $\gcd(x, y) = 1$. I have some hope to reduce the general case to this one.

So in this case we have $\lcm(x, y) = x y$, and thus $$ \gcd(z, x y) = \gcd(z, x) \gcd(z, y), $$ where the two factors are coprime. One sees thus that it is enough to prove that $$ \gcd(y, z) = \gcd(y, \frac{x z}{a^{2}}) $$ for each $a \mid \gcd(x, z)$.

Clearly if $t \mid \gcd(y, \dfrac{x z}{a^{2}})$, then $t \mid x z$, and as $t \mid y$, we have $\gcd(t, x) = 1$, so that $t \mid z$. It follows that $t \mid \gcd(y, z)$, or $$ \gcd(y, \frac{x z}{a^{2}}) \mid \gcd(y, z). $$

Conversely, suppose $t \mid \gcd(y, z)$. Since $a \mid x$, we will have $\gcd(t, a) = 1$. Now $$ t \mid x z = a^{2} \cdot \frac{x z}{a^{2}}, $$ and thus $t \mid \dfrac{x z}{a^{2}}$. This shows that $$ \gcd(y, z) \mid \gcd(y, \frac{x z}{a^{2}}). $$

I do not have a proof, but I will try to classify this problem in terms of linear operators. Maybe someone later can prove the claimed identity using this idea.

Let $f:\mathbb{N}\to \mathbb{R}$ be a multiplicative arithmetic function.

Define a linear operator $T(n)$ on $f$ as follows: $$ (T(n)f)(N):=\sum_{a\mid \gcd(n,N)}f(\frac{Nn}{a}). $$ It is clear that $(T(1)f)(N)=f(N)$.

For two positive integers $n_1,n_2$ one has $$ \big(T(n_2)\circ T(n_1)f\big)(N)= \sum_{a\mid \gcd(n_1,N)}\sum_{b\mid \gcd(n_2,\frac{Nn}{a^2})}f(\frac{Nn_1n_2}{ab}). $$ On the other hand, one has $$ \Big(\sum_{g\mid \gcd(n_1,n_2)}(T(\frac{n_1n_2}{g^2})f)(N)\Big)= \sum_{g\mid \gcd(n_1,n_2)}\sum_{d\mid \gcd(\frac{n_1n_2}{g^2},N)}f(\frac{n_1n_2}{g d}). $$ Comparing the last two equations with the original identity in the first post, we can now say that the proof of this identity is equivalent to proving that

$$ T(n_2)\circ T(n_1)=\sum_{g\mid \gcd(n_1,n_2)}T(\frac{n_1n_2}{g^2}). $$

Remark Some functions like $\sigma$ divisor function, and $\tau$ function satisfy similar identities.