A derivation of a Riemann zeta function identity

Yesterday, I saw the following Riemann zeta function identity:

\displaystyle\sum_{n=1}^{\infty} \frac{\sigma_a(n)\sigma_b(n)}{n^s} = \frac{\zeta(s)\zeta(s-a)\zeta(s-b)\zeta(s-a-b)}{\zeta(2s-a-b)}.

I took some time to try to derive it myself and to my great pleasure, I succeeded.

Eventually, I realized that it suffices to show that

\{(dd_a, dd_b, d^2 n) : d_a | n, d_b | n : d, d_a, d_b, n \in \mathbb{Z}\}

and

\{(dd_a, dd_b, n) : dd_a d_b | n : d, d_a, d_b, n \in \mathbb{Z}\}

are equal as multisets. As sets, they are both representations of the set of 3-tuples of positive integers such that the third is a multiple of the least common multiple of the first two. In the latter one, the frequency of (a,b,c) is the number of d that divides both a and b such that ab | cd. In the other one, if we write (a,b,c) as (d_1 d_2 a', d_1 d_2 b', c) where \mathrm{gcd}(a', b') = 1, the ab | cd condition equates to d_1^2 d_2 a'b' | c, which corresponds to the number of d_1 dividing a and b and such that d_1^2 | c and with that, d_2a', d_2b' both dividing d_1^2 / c, which is the frequency of (a,b,c) via the former representation.

The coefficients \{a_n\} of the Dirichlet series of the LHS of that identity can be decomposed as follows:

a_n = \displaystyle\sum_{d^2 | n, d_a | \frac{n}{d^2}, d_b | \frac{n}{d^2}} (dd_a)^a (dd_b)^b.

The coefficients \{b_n\} of the Dirichlet series of the RHS of that identity are

b_n = \displaystyle\sum_{dd_a d_b | n} (dd_a)^a (dd_b)^b.

Observe how both are equivalent in that via the multiset equivalence proved above, n determines the same multiset of (dd_a, dd_b) for both and across that, the values of the same function (dd_a)^a (dd_b)^b are summed. Hence the two series are equal.