Where does the identity $\tau(nm) = \sum_{d \mid (m, n)} \mu(d) \tau(n/d) \tau(m/d)$ come from? I was reading an article which states (without reference) a recurrence relation for a character sum associated to a quadratic extension of number fields. It is a bit complicated to state, but the analogue for the trivial extension $\mathbb Q / \mathbb Q$ and the trivial character would be that
$$\tau(nm) = \sum_{d \mid (m, n)} \mu(d) \tau(n/d) \tau(m/d)$$
where $\tau$ denotes the divisor counting function. This seems to be true. I quickly checked this when $mn$ has $1$ or $2$ prime divisors; it then reduces to the identities
$$a+b+1 = (a+1)(b+1) - ab$$
resp.
$$\begin{align*}(a+b+1)(c+d+1)
& = (a+1)(b+1)(c+1)(d+1) \\
& - ab(c+1)(d+1) \\
& - (a+1)(b+1)cd \\
& + abcd
\end{align*}$$
Where does this identity (if true) come from? Can it be understood by comparing the coefficients of two Dirichlet series? Modular forms? It is reminiscent of the recurrence relation for Hecke operators.
 A: This formula is a theorem in the context of arithmetical functions. The divisor function $\tau$ is a multiplicative function which can be represented as convolution of completely multiplicative functions, namely as Dirichlet convolution of the zeta-function 
\begin{align*}
\zeta(n)=1\qquad n\geq 1
\end{align*}
We have for $n\geq 1$
\begin{align*}
\left(\zeta\star\zeta\right)(n)=\sum_{d|n}\zeta(d)\zeta(n/d)=\sum_{d|n}1=\tau(n)
\end{align*}

We find e.g. in Introduction to Arithmetical Functions by P.J. McCarthy:
Theorem 1.12: If $f$ is a multiplicative function then the following statements are equivalent:
  
  
*
  
*(1) $f$ is a convolution of two completely multiplicative functions
  
*(2) There is a multiplicative function $F$ such that for all $m$ and $n$,
\begin{align*}
f(mn)=\sum_{d|(m,n)}f(m/d)f(n/d)F(d)\tag{1}
\end{align*}
  
  
*
  
*(3) There is a completely multiplicative function $B$ such that for all $m$ and $n$
\begin{align*}
f(m)f(n)=\sum_{d|(m,n)}f(mn/d^2)B(d)
\end{align*}
  
  
*
  
*(4) For all primes $p$ and all $\alpha \geq 1$,
  
  
  \begin{align*}
f\left(p^{\alpha+1}\right)=f(p)f\left(p^{\alpha}\right)+f\left(p^{\alpha-1}\right)\left[f\left(p^2\right)-f\left(p\right)^2\right]
\end{align*}

In the current situtation we have $F(d)=\mu(d)$ and get from (1)
\begin{align*}
\tau(mn)=\sum_{d|(m,n)}\tau(m/d)\tau(n/d)\mu(d)
\end{align*}
Multiplicative functions which have the property stated in Theorem 1.12 are called specially multiplicative functions.
A: For the proof of your identity let
$$n =N \prod_j p_j^{a_j+1}, \qquad m =M \prod_j p_j^{b_j+1}, \qquad \gcd(N,M)=\gcd(NM,\prod_j p_j)=1$$
Then $$\tau(nm) = \tau(N)\tau(M) \prod_j \tau(p_j^{a_j+b_j+2})=  \tau(N)\tau(M) \prod_j  ( \tau(p_j^{a_j+1})\tau(p_j^{b_j+1})- \tau(p_j^{a_j})\tau(p_j^{b_j}))$$
$$=\tau(N)\tau(M) \prod_{p_j^r \| (n,m)} \sum_{s=0}^r \mu(p_j^s)\tau(p_j^{a_j+1-s})\tau(p_j^{b_j+1-s})$$
$$ = \tau(N)\tau(M)\sum_{d |(n,m)} \mu(d) \tau(\frac{n}{Nd})\tau(\frac{m}{Md})$$ $$=\sum_{d |(n,m)} \mu(d) \tau(n/d)\tau(m/d)$$
