It is a basic result that $$ \log n \sim H_n \equiv \sum_{k = 1}^n \frac{1}{n}. $$ We also have $$ \log(mn) = \log m + \log n. $$ Based on this, it seems like we should be able to prove something along the lines of: for every $\epsilon > 0$ there exists $N$ such that when $n,m > N$ $$ \mid H_{mn} - (H_{m} + H_n) \mid < \epsilon. $$ Is something like this true? If so, how does one go about proving it?
Edit
Using the estimate $$ H_n = \log n + \gamma + O\left(\frac{1}{n}\right) $$ We find \begin{align} H_{mn} - (H_m + H_n) &= \left( \log (mn) + \gamma + O\left(\frac{1}{mn}\right) \right) - \left(\log m + \gamma + O\left(\frac{1}{m}\right) + \log n + \gamma + O\left(\frac{1}{n}\right)\right) \\&= -\gamma + O\left(\frac{1}{m}\right) + O\left(\frac{1}{n}\right) \end{align} Or, less formally, for large $m$ and large $n$ $$ H_{mn} \approx H_m + H_n - \gamma $$ Interestingly, this gives us yet another way to calculate $\gamma$.
In my opinion, you missed an important term in the asymptotics of $H_n$
$$H_n=\gamma +\log \left({n}\right)+\frac{1}{2 n}+O\left(\frac{1}{n^2}\right)$$ which would make $$H_{mn}-H_m-H_n\sim-\gamma -\frac{1}{2 m}-\frac{1}{2 n}+\frac{1}{2 mn}+\cdots$$