As this hasn't attracted an answer yet, I will make my comments.
IMO the authors simply claimed this, and if I am guessing their intended method correctly, it doesn't actually work as easily as presumed.
Let's ignore the $m_1$, drop the $Y$, and use $\log (X/m)\le\log X$ in the first sum.
Then we want to show (here $n$ is $m_2$)
$$\sum_{n<X} {|a(n^2)|\over n^2}\ll \log X.$$
My guess is that the authors basically assumed that this was similar to $\sum_n |a(n)|^2/n^2$, which is indeed the classical Rankin convolution. The routine argument says (e.g.)
$$\sum_{n\le X}{|a(n)|^2\over n^2}
\ll\sum_n{|a(n)|^2\over n^2}e^{-n/X}
=\int X^s\Gamma(s)L(f\otimes f,s+2){ds\over 2\pi i}$$
and the self-convolution $L$-function has a simple pole at $s=2$, thus a double pole for the integrand at $s=0$, giving $\ll\log X$. (You can include the $m_1$ as a finite Euler product that is nice somewhat to the left of $\sigma=0$.) A similar argument presumably works with $\tau(n)$ included (though not $\tau(n^2)$ if I am correct), as that just appends another $\log X$.
But the situation with $|a(n^2)|$ is different. Again I can't know their thoughts, and the method I propose is (much) too complicated to go without comment. Here goes: note $a(n^2)$ gives (up to bad Euler factors and a tame $\zeta$-factor) the symmetric square coefficients (see middle of page 450). Then the $L^1$-norm of these symmetric square coefficients can be bounded, similar to the paper's Lemma 17 (Rankin's bound for $\sum_n |a(n)|$), with a savings of a power of log. (I'm not actually sure how much one needs holomorphy of the self-convolution of the symmetric square, or indeed whether the details of this were known at the time.) Then by partial summation one gets the desired bound. The second term is similar (again one can use smoothing as above), but again I think $\tau(m_2^2)$ is going to lose an extra logarithm.