\section{Introduction}
Moments of $L$-functions in families encode deep properties about the family and have a wide range of applications. Estimating moments of families of $L$-functions, especially their asymptotic formulae, is regarded as a central problem in number theory. There is a very general conjecture on asymptotic formulae for integral moments of $L$-functions by Conrey, et.al., \cite{CFK+05}. However, we can hardly obtain an asymptotic formula for a moment higher than fourth by now.

Moments of the Riemann zeta-function were first concerned and most extensively studied by researchers. An asymptotic formula for the fourth moment of the Riemann zeta-function may go back to Ingham \cite{Ing26}, who proved
\[
\frac1T\int_0^T|\zeta(\tfrac12+it)|^4dt=a_4(\log T)^4+O\left((\log T)^3\right)
\]
with $a_4=(2\pi^2)^{-1}$. Ingham's result only gave the leading term of the asymptotic formula.
Heath-Brown \cite{HB79} improved this result by obtaining all lower-order main terms that
\begin{align}\label{eqHeath-Brown}
\frac1T\int_0^T|\zeta(\tfrac12+it)|^4dt=P_4(\log T)+O\left(T^{-\frac18+\ve}\right),
\end{align}
where $P_4(x)$ is a computable absolute polynomial of degree $4$ (see also \cite[formula (5.1.4)]{CFK+05}). A power saving in the error term was first obtained in Heath-Brown's result. Distinguishing all lower-order main terms with a power saving in the error is a significant challenge and requires a difficult analysis on off-diagonal terms.

The analysis on off-diagonal terms is related to the divisor problem
\begin{align}\label{eqdivisor1}
\sum_{n\le x}d(n)d(n+r)
\end{align}
uniformly for $r$ as large as possible.
A divisor problem of this kind has been widely studied, and one can also refer to \cite{DFI94} for example. The treatment of the divisor problem is based on estimations on Kloosterman sums. For Heath-Brown's result \eqref{eqHeath-Brown}, Weil bound for Kloosterman sums works well, but to obtain a further power saving in the error, we require the spectral theory of automorphic forms, in particular, the Kuznetsov formula to treat sums of Kloosterman sums.  Then Zavorotnyi \cite{Zav89} improved the result to
\begin{align}\label{eqZav}
\frac1T\int_0^T|\zeta(\tfrac12+it)|^4dt=P_4(\log T)+O\left(T^{-\frac13+\ve}\right).
\end{align}
furthermore, a beautiful explicit formula for a smoothed version of the fourth moment of the Riemann zeta-function has been obtained by Monohashi
\cite{Mot93} (see also \cite[Theorem 4.2]{Mot97}), which expresses the moment in terms of the central values of certain automorphic $L$-functions. Based on this explicit formula, Ivi\'{c} and Motohashi \cite{IM95} have replaced the factor $T^\ve$ in the error of \eqref{eqZav} with a power of $\log T$, which is the current best-known result. A generalization of Motohashi's formula to the fourth moment of Dirichlet $L$-functions modulo $q$ weighted by a non-archimedean test function has also been obtained by Blomer, et.al., \cite{BHK+20}, which proceeds differently, as well as some interesting applications.

The fourth moment of the family of Dirichlet $L$-functions averaged over all primitive characters modulo $q$ turns out to be much hard to study. Although it has some analogous places to the Riemann zeta-function in $t$-aspect, as we will see, there are significant differences which cause the fourth moment of Dirichlet $L$-functions to be much more difficult to treat. It has ever be conjectured by \cite{Ram79} that, for prime $q$,
\[
q^{-1}\sum_{\chi(\bmod q)}|L(\tfrac12,\chi)|^4\sim (2\pi^2)^{-1}(\log q)^4.
\]
For $q\nequiv 2 \pmod4$, Heath-Brown \cite{HB81} proved
\[
\frac1{\varphi^*(q)}\mathop{\sum\nolimits^*}_{\chi(\bmod q)}|L(\tfrac12,\chi)|^4=\frac1{2\pi^2}\prod_{p\mid q}\frac{(1-p^{-1})^3}{(1+p^{-1})}(\log q)^4+O\left(2^{\omega(q)}\frac{q}{\varphi^*(q)}(\log q)^3\right),
\]
where the sum is over all primitive characters modulo $q$, $\varphi^*(q)$ denotes the number of such primitive characters, and $\omega(q)$ denotes the number of distinct prime factors of $q$. Note that the condition $q\nequiv2 \pmod 4$ is reasonable since there are no primitive characters modulo $q$ if $q\equiv2 \pmod 4$. This gives an asymptotic formula if $q$ is not divided by too many primes. Soundararajan \cite{Sou07} filled this gap by improving the error, so that the leading term of the asymptotic formula has been obtained for all $q$.

However, due to the conjecture for integral moments by \cite{CFK+05}, it is believed that
\[
\frac1{\varphi^*(q)}\mathop{\sum\nolimits^*}_{\chi(\bmod q)}|L(\tfrac12,\chi)|^4=\prod_{p\mid q}\frac{(1-p^{-1})^3}{(1+p^{-1})}P_4(\log q)+O\left(q^{-\frac12+\ve}\right)
\]
with $P_4(x)$ being a computable absolute polynomial of degree $4$. When $q$ is a prime, Young \cite{You11} pushed the result much close to the conjecture, more precisely, he proved
\[
\frac1{\varphi^*(q)}\mathop{\sum\nolimits^*}_{\chi(\bmod q)}|L(\tfrac12,\chi)|^4=P_4(\log q)+O\left(q^{-\frac1{80}+\frac{\theta}{40}+\ve}\right),
\]
where $\theta=7/64$ is the current best-known bound on the size of the Hecke eigenvalue of a Maass form.

Young's work also motivates the study of another difficult and interesting problem, to obtain an asymptotic formula with a power saving for the second moment
\begin{align}\label{eqMLF}
\mathop{\sum\nolimits^*}_{\chi(\bmod q)}L(\tfrac12, f\varotimes \chi)\ol{L(\tfrac12, f\varotimes \chi)}
\end{align}
of any fixed cusp forms $f$ twisted by all primitive characters modulo $q$.
Blomer and Mili\'cevi\'c \cite{BM15} solved this problem for all sufficiently factorable $q$ including $99.9\%$ of all admissible moduli. Exceptions are for those $q$ that are essentially a prime or the product of two primes of almost equal size. As observed by Fouvry, Kowalski, and Michel, the method in \cite{BM15} can be applied to improve Young's result, especially to remove the dependence of Ramanujan-Petersson conjecture by average results on Hecke eigenvalues. A remaining case for \eqref{eqMLF} with $q$ prime was then solved by Blomer, et.al., \cite{BFK+17a}, where estimates on bilinear forms in Kloosterman sums were introduced in solving the problem. Along this way, Young's exponent $-1/80+\theta/40$ was then improved to $-1/20$ by Blomer, et.al., \cite{BFK+17a, BFK+17b}. This improvement is based on two ingredients applied, one of which is some powerful results concerning bilinear forms in Kloosterman sums \cite{FKM14, KMS17, SZ16}, and the other is an average result concerning Hecke eigenvalues to remove the dependence of the Ramanujan-Petersson conjecture.

However, there are essential differences between treatments for Dirichlet $L$-functions and \eqref{eqMLF}, e.g., main terms for \eqref{eqMLF} are only contributed by diagonal terms, while an asymptotic formula of the fourth moment of Dirichlet $L$-function relies on a difficult analysis of off-diagonal terms to deduce all lower-order main terms. For the fourth moment of Dirichlet $L$-functions, it has been more than ten years since Young announced his classical result for prime moduli, while the case for general moduli is still open.
This paper is devoted to this problem, and we specify our main result in the following theorem.
\begin{theorem}\label{thmmain}
For any $q\nequiv 2 \pmod 4$, we have
\begin{align}\label{eqmain}
\frac1{\varphi^*(q)}\mathop{\sum\nolimits^*}_{\chi(\bmod q)}|L(\tfrac12,\chi)|^4=\prod_{p\mid q}\frac{(1-p^{-1})^3}{(1+p^{-1})}P_4(\log q)+O\left(q^{-\frac1{80}+\frac{\theta}{40}+\ve}\right)
\end{align}
for a certain computable absolute polynomial $P_4(x)$.
\end{theorem}

Note that Theorem \ref{thmmain} gives an asymptotic formula with the same error as the prime cases obtained in \cite{You11}. The two ingredients introduced in \cite{BFK+17a, BFK+17b} to improve the error term for prime moduli do not seem to work well for general moduli, while we have made no attempt to optimize the error term.

If only considers the leading main term such as results of Heath-Brown and Soundararjan, one can just bound all off-diagonal terms. However, to distinguish all lower-order main terms with a power saving in the error, we should prove an asymptotic formula for these off-diagonal terms. Comparing to the divisor problem \eqref{eqdivisor1} for the Riemann zeta-function, we should treat the following divisor sum
\begin{align}\label{eqdivisor2}
\frac1{\varphi^*(q)}\sum_{d\mid q}\varphi(d)\mu(q/d)\mathop{\sum\sum}_{\substack{m\equiv n(\bmod d)\\ (mn,q)=1}}\frac{d(m)d(n)}{\sqrt{mn}}V\left(\frac{mn}{q^2}\right),
\end{align}
where $V(x)$ is a smooth function with rapid decay on $x$, satisfying $V(x)\sim 1$ for small $x$. The diagonal terms $m=n$ contribute the leading term of the polynomial in \eqref{eqmain}.

The treatment of off-diagonal terms in \eqref{eqdivisor2} is much more difficult. One of the primary difficulties is the large ranges of the sum over $m$ and $n$. For the fourth moment of the Riemann zeta-function, the averaging over $t$-aspect forces $m$ and $n$ to be very close to each other. A similar phenomenon also arises for Dirichlet $L$-functions, which is mainly due to the cancelation of the M\"{o}bius function in \eqref{eqdivisor2} but is much weaker. Another primary difficulty is the coprime condition in the sum over $m$ and $n$, due to which some classical tools such as the Voronoi summation formula and the Estermann $D$-function are unavailable or less power.

Difficulties of these two aspects may be a little bit clear or can be removed when $q$ is prime, in which case $d$ just takes values $1$ or $q$ in \eqref{eqdivisor2}. When $m$ and $n$ are far away from each other, one may expect the cancelation due to the uniform distribution of the divisor function in arithmetic progressions, that is if $(m,q)=1$,
\begin{align}\label{eqdivisor3}
\sum_{\substack{n\le x\\ n\equiv m (\bmod q)}}d(n)=\frac1{\varphi(q)}\sum_{\substack{n\le x\\ (n,q)=1}}d(n)+(\text{error}).
\end{align}
The problem of estimating \eqref{eqdivisor3} is a challenging open problem, which has already been studied such as in \cite{Fou85, FI92}, but known results are still insufficient in estimating \eqref{eqdivisor2} for prime $q$. However, due to a more complicate cancelation of the M\"{o}bius function, it does not seem that any further progression on \eqref{eqdivisor3} may benefit the estimation \eqref{eqdivisor2} for general $q$. On the other hand, when $q$ is prime, the coprime condition $(mn,q)=1$ does not really work since it can be removed directly with a negligible error $q^{-1+\ve}$.

Let $B_{M,N}$ denote the same sum as in \eqref{eqdivisor2} but with $m, n$ restricted to the dyadic segments $M<m\le 2M$, $N<n\le 2N$. The treatment of $B_{M,N}$ proceeds essentially in two different ways, according to the relative size of $M$ with respect to $N$. Estimating $B_{M,N}$ for general $q$ is much more difficult, and many essential difficulties, especially from the coprime condition, should be handled with new ingredients. Nevertheless, the basic strategy closely follows Young.

Off-diagonal terms with  $M$ and $N$ relatively close contribute lower-order main terms. When $q$ is prime, the theory of the Estermann $D$-function plays a key role in the deduction of these lower-order main terms, but it does not work for general $q$, due to the coprime condition. Motivated by \cite{Wu19}, we deduce the theory for such special divisor sum function
\[
\mathcal{D}_{q}\left(s,\lambda,\frac h{l},r\right)=\sum_{\substack{(n,q)=1\\ (n+r,q)=1}}\frac{\sigma_\lambda(n)} {n^s} e\left(n\frac h{l}\right)
\]
with $\sigma_{\lambda}(n)=\sum_{d\mid n}d^\lambda$.
The deduction is elaborate and presented in Section \ref{secDq}. In addition, an application of the Kuznetsov formula for classical Kloosterman sums is insufficient to treat the error arising from the functional equation, instead, we should apply the Kuznetsov trace formula for general Kloosterman sums for congruence groups after the evolution of the original Kloosterman sum.

When $M$ and $N$ are far away from each other, such off-diagonal terms just contribute to the error term, and Weil bound for Kloosterman sums controls them well for $N/M>q^{1+\ve}$. To extend the range, we open the divisor function $d(n)$ and let $n=n_1n_2$, where $n_1\asymp N_1$ and $n_2\asymp N_2$ with $N_1\ll N_2$ and $N_1N_2\asymp N$. For large $N_1\gg q^{\frac12+\ve}$, we bound them by treating the following double sum in Kloosterman sums
\[
\mathfrak{D}(L,K;q)=\sum_{l\le L}\left|\mathop{\sum\nolimits^*}_{k\le K} \alpha_k e\left(\frac{al\ol{k}}{q}\right)\right|,
\]
which is of independent interest.
High power factors in $q$ cause difficulties in the treatment of such double sum, and the new ingredient is our applying exact formulae of Kloosterman sums for prime power moduli, by which we may reduce the sum of Kloosterman sums to a Gauss sum at last. The treatment relies on careful case analysis, and we present it in Section \ref{secdoublesum}. For small $N_1$, the functional equation for generalized Estermann $D$-function with a character allows for an extra saving by averaging over $m$.

\begin{remark}
In combination with \cite{Wu20}, we may extend the way to treat the fourth moments with $t$ also enclosed. This, in combination with \cite[Theorem 1.1]{Wu20}, may give an asymptotic formula with a power saving for the fourth moments with $t$ also enclosed, which will be shown in detail in a later version of \cite{Wu20}.
\end{remark}

\subsection{The shifted fourth moment}
We have chosen to calculate the shifted fourth moment of Dirichlet $L$-functions including the shifts $\alpha,\beta,\gamma,\delta$, because the shifts can split higher-order poles into several simple poles and thus allow for a clearer structure of the main terms. As mentioned by Young \cite{You11}, we may focus on the treatment of the even characters because the case of the odd characters is similar (see also \cite[Section 8.3]{You11} for necessary changes).
Thus, what we treat directly is
\begin{align}
M(\alpha,\beta,\gamma,\delta)=\frac2{\vp^*(q)}\mathop{{\sum}^+}_{\chi(\bmod q)} L\left(\tfrac12+\alpha,\chi\right) L\left(\tfrac12+\beta,\chi\right) L\left(\tfrac12+\gamma,\ol{\chi}\right) L\left(\tfrac12+\delta,\ol{\chi}\right),\notag
\end{align}
where the ``$+$" indicates that the summation is over all primitive even characters.

For convenience, we introduce some notations
\[
Z_q(\alpha,\beta,\gamma,\delta)=\frac{\zeta_q(1+\alpha+\gamma)\zeta_q(1+\alpha+\delta)\zeta_q(1+\beta+\gamma)\zeta_q(1+\beta+\delta)}{\zeta_q(2+\alpha+\beta+\gamma+\delta)}
\]
and
\begin{align}\label{eqX}
X_{\alpha,\gamma}=X_{\alpha}X_{\gamma},\ \ \ \ \ \ \ \ X_{\alpha,\beta,\gamma,\delta}=X_{\alpha}X_{\beta}X_{\gamma}X_{\delta}
\end{align}
with
\begin{align}\label{eqdefX}
X_{\alpha}=\left(\frac{q}{\pi}\right)^{-\alpha}\frac{\Gamma\left(\frac{\frac12-\alpha+\mathfrak{a}}{2}\right)} {\Gamma\left(\frac{\frac12+\alpha+\mathfrak{a}}{2}\right)} \ \ \ \ \ \text{for}\ \ \ \mathfrak{a}=0,1.
\end{align}
Note that $Z_q(\alpha,\beta,\gamma,\delta)$ is invariable under transpositions $\alpha\leftrightarrow\beta$ as well as $\gamma\leftrightarrow\delta$, and it is also invariable under the two transpositions $\alpha\leftrightarrow\gamma$, $\beta\leftrightarrow\delta$ simultaneous (e.g., $Z_q(-\gamma,-\delta,-\alpha,-\beta)=Z_q(-\alpha,-\beta,-\gamma,-\delta)$). In addition, $X_{\alpha,\gamma}$ and $X_{\alpha,\beta,\gamma,\delta}$ are invariable under any transpositions of the shifts.
\begin{conjecture}[CFK+05]\label{conjecture}
For any integer $q\nequiv 2 \pmod{4}$, and with shifts $\ll(\log q)^{-1}$, we have
\begin{align}\label{eqconjecture}
M(\alpha,\beta,\gamma,\delta)&=Z_{q}(\alpha,\beta,\gamma,\delta) +X_{\alpha,\beta,\gamma,\delta}Z_{q}(-\gamma,-\delta,-\alpha,-\beta)+X_{\alpha,\gamma}Z_{q}(\beta,-\gamma,\delta,-\alpha)\\
+&X_{\beta,\gamma}Z_{q}(\alpha,-\gamma,\delta,-\beta)
+X_{\alpha,\delta}Z_{q}(\beta,-\delta,\gamma,-\alpha) +X_{\beta,\delta}Z_{q}(\alpha,-\delta,\gamma,-\beta)+O\left(q^{-\frac12+\ve}\right)\notag
\end{align}
with $\mathfrak{a}=0$.
\end{conjecture}
There is a similar conjecture for the odd character case, and the only difference is to specify $\mathfrak{a}=1$ in \eqref{eqdefX}, while $\mathfrak{a}=0$ for the even character case.

\begin{theorem}\label{thmmains}
Conjecture \ref{conjecture} holds but with an error of size $O\left(q^{-\frac1{80}+\frac{\theta}{40}}\right)$.
\end{theorem}

The symmetry implies that the right-hand side of \eqref{eqconjecture} is holomorphic concerning the shifts, which has been proved in a more general setting in \cite[Lemma 2.5.5]{CFK+05}. Thus, taking the limit as the shifts go to zero in \eqref{eqconjecture} gives Theorem \ref{thmmain}.

Note that $M(\alpha,\beta,\gamma,\delta)$ and all main terms in the right-hand side of \eqref{eqconjecture} are holomorphic concerning the shifts, and so is the error. Thus, the maximum modulus principle implies that it is sufficient to prove Theorem \ref{thmmains} for such special shifts that,
each of shifts lies in a fix annulus with inner and outer radii $\asymp (\log q)^{-1}$, which are separated enough so that $|\alpha\pm\beta|\gg(\log q)^{-1}$, etc.

\subsection{Notation}
We present some notations here for ease of reference.
As usual, we use $\ve$ to denote an arbitrarily small positive constant that may vary from line to line, and $e(x)$ means the exponential function $e^{2\pi ix}$. $\varphi(n)$ is the Euler function, which denotes the number of positive integers less than $n$ and coprime with $n$. The notations $\sigma_\lambda(n)$ and $\sigma_{\alpha,\beta}(n)$ are defined as
\[
\sigma_\lambda(n)=\sum_{d\mid n}d^\lambda,\ \ \ \ \sigma_{\alpha,\beta}(n)=\sum_{ad=n}a^{-\alpha} d^{-\beta},
\]
and it is obvious
\[
\sigma_{\alpha,\beta}(n)=n^{-\alpha}\sigma_{\alpha-\beta}(n).
\]

The symbol $\theta$ always denotes an admissible exponent for the Ramanujan-Petersson conjecture for Maass newforms, and the current best-known result is $\theta=7/64$.

In this paper, $W(x)$ always denotes a smooth non-negative function compactly supported in $[1,2]$, which may have different expressions for each occurrence, and we apply $\wt{W}(u)$ and $\wh{W}(u)$ to denote its Mellin transform and Fourier transform respectively.

For notational convenience, we apply two symbols related to $q$. We use $q^*$ to denote the largest square-free divisor of $q$, and $q_k$ to be the largest divisor of $q$ coprime with $k$. Note that $q_0, q_1, q_2, \dots$ are not the case, where the subscripts are just serial numbers.

In the following, functions we face are often on the shifts $\alpha,\beta,\gamma,\delta$, such as $M(\alpha,\beta,\gamma,\delta)$. In the case of no ambiguity, we frequently omit the shifts for convenience when they are in the order $\alpha,\beta,\gamma,\delta$.

We apply the following especial summation notations
\[
\mathop{\sum\nolimits^*}_{n (\bmod q)}\ \ \ \ \text{and} \ \ \ \ \mathop{\sum\nolimits^*}_{\chi (\bmod q)},
\]
where the first is a sum over such $n$ coprime with $q$, and the second is over all primitive characters modulo $q$.

\section{Background and auxiliary lemmas}
\subsection{The Dirichlet $L$-function and the Hurwitz zeta-function}
Let $q$ be a positive integer, and $\chi$ be a Dirichlet character modulo $q$. For $\text{Re}(s)>1$, the Dirichlet $L$-function may be defined as
\[
L(s,\chi)=\sum_{n=1}^\infty\chi(n)n^{-s}.
\]
The Dirichlet $L$-function satisfies the functional equation
\[
\left(\frac q\pi\right)^{\frac s2}\Gamma\left(\frac{\frac12+s+\mathfrak{a}}{2}\right)L\left(\tfrac12+s,\chi\right)=i^{-\mathfrak{a}}q^{-\frac12}\tau(\chi)\left(\frac q\pi\right)^{-\frac s2}\Gamma\left(\frac{\frac12-s+\mathfrak{a}}{2}\right) L\left(\tfrac12-s,\ol{\chi}\right),
\]
where
\begin{align}
\tau(\chi)=\sum_{x(\bmod q)}\chi(x)e\left(\frac{x}q\right),\ \ \ \ \ \ \ \mathfrak{a}=\left\{
\begin{aligned}
&0, \ \ \ \ \chi(-1)=1,\\
&1, \ \ \ \ \chi(-1)=-1.\notag
\end{aligned}
\right.
\end{align}

For $\text{Re}(s)>1$ and $0<x\le1$, the Hurwitz zeta-function is defined as
\[
\zeta(s,x)=\sum_{n=0}^\infty(n+x)^{-s}.
\]
The functional equation for the Hurwitz zeta-function is
\[
\zeta(1-s,x)=(2\pi)^{-s}\Gamma(s)\left(e\left(-\frac{\pi s}{2}\right)F(s,x)+e\left(\frac{\pi s}{2}\right)F(s,-x)\right),
\]
where
\[
F(s,x)=\sum_{n=1}^{\infty} e(nx)n^{-s}
\]
for $\text{Re}(s)>1$.
\subsection{Approximate functional equation}
To expand $M(\alpha,\beta,\gamma,\delta)$, we apply the following approximate functional equation, deduced from the functional equation of $L(s,\chi)$ standardly, see also \cite[Proposition 2.4]{You11}.
\begin{lemma}[Approximate functional equation]
\label{lemafe}
Let $G(s)$ be an even, entire function of exponential decay in any strip $|\text{Re}(s)|<C$ satisfying $G(0)=1$, and let
\begin{align}
\label{defV}
V_{\alpha,\beta,\gamma,\delta}(x)=\frac1{2\pi i}\int_{(1)}\frac{G(s)}{s}g_{\alpha,\beta,\gamma,\delta}(s)x^{-s}ds,
\end{align}
where
\begin{align}
\label{defg}
g_{\alpha,\beta,\gamma,\delta}(s)=\pi^{-2s}\frac{\Gamma\left(\frac{ \frac12+\alpha+s+\mathfrak{a}}2\right) \Gamma\left(\frac{ \frac12+\beta+s+\mathfrak{a}}2\right) \Gamma\left(\frac{ \frac12+\gamma+s+\mathfrak{a}}2\right) \Gamma\left(\frac{ \frac12+\delta+s+\mathfrak{a}}2\right)}{\Gamma\left(\frac{\frac12+\alpha+\mathfrak{a}}2\right) \Gamma\left(\frac{\frac12+\beta+\mathfrak{a}}2\right) \Gamma\left(\frac{\frac12+\gamma+\mathfrak{a}}2\right) \Gamma\left(\frac{\frac12+\delta+\mathfrak{a}}2\right)}.
\end{align}
Furthermore, set
\begin{align}
\widetilde{V}_{\alpha,\beta,\gamma,\delta}(x)= X_{-\gamma,-\delta,-\alpha,-\beta}V_{\alpha,\beta,\gamma,\delta}(x)\notag
\end{align}
with $X$ defined by \eqref{eqX}.
Then
\begin{align}
L\left(\tfrac12+\alpha,\chi\right) &L\left(\tfrac12+\beta,\chi\right) L\left(\tfrac12+\gamma,\ol{\chi}\right) L\left(\tfrac12+\delta,\ol{\chi}\right)\notag\\
&=\sum_{m,n}\frac{\sigma_{\alpha,\beta}(m)\sigma_{\gamma,\delta}(n) \chi(m)\ol{\chi}(n)} {(mn)^{\frac12}}V_{\alpha,\beta,\gamma,\delta}\left(\frac{mn}{q^2}\right)\notag\\ &\ \ \ + \sum_{m,n}\frac{\sigma_{-\gamma,-\delta}(m)\sigma_{-\alpha,-\beta}(n) \chi(m)\ol{\chi}(n)} {(mn)^{\frac12}}\wt{V}_{-\gamma,-\delta,-\alpha,-\beta}\left(\frac{mn}{q^2}\right).\notag
\end{align}
\end{lemma}

This approximate functional equation holds for general $G$.  However, we will apply it with a special $G$ in the following.
\begin{definition}[Definition of $G(s)$]\label{DefG}
Let $G(s)=P_{\alpha,\beta,\gamma,\delta}(s)\exp(s^2)$, where $P_{\alpha,\beta,\gamma,\delta}(s)$ is a even polynomial in $s$ satisfying some properties: it is rational, symmetric, and even in the shifts; $P_{\alpha,\beta,\gamma,\delta}(0)=1$ and $P_{\alpha,\beta,\gamma,\delta}\left(-\frac{\alpha+\gamma}2\right)=0$, $P_{\alpha,\beta,\gamma,\delta}\left(\frac12\pm\alpha\right)=0$ (as well as $\frac{\beta+\delta}2, \frac12\pm\beta$, etc., by symmetry).
\end{definition}




\subsection{The orthogonality formula}
We apply the following orthogonality formula to remove the sum over characters, whose proof may refer to many sources, such as \cite{HB81} and \cite{Sou07}.
\begin{lemma}[The orthogonality formula]
\label{lemof}
For $(mn,q)=1$, we have
\begin{align}
\label{lemof2}
\mathop{{\sum}^+}_{\chi(\bmod q)}\chi(m)\ol{\chi}(n) =\frac12\left(\sum_{d\mid(q,m-n)}\vp(d)\mu\left(\frac qd\right) +\sum_{d\mid(q,m+n)}\vp(d)\mu\left(\frac qd\right)\right).
\end{align}
\end{lemma}




\subsection{Additive characters}
Frequently, we use additive characters to remove the coprime condition as well as the divisibility condition in sums, which we present here for ease of reference.
\begin{lemma}\label{lemcoprimec}
Let $a(n)$ be an arithmetic function. For any integer $k$, we have
\begin{align}\label{eqdivisor}
\sum_{k\mid n}a(n)=\frac1k\sum_{d\mid k}\sum_nc_{d}(n)a(n),
\end{align}
where $c_d(n)$ is the Ramanujan sum.
\end{lemma}
\begin{proof}
Applying the sum of divisors formula for the Ramanujan sum and reversing the order of the summation, we have
\begin{align}
\sum_{d\mid k}c_{d}(n)=&\sum_{d\mid k}\sum_{l\mid (d,n)}l\mu\left(\frac {d}{l}\right)\notag\\
=&\sum_{l\mid (k,n)}l\sum_{m\mid \frac kl}\mu(m)=\left\{
\begin{aligned}
&k \ \ \ \text{for}\ \ \ k\mid n,\\
&0 \ \ \ \text{otherwise}.\notag
\end{aligned}
\right.
\end{align}
Then \eqref{eqdivisor} follows directly from this.
\end{proof}


\begin{lemma}
Let $a(n)$ be an arithmetic function. For any integer $q$, we have
\begin{align}\label{eqcoprimee}
\sum_{(n,q)=1}a(n)=\sum_{d\mid q}\frac{\mu(d)}{d}\sum_{r (\bmod d)} e\left(\frac{rn}{d}\right)a(n)
\end{align}
and
\begin{align}\label{eqcoprime}
\sum_{(n,q)=1}a(n)=\sum_{d\mid q}\frac{\mu(d)}{d}\sum_{d_1\mid d}\sum_nc_{d_1}(n)a(n).
\end{align}
\end{lemma}
\begin{proof}
We first come to \eqref{eqcoprimee}. Note that the sum over $r$ is $0$ unless $d\mid n$, in which case it equals $d$. Hence
\begin{align}\label{eqq_0mid}
\sum_{r (\bmod d)} e\left(\frac{rn}{d}\right)a(n)=d\sum_{d\mid n}a(n).
\end{align}
Applying this into the right-hand side of \eqref{eqcoprimee}, we find that it evolves into
\begin{align}
\sum_{d\mid q}\mu(d)\sum_{d\mid n}a(n)=\sum_na(n)\sum_{d\mid(q,n)}\mu(d)=\sum_{(n,q)=1}a(n).\notag
\end{align}
The same discussion but with \eqref{eqq_0mid} replaced by \eqref{eqdivisor} gives \eqref{eqcoprime}.
\end{proof}


\subsection{Multiplicative characters}
Sometimes, especially in the deduction of an asymptotic functional equation for the $\mathcal{D}_q$-function, we need to apply multiplicative characters to remove the coprime condition in the sum.
\begin{lemma}\label{lemcoprimechie}
Let $a(n)$ be an arithmetic function. For any integers $r$ and $q$, we have
\begin{align}\label{eqcoprimechi}
\sum_{\substack{(n,q)=1\\ (n+r,q)=1}}a(n)=\sum_{b\mid q}\frac{\mu(b)}{\vp(b)} \sum_{\chi(\bmod b)}\ol{\chi}(-r)\sum_{n}\chi(n)\chi_0(n)a(n),
\end{align}
where $\chi_0$ is the principal character modulo $q$.
\end{lemma}
\begin{proof}
Note that $\ol{\chi}(-r)=0$ when $(b,r)\neq1$. The right-hand side of \eqref{eqcoprimechi} equals
\[
\sum_{\substack{b\mid q\\ (b,r)=1}}\frac{\mu(b)}{\vp(b)} \sum_{\chi(\bmod b)}\ol{\chi}(-r)\sum_{n}\chi(n)\chi_0(n)a(n).
\]
For $(b,r)=1$, it is well-known that
\begin{align}
\sum_{\chi(\bmod b)}\ol{\chi}(-r)\chi(n)=
\left\{
\begin{aligned}
&\vp(b)\ \ \ \ \text{for} \ \ \ \ b\mid n+r,\\
&0 \ \ \ \ \ \ \ \ \ \text{otherwise}.\notag
\end{aligned}
\right.
\end{align}
So the right-hand side of \eqref{eqcoprimechi} evolves into
\[
\sum_{\substack{b\mid q\\ (b,r)=1}}\mu(b)\sum_{b\mid n+r}\chi_0(n)a(n)=\sum_{(n,q)=1}a(n)\sum_{\substack{b\mid q\\ b\mid n+r}}\mu(b)=\sum_{\substack{(n,q)=1\\ (n+r,q)=1}}a(n),
\]
which establishes the lemma.
\end{proof}
\subsection{The reciprocity law}
To split and combine the exponential function, we will frequently apply the reciprocity law, which can be checked directly by the Chinese Remainder Theorem.
\begin{lemma}\label{lemmaspe(x)}
Let $a,b$ be two coprime integers. For
\begin{align}
\left\{
\begin{aligned}
&x\equiv y \pmod a,\\
&x\equiv z \pmod b,\\
\end{aligned}
\right.\notag
\end{align}
we have
\begin{align}
e\left(\frac x{ab}\right)=e\left(\frac{y\ol{b}}{a}\right)e\left(\frac{z\ol{a}}{b}\right).\notag
\end{align}
\end{lemma}

\subsection{An arithmetic formula}
In some ways, the following arithmetic formula may reflect the cancelation of M\"{o}bius function in \eqref{eqdivisor2}.
\begin{lemma}\label{lemsumd}
For any integer $q$, we have
\begin{align}
\sum_{d\mid q}\frac{\varphi(d)}{d}\mu\left(\frac qd\right)=\frac{\mu(q)}q.\notag
\end{align}
\end{lemma}
\begin{proof}
Note that both sides are multiplicative functions on $q$, so it is sufficient to check the identity for prime powers.
For $q=p^m$ with $m\ge1$, we have
\begin{align}
\sum_{d\mid p^m}\frac{\varphi(d)}{d}\mu\left(\frac {p^m}d\right)&=\frac{\varphi(p^m)}{p^m}-\frac{\varphi(p^{m-1})}{p^{m-1}}\notag\\
&=\left\{
\begin{aligned}
&0\ \ \ \ \ \ \ \ \ \text{for}\ \ \ \ m>1,\\
&-\frac1p\ \ \ \ \text{for}\ \ \ \ m=1,\notag
\end{aligned}
\right.\notag
\end{align}
completing the proof of the lemma.
\end{proof}

\subsection{The Estermann  $D$-function}
The Estermann $D$-function is defined as
\begin{align}
D\left(s,\lambda,\frac hl\right)=\sum_{n}\frac{\sigma_\lambda(n)}{n^s}e\left(n\frac hl\right).\notag
\end{align}


\begin{lemma}[The Estermann $D$-function]\label{lemED}
Let $h,l$ be two coprime integers. For any given $\lambda\in\mathbb{C}$, $D(s,\lambda,\frac hl)$ is meromorphic as a function of $s$, and satisfies the functional equation
\begin{align}
D\left(\tfrac12+s,\lambda,\frac hl\right)=&2(2\pi)^{-1-\lambda+2s}\Gamma\left(\tfrac12-s\right) \Gamma\left(\tfrac12+\lambda-s\right)l^{\lambda-2s}\notag\\
&\times\left[D\left(\tfrac12-s,-\lambda,-\frac{\ol{h}}{l}\right)\sin\left(\pi\left(s-\frac \lambda 2\right)\right) +D\left(\tfrac12-s,-\lambda,\frac{\ol{h}}{l}\right)\cos\left(\frac{\pi\lambda}{2}\right)\right].\notag
\end{align}
If $\lambda\neq0$ then $D$ has simple poles at $s=1$ and $s=1+\lambda$ with respective residues
\begin{align}
l^{-1+\lambda}\zeta(1-\lambda),\ \ \ \ \ \ \ l^{-1-\lambda}\zeta(1+\lambda).\notag
\end{align}
\end{lemma}

The lemma was observed by Estermann \cite{Est30}, whose proof may be found in \cite[Lemma 3.7]{Mot97}.







\subsection{A dyadic partition of unity}
We also require a dyadic partition of unity. Let $W$ be a smooth non-negative function compactly supported in $[1,2]$ such that, for any $x\in\mathbb{R}^+$,
\begin{align}
\sum_{M}W\left(\frac xM\right)=1.\notag
\end{align}
Here $M$ varies over a set of positive real numbers, with the number of such $M$ less than $X$ being $O(\log X)$. The $W$ function has the Mellin pair
\begin{align}
\left\{
\begin{aligned}
&\wt{W}(u)=\int_0^\infty W(x)x^{u-1}dx,\\
&W(x)=\frac1{2\pi i}\int_{(c_u)}\wt{W}(u)x^{-u}du.
\end{aligned}
\right.\notag
\end{align}
By taking products of the $W\left(\frac xM\right)$, we then create the partition of unity $\left\{W\left(\frac xM\right)W\left(\frac yN\right)\right\}$ on $\mathbb{R}^+\times\mathbb{R}^+$.

\subsection{The classical Kloosterman sums}
In this section, we present some materials for the classical Kloosterman sum, most of which can be found in Section 4.3 of \cite{Iwa97}. The classical Kloosterman sum is defined as
\begin{align}\label{eqdefKSC}
S(m,n;q)=\mathop{\sum\nolimits^*}_{d(\bmod q)}e\left(\frac{m\ol{d}+nd}{q}\right),
\end{align}
which meets the following elementary properties:
\begin{align}\label{eqKSjh}
S(m,n;q)=S(n,m;q),
\end{align}
\begin{align}
S(cm,n;q)=S(m,cn;q)\ \ \ \ \text{if}~(c,q)=1,
\end{align}
\begin{align}\label{eqexpandKS}
S(m,n;q)=\sum_{d\mid (m,n,q)}dS(mnd^{-2},1;qd^{-1}).
\end{align}
The Kloosterman sum reduces to the Ramanujan sum if one of $m, n$ is equal to $0 \pmod q$.
It is well-known by Weil bound that
\[
S(m,n;q)\ll (m,n,q)^{\frac12}q^{\frac 12}d(q).
\]

When $q=p^j$ with $j\ge2$, we can deduce exact expressions for Kloosterman sums.
\begin{lemma}\label{lemKS0}
If $q=p^j$ with $j\ge2$ and $q\neq(m,q)\neq (n,q)\neq q$, then
\[
S(m,n;q)=0.
\]
\end{lemma}
\begin{proof}
We first treat the special case with $(m,q)>1$ and $n=1$. Put $d=bp^{j-1}+a$ in \eqref{eqdefKSC}, where $a$ ranges $\bmod p^{j-1}$, $a$ is prime to $p$, and $b$ ranges freely $\bmod p$. Then
\begin{align}
S(m,1;q)=\mathop{\sum\nolimits^*}_{a(\bmod p^{j-1})}\sum_{b(\bmod p)}e\left(\frac{m\ol{(bp^{j-1}+a)}}{p^j}\right)e\left(\frac{bp^{j-1}+a}{p^j}\right).\notag
\end{align}
Note that
\[
e\left(\frac{m\ol{(bp^{j-1}+a)}}{p^j}\right)=e\left(\frac{m\ol{a}}{p^j}\right)
\]
for $p\mid m$.
Consequently, the sum over $b$ vanishes, and thus $S(m,1;q)=0$. Then the lemma follows easily from this in combination with properties \eqref{eqKSjh} and \eqref{eqexpandKS}.
\end{proof}

For the case $(m,q)=(n,q)=(m,n,q)$, we can extract $(m,n,q)$ from $S(m,n;q)$ by the following lemma.
\begin{lemma}\label{lemKSrelation}
If $q=p^j$ with $j\ge2$ and $(m,n,q)\neq q$, we have
\begin{align}\label{eqKSdivisor}
S(m,n;q)=(m,n,q)S(m',n';q'),
\end{align}
where $m'=m/(m,n,q)$, $n'=n/(m,n,q)$, and $q'=q/(m,n,q)$.
\end{lemma}
\begin{proof}
We apply \eqref{eqexpandKS} to extract $(m,n,q)$ from $S(m,n;q)$, and then by Lemma \ref{lemKS0}, all terms except the one with $d=(m,n,q)$ vanish. This only term gives the right-hand side of \eqref{eqKSdivisor}, establishing the lemma.
\end{proof}

Thus, to obtain an exact expression for $S(m,n,q)$ with $q=p^j$, $j\ge2$, it is sufficient to treat the case $(mn,q)=1$, which has been deduced by Iwaniec \cite{Iwa97}.
\begin{lemma}[Iwaniec, {\cite[p. 60]{Iwa97}}]\label{lemKScoprime}
If $q=p^j$ with $j\ge2$ and $(q,2mn)=1$, we have $S(m,n;q)=0$ unless $m\equiv l^2n \pmod q$, in which case there is
\begin{align}\label{eqKSexact}
S(m,n;q)=\left(\frac{ln}{q}\right)q^{\frac12}\left\{\ve_q e\left(\frac{2ln}{q}\right)+\ol{\ve}_q e\left(-\frac{2ln}{q}\right)\right\},
\end{align}
where $\left(\frac{ln}{q}\right)$ is the Legendre-Jacobi symbol and $\ve_q=1,i$ according to whether $q\equiv 1$ or $-1 \pmod 4$. There is a similar result for $q=2^j$ with $j\ge2$.
\end{lemma}

\subsection{The Gauss sum}
For any real number $\alpha, \beta$ with $\alpha\neq0$, a general Gauss sum is given by
\[
\sum_{1\le n\le N}e(\alpha n^2+\beta n).
\]
In the special case where $a,b$ are integers with $(2a,N)=1$, we have
\[
\left|\sum_{1\le n\le N}e\left(\frac{an^2+bn}{N}\right)\right|=\sqrt{N}.
\]
Though there is no simple expression for a general Gauss sum, it has the following bound.
\begin{lemma}[Iwaniec and Kowalski, {\cite[Theorem 8.1]{IK04}}]\label{lemGauss}
We have
\[
\sum_{1\le n\le N}e(\alpha n^2+\beta n)\le 2Nq^{-\frac12}+q^{\frac12}\log q,
\]
where $q$ is decided by the rational approximation to $2\alpha$ of type
\[
\left|2\alpha-\frac aq\right|\le\frac{1}{2Nq}
\]
with $(a,q)=1$ and $1\le q\le 2N$.
\end{lemma}


\subsection{Automorphic preliminaries}
We sketch some elementary materials we require on automorphic forms to apply the Kuznetsov formula in this section, which mainly follows from \cite{BM15}, \cite{Iwa95}, and \cite{DI82}. We write the Fourier expansion of a holomorphic modular form $f$ of level $Q$ and weight $k$ as
\[
f(\sigma_{\mathfrak{a}}z)=\sum_{n\ge1}\rho_f(\mathfrak{a},n)(4\pi n)^{\frac k2}e(nz),
\]
and similarly write for a Maass  form $f$ of level $Q$ and spectral parameter $\kappa=\kappa_f\in\mathbb{R}\cup[-i\theta,i\theta]$ as
\[
f(\sigma_{\mathfrak{a}}z)=\sum_{n\neq0}\rho_f(\mathfrak{a},n)W_{0,i\kappa}(4\pi|n|y)e(nx),
\]
where $W_{0,i\kappa}(y)=(y/\pi)^{1/2}K_{i\kappa}(y/2)$ is a Whittaker function. For cusps $\mathfrak{a}, \mathfrak{c}$ of $\Gamma_0(Q)$, the Fourier expansion of the Eisenstein series $E_{\mathfrak{c}}(\sigma_\mathfrak{a} z,s)$ at $s=\frac12+i\kappa$ is written as
\[
E_{\mathfrak{c}}(\sigma_\mathfrak{a} z, \tfrac12+i\kappa)=\delta_{\mathfrak{a},\mathfrak{c}}y^{1/2+i\kappa}+\varphi_{\mathfrak{a},\mathfrak{c}}(\tfrac12+i\kappa)y^{1/2-i\kappa} +\sum_{n\neq0}\rho_{\mathfrak{a},\mathfrak{c}}(n,\kappa)W_{0,i\kappa}(4\pi|n|y)e(nx).
\]
For the special case $\mathfrak{a}=\infty$, we denote the Fourier coefficients as
\[
\rho_f(n)=\rho_f(\infty,n), \ \ \ \ \ \ \ \ \ \ \ \ \rho_{\mathfrak{c}}(n,\kappa)=\rho_{\infty,\mathfrak{c}}(n,\kappa)
\]
respectively.
If $f$ is a cuspidal newform with normalized Hecke eigenvalues denoted by $\lambda_f(n)$, there is
\[
\lambda_{f}(n)\rho_f(1)=\sqrt{n}\rho_f(n).
\]
Estimating the size of $\lambda_f(n)$ is a significant problem, and is generally believed as
\[
|\lambda_f(n)|\le d(n),
\]
which is known as the Ramanujan-Petersson conjecture. The best result known now is
\begin{align}
\lambda_f(n)\ll n^{\theta+\ve}\notag
\end{align}
with $\theta=7/64$, due to Kim and Sarnak \cite{Kim03}.

Since not all cusp forms are newforms, a $L^2$-basis $\mathcal{B}(Q)$ for $\mathcal{A}(Q)$, the space of Maass forms of level $Q$, will in general also include oldforms. The same phenomenon also occurs in a $L^2$-basis $\mathcal{B}_k(Q)$ for $S_k(Q)$, the space of holomorphic cusp forms of weight $k$ and level $Q$. However, one may apply the newform theory to deduce an orthogonal basis based on normalized newforms. Taking Maass forms as an example,
one may find a special orthogonal basis $\mathcal{B}(Q)$ that, for each $f\in \mathcal{B}(Q)$, there is a normalized newform $f^*$ of level $Q_1$ with $Q_1\mid Q$ and of the same spectral parameter (see also \cite[Section 5]{BM15}). In particular,
if $m=bm'$ with $(m',b)=1$, then
\begin{align}\label{eqMrho}
\sqrt{m}\rho_f(m)=\sum_{d\mid (Q,b/(b,Q))}\mu(d)\chi_0(d)\lambda_{f^*}\left(\frac{b}{d(b,Q)}\right)\left(\frac{(b,Q)m'}{d}\right)^{\frac12}\rho_f\left(\frac{(b,Q)m'}{d}\right),
\end{align}
where $f^*$ is the underlying newform and $\chi_0$ is the trivial character modulo $Q_1$ (see also \cite[formula (5.2)]{BM15}).
Moreover, if $f$ is a holomorphic forms, $f^*$ satisfies the Ramanujan-Petersson conjecture. If $b\in\mathbb{N}$ and $a_m$ is any finite sequence of complex numbers supported on the integers $m=bm'$ with $(m',b)=1$, then
\begin{align}\label{eqhrho}
\left|\sum_m a_m\sqrt{m}\rho_f(m)\right|^2\le d(b)^2\sum_{d\mid(b,Q)}\left|\sum_{m'}a_{bm'}\sqrt{dm'}\rho_f(dm')\right|^2
\end{align}
(see also \cite[formula (5.4)]{BM15}).
There is a similar result for the coefficients of Eisenstein series that
\begin{align}\label{eqErho}
\sum_{\mathfrak{c}}\bigg|\sum_m a_m\sqrt{m}\rho_{\mathfrak{c}}(m,\kappa)\bigg|^2\le 9d(Q)^3d(b)^4\sum_{d\mid(b,Q)}\sum_{\mathfrak{c}}\bigg|\sum_{m'} a_{bm'}\sqrt{dm'}\rho_{\mathfrak{c}}(dm',\kappa)\bigg|^2
\end{align}
(see also \cite[ formula (5.5)]{BM15}).


\subsection{Kloosterman sums for congruence groups and the Kuznetsov formula}
In this section, we present some materials we require on Kloosterman sums for congruence groups, which can refer to \cite{DI82}.
Let $\mathfrak{a}$ and $\mathfrak{b}$ be two cusps of the Hecke congruence group $\Gamma_0(Q)$ and let $\sigma_{\mathfrak{a}}$ in $PSL(2,\mathbb{R})$ be a scaling matrix such that
\[
\sigma_{\mathfrak{a}}\infty =\mathfrak{a}\ \ \ \ \ \text{and}\ \ \ \ \sigma_{\mathfrak{a}}^{-1}T_{\mathfrak{a}}\sigma_{\mathfrak{a}}=\Gamma_{\infty},
\]
and similarly $\sigma_{\mathfrak{b}}$, where $\Gamma_{\mathfrak{a}}$ is the stabilizer of $\mathfrak{a}$, defined as $\Gamma_{\mathfrak{a}}=\{\sigma\in\Gamma_0(Q);\sigma\mathfrak{a}=\mathfrak{a}\}$. For any $m,n\in\mathbb{Z}$ and $\gammaup\in\mathbb{R}$, for which there exists
$\begin{pmatrix}
     \alphaup & \betaup \\
     \gammaup & \deltaup
\end{pmatrix}$
in $\sigma_{\mathfrak{a}}^{-1}\Gamma_0(Q)\sigma_{\mathfrak{b}}$, we define the Kloosterman sum as
\begin{align}\label{eqdefS}
S_{\mathfrak{a},\mathfrak{b}}(m,n;\gamma)=\mathop{\sum\nolimits^{'}}_{\deltaup(\bmod {\gammaup\mathbb{Z}})}e\left(m\frac{\alphaup}{\gammaup}+n\frac{\deltaup}{\gammaup}\right),
\end{align}
where the sum is over the $\deltaup$'s, taken modulo $\gamma\mathbb{Z}$, for which there exist $\alpha$ and $\beta$ such that $\begin{pmatrix}
     \alphaup & \betaup \\
     \gammaup & \deltaup
\end{pmatrix}$
in $\sigma_{\mathfrak{a}}^{-1}\Gamma_0(Q)\sigma_{\mathfrak{b}}$. This definition does not rely on choices of the scaling matrices, more precisely, different scaling matrices only affect the sum by multiplying a complex number with modulus $1$. For the special case $Q=1$, $\mathfrak{a}=\mathfrak{b}=\infty$ and $\sigma_{\infty}=1$, the Kloosterman sum $S_{\infty,\infty}(m,n;\gammaup)$ reduces to a classical Kloosterman sum.

An interesting example introduced in \cite{DI82} is a particular we focus on. Let $Q=\tau s$ with $(\tau,s)=1$, and we consider $S_{\infty,1/s}(m,n;\gammaup)$ with $\sigma_{1/s}=\begin{pmatrix}
     \sqrt{\tau} & 0 \\
     s\sqrt{\tau} & \frac1{\sqrt{\tau}}
\end{pmatrix}$.
From the definition \eqref{eqdefS}, one notes that
$S_{\infty,1/s}(m,n;\gammaup)$ is defined if and only if, $\gammaup$ may be written as $\gammaup=s\sqrt{\tau}C$ with integer $C$ coprime with $\tau$, and in this case there is
\begin{align}\label{eqSS}
S_{\infty,1/s}(m,n;\gammaup)=e\left(n\frac{\ol{s}}{\tau}\right)S(m\ol{\tau},n;sC),
\end{align}
where $s\ol{s}\equiv1 \pmod{\tau}$ and $\tau\ol{\tau}=1 \pmod {sC}$
(see also \cite[formula (1.6)]{DI82}).

To estimate a sum of Kloosterman sums of type \eqref{eqdefS}, we require the following Kuznetsov trace formula (see also \cite[Theorems 9.4, 9.5, 9.7]{Iwa95}).
\begin{lemma}[Kuznetsov formula]
Let $m,n$ be two positive integers and $\phi$ be a function $\mathcal{C}^2$ class on $[0,\infty)$, satisfying $\phi(0)=0$, $\phi^{(j)}(x)\ll(1+x)^{-2-\ve}$ for $j=0,1,2$; let $\mathfrak{a}$ and $\mathfrak{b}$ be two cusps of  $\Gamma=\Gamma_0(Q)$; denoting by $\sum\limits^{\Gamma}$ a summation performed over the positive real numbers $\gammaup$ for which $S_{\mathfrak{a},\mathfrak{b}}(m,n;\gammaup)$ is defined, one has
\begin{align}
\sum^{\Gamma}\frac1{\gammaup}S_{\mathfrak{a},\mathfrak{b}}(m,n;\gammaup)\phi\left(\frac{4\pi\sqrt{mn}}{\gammaup}\right)=&\sum_{2\le k\equiv0(\bmod 2)}\sum_{f\in \mathcal{B}_k(Q)}\Gamma(k)\tilde{\phi}(k)\sqrt{mn}\ol{\rho}_f(\mathfrak{a},m)\rho_f(\mathfrak{b},n)\notag\\
&+\sum_{f\in\mathcal{B}(Q)}\hat{\phi}(\kappa_f)\frac{\sqrt{mn}}{\cosh(\pi\kappa_f)}\ol{\rho}_f(\mathfrak{a},m)\rho_f(b,n)\notag\\
&+\frac1{4\pi}\sum_{\mathfrak{c}}\int_{-\infty}^{\infty}\hat{\phi}(\kappa)\frac{\sqrt{mn}}{\cosh(\pi\kappa)}\ol{\rho}_{\mathfrak{a},\mathfrak{c}}(m,\kappa) \rho_{\mathfrak{b},\mathfrak{c}}(n,\kappa) d\kappa\notag
\end{align}
and
\begin{align}
\sum^{\Gamma}\frac1{\gammaup}S_{\mathfrak{a},\mathfrak{b}}(m,-n;\gammaup)\phi\left(\frac{4\pi\sqrt{mn}}{\gammaup}\right)= &\sum_{f\in\mathcal{B}(Q)}\breve{\phi}(\kappa_f)\frac{\sqrt{mn}}{\cosh(\pi\kappa_f)}\ol{\rho}_f(\mathfrak{a},m)\rho_f(b,-n)\notag\\
&+\frac1{4\pi}\sum_{\mathfrak{c}}\int_{-\infty}^{\infty}\breve{\phi}(\kappa)\frac{\sqrt{mn}}{\cosh(\pi\kappa)}\ol{\rho}_{\mathfrak{a},\mathfrak{c}}(m,\kappa) \rho_{\mathfrak{b},\mathfrak{c}}(-n,\kappa) d\kappa,\notag
\end{align}
where the Bessel transforms are defined by
\begin{align}
&\tilde{\phi}(k)=4i^k\int_0^\infty\phi(x)J_{k-1}(x)\frac{dx}x,\notag\\
&\hat{\phi}(\kappa)=2\pi i\int_0^\infty\phi(x)\frac{J_{2i\kappa}(x)-J_{-2i\kappa}(x)}{\sinh(\pi \kappa)}\frac{dx}x,\notag\\
&\breve{\phi}(\kappa)=8\int_0^\infty\phi(x)\cosh(\pi\kappa)K_{2i\kappa}(x)\frac{dx}x.\notag
\end{align}
\end{lemma}

The Kuznetsov formula is often used together with the spectral large sieve inequalities, proved by Deshouillers and Iwaniec \cite[Theorem 2]{DI82}.
\begin{lemma}[Spectral large sieve]\label{lemsls}
Let $K\ge 1$, $N\ge1$, and $(a_n)$ be a sequence of complex numbers. Let $\mathfrak{a}$ be a cusp of $\Gamma_0(Q)$, which is equivalent to some $\frac u w$ with positive coprime $u$ and $w$ such that $w\mid Q$. Then all three quantities
\[
\sum_{\substack{2\le k\le K\\ k~ \text{even}}}\Gamma(k)\sum_{f\in\mathcal{B}_k(Q)}\left|\sum_{n}a(n)\sqrt{n}\rho_f(\mathfrak{a},n)\right|^2,  \ \ \ \ \ \ \ \ \ \ \ \ \sum_{\substack{f\in\mathcal{B}(Q)\\ |\kappa_f|\le K}}\frac1{\cosh(\pi\kappa_f)}\left|\sum_n a_n\sqrt{n}\rho_f(\mathfrak{a},\pm n)\right|^2,
\]
\[
\sum_{\mathfrak{c}}\int_{-K}^K\frac1{\cosh(\pi\kappa)}\left|\sum_n a_n\sqrt{n}\rho_{\mathfrak{a},\mathfrak{c}}(\pm n,\kappa)\right|^2d\kappa
\]
are bounded by
\[
(K^2+\muup(\mathfrak{a})N^{1+\ve})\sum_n|a_n|^2,
\]
where $\muup(\mathfrak{a})$ is defined as
\[
\muup(\infty)=Q^{-1}\ \ \text{and} \ \ \muup(\mathfrak{a})=\left(w,\frac Qw\right)Q^{-1}.
\]
\end{lemma}









\section{Sketch of the proof}\label{secT}
Let us recall that
\begin{align}
M(\alpha,\beta,\gamma,&\delta)=\frac2{\vp^*(q)}\mathop{{\sum}^+}_{\chi(\bmod q)}L\left(\tfrac12+\alpha,\chi\right) L\left(\tfrac12+\beta,\chi\right) L\left(\tfrac12+\gamma,\ol{\chi}\right) L\left(\tfrac12+\delta,\ol{\chi}\right).\notag
\end{align}
By applying the approximate function equation in Lemma \ref{lemafe} and the orthogonality formula in Lemma \ref{lemof}, we split $M$ as
\begin{align}
M(\alpha,\beta,\gamma,\delta)=A_{1}(\alpha,\beta,\gamma,\delta) +A_{-1}(\alpha,\beta,\gamma,\delta),\notag
\end{align}
where
\begin{align}
A_1(\alpha,\beta,\gamma,\delta)=&\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)\sum_{\substack{(mn,q)=1\\ m\equiv \pm n(\bmod d)}}\frac{\sigma_{\alpha-\beta}(m) \sigma_{\gamma-\delta}(n)} {m^{\frac12+\alpha}n^{\frac12+\gamma}}
V_{\alpha,\beta,\gamma,\delta}\left(\frac{mn}{q^2}\right)\notag
\end{align}
and
\begin{align}\label{eqrA_1A_{-1}}
A_{-1}(\alpha,\beta,\gamma,\delta)=X_{\alpha,\beta,\gamma,\delta}A_1(-\gamma,-\delta,-\alpha,-\beta).
\end{align}
Thus, it is sufficient to deal with $A_1$.

To treat $A_1$, we split it into diagonal terms and off-diagonal terms, and then classify the off-diagonal terms into two categories according to $m\equiv n \pmod d$ and $m\equiv -n\pmod d$. Furthermore, we split $m\equiv n \pmod d$ according to $m<n$ and $m>n$. Consequently, we rewrite $A_{1}$ as
\begin{align}\label{eqA_1}
A_{1}(\alpha,\beta,\gamma,\delta)=A_{D}(\alpha,\beta,\gamma,\delta) +A_{O}(\alpha,\beta,\gamma,\delta)+A^*_{O}(\alpha,\beta,\gamma,\delta) +A_{\ol{O}}(\alpha,\beta,\gamma,\delta),
\end{align}
where $A_{D}$ and $A_{\ol{O}}$ denote respective contributions from the diagonal terms and the off-diagonal terms with $m\equiv -n \pmod d$, and where $A_{O}$ and $A^*_{O}$ denote contributions from such off-diagonal terms with $m\equiv n \pmod d$, according to $m<n$ and $m>n$.

The diagonal terms mainly contribute to the main term, and we evaluate $A_D$ in the next section, which is not complicated.
For the off-diagonal terms, there is an obvious relationship
\begin{align}\label{A*OAO}
A^*_O(\alpha,\beta,\gamma,\delta)=A_O(\gamma,\delta,\alpha,\beta),
\end{align}
and thus we just need to treat $A_O$ and $A_{\ol{O}}$, which are specifically defined by
\begin{align}
A_O(\alpha,\beta,\gamma,\delta)=&\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right) \sum_{\substack{(mn,q)=1, m<n\\ m\equiv  n(\bmod d)}}\frac{\sigma_{\alpha-\beta}(m) \sigma_{\gamma-\delta}(n)} {m^{\frac12+\alpha}n^{\frac12+\gamma}}
V_{\alpha,\beta,\gamma,\delta}\left(\frac{mn}{q^2}\right)\notag
\end{align}
and
\begin{align}
A_{\ol{O}}(\alpha,\beta,\gamma,\delta)=&\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right) \sum_{\substack{(mn,q)=1\\ m\equiv  -n(\bmod d)}}\frac{\sigma_{\alpha-\beta}(m) \sigma_{\gamma-\delta}(n)} {m^{\frac12+\alpha}n^{\frac12+\gamma}}
V_{\alpha,\beta,\gamma,\delta}\left(\frac{mn}{q^2}\right).\notag
\end{align}
Applying the dyadic partition of unity to the sums over $m, n$ gives
\begin{align}\label{eqAOBMN}
A_{O}=\sum_{M, N} B_{M,N}
\end{align}
and
\begin{align}\label{eqAolOBMN}
A_{\ol{O}}=\sum_{MN} B_{\ol{M,N}}=\sum_{M\le N} B_{\ol{M,N}}+\sum_{M>N} B_{\ol{M,N}},
\end{align}
where $B_{M,N}$ and $B_{\ol{M,N}}$ are the same expressions as $A_{O}$ and $A_{\ol{O}}$ respectively, but weighted by $W\left(\frac mM\right) W\left(\frac nN\right)$.

By swapping $M$ and $N$, one notes that
\[
B_{\ol{M,N}}(\alpha,\beta,\gamma,\delta)=B_{\ol{N,N}}(\gamma,\delta,\alpha,\beta).
\]
Thus we can treat the $M>N$ case by transforming it to the $M\le N$ case. Doing so allows for the treatment of the error term with $M\le N$ only, but there is one thing very important to note: we must insure the same expression for the main term of the two cases before restoring the partition.


Since we focus on $B_{\ol{M,N}}$ with $M\le N$ and only terms with $m<n$ are counted in $B_{M,N}$, we may assume
\[
M\ll Nq^\ve,
\]
a convention that holds throughout the work. We may also assume $MN\le q^{2+\ve}$ due to the rapid decay of $V$.

We primarily focus on the treatment of $B_{M,N}$ since the one for $B_{\ol{M,N}}$ is much the same with some variations, and we will point out their differences especially.
The treatment of off-diagonal terms in $B_{M,N}$ is extremely complicated, which takes up most of the space.
We classify these off-diagonal terms into two categories, corresponding to $M, N$ close and far away from each other.
As we will see, the first category contributes to the main term, while the second category contributes to the error only.

To evaluate $B_{M,N}$, we write it as
\begin{align}
B_{M,N}=(\text{Main term})_{M,N}+E_{M,N}
\end{align}
for a certain main term that is too complex to explicitly write here.
We deduce an exact expression for $(\text{Main term})_{M,N}$ in Section \ref{secmainAO} \footnote{Young have chosen to directly deduce the summation of main terms over $M,N$ instead, because there is an obstruction in \cite{You11} to deduce an exact expressions for a single $(\text{Main term})_{M,N}$.}. Deduced from the exact expression, the following upper bound
\begin{align}\label{eqbmain}
(\text{Main term})_{M,N}\ll M^{\frac12}N^{-\frac12}q^\ve
\end{align}
holds, which turns out to be very small when $M,N$ are not close enough.

To bound the error term, we prove the following estimate for $E_{M,N}$ in Section \ref{secError}.
\begin{theorem}\label{thmEMN}
For $M\ll N$, $MN\ll q^{2+\ve}$, we have
\begin{align}\label{eqEMN}
E_{M,N}\ll q^{-\frac12+\theta+\ve}M^{-\frac12}N^{\frac12}.
\end{align}
\end{theorem}


The bounds \eqref{eqbmain} and \eqref{eqEMN} also holds for $B_{\ol{M,N}}$, and the treatment proceeds much the same as $B_{M,N}$, but with some variations, which are mainly because the different Gamma factors, caused by $m\equiv-n$.


In Section \ref{secmain}, we combine all various main terms from $A_O$, $A_{\ol{O}}$, etc., to deduce the main term of Theorem \ref{thmmains}. In Section \ref{secfar}, we complete our proof of Theorem \ref{thmmains} by bounding off-diagonal terms with $M,N$ far away from each other, where $B_{M,N}$ and $B_{\ol{M,N}}$ are treated identically.  We will also need the trivial bound
\begin{align}\label{eqtrivalbound}
B_{M,N}, B_{\ol{M,N}}\ll q^{-1+\ve}(MN)^{\frac12}.
\end{align}



\section{Diagonal terms}
For the diagonal terms with $m=n$, we have
\begin{align}
A_D(\alpha,\beta,\gamma,\delta)=&\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)\sum_{(n,q)=1}\frac{\sigma_{\alpha-\beta}(n) \sigma_{\gamma-\delta}(n)} {n^{1+\alpha+\gamma}}
V_{\alpha,\beta,\gamma,\delta}\left(\frac{n^2}{q^2}\right)\notag\\
=&\frac1{2\pi i}\int_{(1)}\frac{G(s)}{s}g_{\alpha,\beta,\gamma,\delta}(s)q^{2s} \sum_{(n,q)=1}\frac{\sigma_{\alpha-\beta}(n) \sigma_{\gamma-\delta}(n)} {n^{1+\alpha+\gamma+2s}}ds.\notag
\end{align}
By the Ramanujan identity, the sum over $n$ is
\begin{align}
\frac{\zeta_q(1+\alpha+\gamma+2s) \zeta_q(1+\alpha+\delta+2s)\zeta_q(1+\beta+\gamma+2s)\zeta_q(1+\beta+\delta+2s)} {\zeta_q(2+\alpha+\beta+\gamma+\delta+4s)},\notag
\end{align}
which has simple poles at $2s=-\alpha-\gamma$, etc., while $G(s)$ vanishes at these poles. We move $\text{Re}(s)$ to $-\frac14+\ve$, passing the pole at $s=0$ only. Then the integral along the new line is
bounded by
$\ll q^{-\frac12+\ve}$,
and the pole at $s=0$ gives
$Z_q(\alpha,\beta,\gamma,\delta)$.
We summarize this calculation in the following:
\begin{lemma}
\label{lemD}
We have
\begin{align}
A_D(\alpha,\beta,\gamma,\delta)=Z_q(\alpha,\beta,\gamma,\delta) +O\left(q^{-\frac12+\ve}\right),\notag
\end{align}
and similarly the contribution of the diagonal terms to $A_{-1}$ is
\begin{align}
A_{-D}(\alpha,\beta,\gamma,\delta)=X_{\alpha,\beta,\gamma,\delta} Z_q(-\gamma,-\delta,-\alpha,-\beta) +O\left(q^{-\frac12+\ve}\right).\notag
\end{align}
\end{lemma}



\section{The $\mathcal{D}_{q}$-function}\label{secDq}
For off-diagonal terms, we need to treat the divisor sum of type
\[
\sum_{n\in S_q}\frac{\sigma_\lambda(n)} {n^s} e\left(n\frac h{l}\right)
\]
with $S_q$ being a set of positive integers meeting some coprime conditions. More precisely, we should deduce the theory for this type divisor sum with $S_q=\{n\in \mathbb{N}: (n,q)=1\}$ and $S_q=\{n\in \mathbb{N}: (n,q)=1, (n+r,q)=1\}$ for any given integer $r$.
\subsection{The generalized Estermann $D$-function with a character}
Suppose that $q, l, h$ are integers, and $\chi$ is a character modulo $q$.  For any given $\lambda\in\mathbb{C}$, the generalized Estermann $D$-function with a character is defined as
\begin{align}
D\left(s,\lambda,\frac h{lq}, \chi\right)=\sum_{n}\frac{\sigma_\lambda(n)\chi(n)}{n^s}e\left(n\frac h{lq}\right).
\end{align}
\begin{lemma}\label{lemGED}
For any fixed $\lambda\in\mathbb{C}$, $D(s,\lambda,\frac h{lq}, \chi)$ is meromorphic as a function of $s$, satisfying the functional equation
\begin{align}\label{eqfeqD}
D\left(s,\lambda,\frac h{lq},\chi\right)&=2(2\pi)^{-2-\lambda+2s}(lq)^{\lambda-2s}\Gamma\left(1-s\right) \Gamma\left(1+\lambda-s\right)\\
&\times\left[-\cos\left(\frac{\pi}2\left(2s-\lambda\right)\right) \mathcal{A}_1\left(1-s,\lambda,\frac h{lq},\chi\right)+\cos\left(\frac{\pi\lambda}2\right)\mathcal{A}_2\left(1-s,\lambda,\frac h{lq},\chi\right)\right],\notag
\end{align}
where
\begin{align}\label{eqA1-2}
\mathcal{A}_1\left(s,\lambda,\frac h{lq},\chi\right)&=\sum_{1\le u,v\le lq}\chi(uv)e\left(\frac{uvh}{lq}\right)F\left(s,\frac{u}{lq}\right) F\left(s+\lambda,\frac{v}{lq}\right),\\
\mathcal{A}_2\left(s,\lambda,\frac h{lq},\chi\right)&=\sum_{1\le u,v\le lq}\chi(uv)e\left(\frac{uvh}{lq}\right)F\left(s,\frac{u}{lq}\right) F\left(s+\lambda,-\frac{v}{lq}\right).
\end{align}
If $\lambda\neq 0$, then $D$ has simple poles at $s=1$ and $s=1+\lambda$ with respective residues
\begin{align}
&(lq)^{-2+\lambda}\sum_{1\le u,v\le lq}\chi(uv)e\left(\frac{uvh}{lq}\right) \zeta\left(1-\lambda,\frac{v}{lq}\right), \notag\\
&(lq)^{-2-\lambda}\sum_{1\le u,v\le lq}\chi(uv)e\left(\frac{uvh}{lq}\right)\zeta\left(1+\lambda,\frac{u}{lq}\right).\notag
\end{align}
\end{lemma}
\begin{proof}
The functional equation has already been obtained by Wu \cite[Lemma 4.3]{Wu19}, whose focus is primitive character. However, the primitive character condition is unused in the proof of the functional equation, and so it applies here.

To the residues, we observe that
\begin{align}\label{eqD}
D\left(s,\lambda,\frac h{lq},\chi\right)=(lq)^{-2s+\lambda}\sum_{1\le u,v\le lq}\chi(uv)e\left(\frac{uvh}{lq}\right)\zeta\left(s,\frac{u}{lq}\right) \zeta\left(s-\lambda,\frac{v}{lq}\right).
\end{align}
Then calculating the residue of the Hurwitz zeta-function at the pole $1$ provides the two residues for the lemma.
\end{proof}

Note that, in the functional equation $\eqref{eqfeqD}$, we do not need the coprime condition $(h,lq)=1$. For the special case with $(h,lq)=1$ and $\chi=\chi_0$, we may simplify $\mathcal{A}_i$ by the following lemma.

\begin{lemma}\label{lemiD}
Let $h,l,q$ be three integers satisfying $(h,lq)=1$ and $\chi_0$ be the principal character modulo $q$. For any given $\lambda\in\mathbb{C}$, we have
\begin{align}
D\left(s,-\lambda,-\frac {\ol{h}}{lq}\right)=\frac1{lq}\sum_{1\le uv\le lq}e\left(\frac{uvh}{lq}\right)\sum_{m}\frac1{m^s}e\left(\frac{mu}{lq}\right) \sum_{n}\frac1{n^{s+\lambda}}e\left(\frac{nv}{lq}\right)\notag
\end{align}
and
\begin{align}
D\left(s,-\lambda,-\frac {\ol{h}}{lq},\chi_0\right)=\frac1{lq}\sum_{1\le uv\le lq}e\left(\frac{uvh}{lq}\right)\sum_{(m,q)=1}\frac1{m^s}e\left(\frac{mu}{lq}\right) \sum_{(n,q)=1}\frac1{n^{s+\lambda}}e\left(\frac{nv}{lq}\right).\notag
\end{align}
\end{lemma}
\begin{proof}
We show the second identity in full details, and the treatment of the first identity is identical.
For the sums on the right-hand side of the second identity, we sum over $v$ first, which is
\begin{align}
\sum_{1\le v\le lq}e\left(\frac{(n+uh)v}{lq}\right)=
\left\{
\begin{aligned}
&lq\ \ \ \ \text{for} \ \ \ lq\mid n+uh,\\
&0\ \ \ \ \text{otherwise}.
\end{aligned}
\right.\notag
\end{align}
Thus, the right-hand side evolves into
\begin{align}
\sum_{1\le u\le lq}\sum_{(m,q)=1}\frac1{m^s}e\left(\frac{mu}{lq}\right) \sum_{\substack{n\equiv -uh (\bmod lq)\\ (n,q)=1}} \frac1{n^{s+\lambda}} =\sum_{(mn,q)=1}\frac{1}{m^sn^{s+\lambda}}e\left(-\frac{mn\ol{h}}{lq}\right)\notag
\end{align}
by an exchange of the summations, and this establishes the lemma.
\end{proof}

Since coprime conditions in the sum do not rely on higher power factors, we may restrict ourselves to square-free $q$.
\begin{proposition}\label{progED1}
Let $q$ be a square-free number, and $\chi_0$ be the principal character modulo $q$. For $q, l, h$ coprime with each other and $\lambda\in\mathbb{C}$, $D(s,\lambda,\frac h{lq}, \chi_0)$ satisfies the functional equation
\begin{align}
D\left(s,\lambda,\frac h{lq},\chi_0\right)&=2(2\pi)^{-2-\lambda+2s}\Gamma\left(1-s\right) \Gamma\left(1+\lambda-s\right)
\sum_{q_1q_2q_3=q}\frac{\mu(q_2q_3)}{q_2q_3}(lq_1q_2)^{1+\lambda-2s}\mathop{\sum\nolimits^*}_{r (\bmod q_2)}\notag\\
&\times\left[-\cos\left(\frac{\pi}2\left(2s-\lambda\right)\right) D\left(1-s,-\lambda,-\frac{\ol{x}_r}{lq_1q_2}\right) +\cos\left(\frac{\pi\lambda}2\right)D\left(1-s,-\lambda,\frac{\ol{x}_r}{lq_1q_2}\right)\right]\notag
\end{align}
with
\begin{align}
\left\{
\begin{aligned}
&\ol{x}_r\equiv \ol{r} \pmod{q_2},\\
&\ol{x}_r\equiv \ol{h}q_3 \pmod{lq_1}.
\end{aligned}
\right.\notag
\end{align}
In addition, the residues at $s=1$ and $s=1+\lambda$ are respectively
\begin{align}
\mu(q)l^{-1+\lambda}q^{-1}\zeta_q(1-\lambda), \ \ \ \ \ \ \ \mu(q)l^{-1-\lambda}q^{-1}\zeta_q(1+\lambda).\notag
\end{align}
\end{proposition}
\begin{proof}
By Lemma \ref{lemGED}, the functional equation may follow from the evolutions of $\mathcal{A}_i$ for $i=1,2$. Let us recall that
\begin{align}
\mathcal{A}_1\left(s,\lambda,\frac h{lq},\chi_0\right)=\sum_{\substack{1\le u,v\le lq\\ (uv,q)=1}}e\left(\frac{uvh}{lq}\right)F\left(s,\frac{u}{lq}\right) F\left(s+\lambda,\frac{v}{lq}\right).\notag
\end{align}
We remove the coprime condition with \eqref{eqcoprimee} to get
\begin{align}\label{eqA1ab}
\mathcal{A}_1=\sum_{q_0q_1=q}\frac{\mu(q_0)}{q_0}\sum_{r(\bmod q_0)}\sum_{1\le u,v\le lq}e\left(\frac{uv(h+rlq_1)}{lq}\right)\sum_{m}\frac1{m^{s}}e\left(\frac{mu}{lq}\right) \sum_{n}\frac1{n^{s+\lambda}}e\left(\frac{nv}{lq}\right).
\end{align}
Since $h+rlq_1$ just runs through the complete system of residue classes modulo $q_0$ as $r$ does, there is
\begin{align}
\sum_{r (\bmod q_0)}e\left(\frac{uv(h+rlq_1)}{lq}\right)=\sum_{r (\bmod q_0)}e\left(\frac{uvh_r}{lq}\right)\notag
\end{align}
with
\begin{align}
\left\{
\begin{aligned}
&h_r\equiv r,\pmod{q_0},\\
&h_r\equiv h,\pmod{lq_1}.\notag
\end{aligned}
\right.
\end{align}
For $(h,q)=1$, it is easy to see that
\[
(h_r,q)=(h_r,q_0)=(r,q_0).
\]
Let $q_0=q_2q_3$ with $q_3=(h_r,q_0)$. We reduce $q_3$ in the exponential function to get
\begin{align}
\sum_{r (\bmod q_0)}e\left(\frac{uvh_r}{l q}\right)=\sum_{q_2q_3=q_0}\mathop{\sum\nolimits^*}_{r(\bmod q_2)}e\left(\frac{uvx_r}{lq_1q_2}\right)\notag
\end{align}
with
\begin{align}
\left\{
\begin{aligned}
&x_r\equiv r \pmod{q_2},\\
&x_r\equiv h\ol{q}_3 \pmod{lq_1}.
\end{aligned}
\right.\notag
\end{align}
Applying this into \eqref{eqA1ab} gives
\begin{align}
\mathcal{A}_1=&\sum_{q_1q_2q_3=q}\frac{\mu(q_2q_3)}{q_2q_3}\mathop{\sum\nolimits^*}_{r (\bmod q_2)}\sum_{1\le u,v\le lq}e\left(\frac{uvx_r}{lq_1q_2}\right) \sum_{m}\frac1{m^{s}}e\left(\frac{mu}{lq}\right) \sum_{n}\frac1{n^{s+\lambda}}e\left(\frac{nv}{lq}\right).\notag
\end{align}

We observe that contribution only comes from such terms with $q_3\mid m$ and $q_3\mid n$, while others vanish if we sum over $u,v$ first. To see this, we sum over $u$ first and find that
\begin{align}
&\sum_{1\le u,v\le lq}e\left(\frac{uvx_r}{lq_1q_2}\right)\sum_{m}\frac1{m^{s}}e\left(\frac{mu}{lq}\right) \sum_{n}\frac1{n^{s+\lambda}}e\left(\frac{nv}{lq}\right)\notag\\
=&\sum_{1\le v\le lq}\sum_n\frac1{n^{s+\lambda}}e\left(\frac{nv}{lq}\right) \sum_{m}\frac1{m^{s}}\sum_{1\le u\le lq}e\left(\frac{vx_r q_3+m}{lq}u\right),\notag
\end{align}
where
\begin{align}
\sum_{1\le u\le lq}e\left(\frac{vx_r q_3+m}{lq}u\right)=\left\{
\begin{aligned}
&lq \ \ \ \text{for} \ \ \ vx_r q_3\equiv -m \pmod{lq},\\
&0\ \ \ \ \text{otherwise}.\notag
\end{aligned}
\right.
\end{align}
This implies that contribution arises just from such terms with
\begin{align}
-m\equiv vx_r q_3~ (\bmod lq),\notag
\end{align}
in particular $q_3\mid m$. The same phenomenon happens for $n$.
Thus, we may add the conditions
\[
q_3\mid m, \ \ \ q_3\mid n
\]
to the sum at no cost.

We extract $q_3$ from the sum to have
\begin{align}
\mathcal{A}_1=&\sum_{q_1q_2q_3=q}\frac{\mu(q_2q_3)}{q_2q_3}q_3^{2-2s-\lambda}\mathop{\sum\nolimits^*}_{r (\bmod q_2)} \sum_{1\le u,v\le lq_1q_2}e\left(\frac{uvx_r}{lq_1q_2}\right) \sum_{m}\frac1{m^{s}}e\left(\frac{mu}{lq_1q_2}\right) \sum_{n}\frac1{n^{s+\lambda}}e\left(\frac{nv}{lq_1q_2}\right).\notag
\end{align}
Then it follows from the first identity of Lemma \ref{lemiD} that
\begin{align}\label{eqA1D}
\mathcal{A}_1\left(s,\lambda,\frac h{lq},\chi_0\right)=\sum_{q_1q_2q_3=q}\frac{\mu(q_2q_3)}{q_2q_3}q_3^{2-2s-\lambda}lq_1q_2\mathop{\sum\nolimits^*}_{r (\bmod q_2)}
D\left(s,-\lambda,-\frac{\ol{x_r}}{lq_1q_2}\right).
\end{align}
Similarly, there is
\begin{align}\label{eqA2D}
\mathcal{A}_2\left(s,\lambda,\frac h{lq},\chi_0\right)=\sum_{q_1q_2q_3=q}\frac{\mu(q_2q_3)}{q_2q_3}q_3^{2-2s-\lambda}lq_1q_2\mathop{\sum\nolimits^*}_{r (\bmod q_2)}
D\left(s,-\lambda,\frac{\ol{x_r}}{lq_1q_2}\right).
\end{align}
Thus, the functional equation follows immediately from \eqref{eqfeqD} with these new expressions for $\mathcal{A}_i$.


We now come to the residues. By Lemma \ref{lemGED}, the residue at $1$ is
\begin{align}
res(1)=(lq)^{-2+\lambda}\sum_{1\le u,v\le lq}\chi_0(uv)e\left(\frac{uvh}{lq}\right) \zeta\left(1-\lambda,\frac{v}{lq}\right).\notag
\end{align}
We first sum over $u$ by writing it as $u=jq+r$ to get
\begin{align}
res(1)=(lq)^{-2+\lambda}\sum_{\substack{1\le v\le lq\\ (v,q)=1}}\zeta\left(1-\lambda,\frac{v}{lq}\right)\mathop{\sum\nolimits^*}_{r (\bmod q)} e\left(\frac{rvh}{lq}\right)\sum_{j(\bmod l)} e\left(\frac{jvh}{l}\right).\notag
\end{align}
The sum over $j$ is $0$ unless $l\mid v$, in which case it equals $l$. Making the change of variables $v\rightarrow vl$, we have
\begin{align}
res(1)&=l^{-1+\lambda}q^{-2+\lambda}\mathop{\sum\nolimits^*}_{v(\bmod q)} c_q(vh) \zeta\left(1-\lambda,\frac vq\right)\notag\\
&=\mu(q)l^{-1+\lambda}q^{-2+\lambda}\mathop{\sum\nolimits^*}_{v(\bmod q)} \zeta\left(1-\lambda,\frac vq\right)\notag\\
&=\mu(q)l^{-1+\lambda}q^{-1}\zeta_q(1-\lambda).\notag
\end{align}
The residue at $1+\lambda$ follows similarly.
\end{proof}



\subsection{The $\mathcal{D}_{q}$-function}
For integers $q, l, r, h$ and any $\lambda\in\mathbb{C}$, we define the $\mathcal{D}_{q}$-function as
\begin{align}
\mathcal{D}_{q}\left(s,\lambda,\frac h{l},r\right)=\sum_{\substack{(n,q)=1\\ (n+r,q)=1}}\frac{\sigma_\lambda(n)} {n^s} e\left(n\frac h{l}\right).\notag
\end{align}
Obviously, there is
\[
\mathcal{D}_{q}\left(s,\lambda,\frac h{l},r\right)=\mathcal{D}_{q^*}\left(s,\lambda,\frac h{l},r\right).
\]


\begin{proposition}\label{lemgED2}
Let $q, l, r, h$ be integers, satisfying $(l,qh)=1$. For any fixed $\lambda\in\mathbb{C}$, $\mathcal{D}_{q}(s,\lambda,\frac h{l},r)$ is meromorphic as a function of $s$, satisfying the functional equation
\begin{align}
\mathcal{D}_{q}&\left(\tfrac12+s,\lambda,\frac h{l},r\right)=2(2\pi)^{-1-\lambda+2s}\Gamma\left(\tfrac12-s\right) \Gamma\left(\tfrac12+\lambda-s\right)\notag\\
&\ \ \ \ \ \ \ \ \ \ \times\sum_{a\mid q}\sum_{b\mid q_r}\frac{\mu(ab)}{ab}\sum_{a_1\mid a}\sum_{b_1\mid b} (la_1b_1)^{\lambda-2s} \mathop{\sum\nolimits^*}_{i(\bmod a_1)}\mathop{\sum\nolimits^*}_{j(\bmod b_1)}
e\left(\frac{jr}{b_1}\right) \notag\\
&\ \ \ \ \ \ \ \ \ \ \times\left[\sin\left(\pi\left(s-\frac\lambda2\right)\right) D\left(\tfrac12-s,-\lambda,-\frac{\ol{h_{i,j}}}{la_1b_1},\chi'_0\right) +\cos\left(\frac{\pi\lambda}2\right) D\left(\tfrac12-s,-\lambda,\frac{\ol{h_{i,j}}}{la_1b_1},\chi'_0\right)\right]\notag
\end{align}
with $\chi'_0$ being the principal character modulo $a_1b_1$ and
\begin{align}
\left\{
\begin{aligned}
&h_{i,j}\equiv ilb_1  \pmod{a_1},\\
&h_{i,j}\equiv jla_1 \pmod{b_1},\\
&h_{i,j}\equiv ha_1b_1 \pmod{l}.
\end{aligned}
\right.\notag
\end{align}
If $\lambda\neq 0$, then $\mathcal{D}_{q}$ has simple poles at $s=1$ and $s=1+\lambda$ with respective residues
\begin{align}
\frac{\vp(q)}{q}\sum_{b\mid q_r}\frac{\mu(b)}{\vp(b)} l^{-1+\lambda}\zeta_q(1-\lambda),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  \frac{\vp(q)}{q}\sum_{b\mid q_r}\frac{\mu(b)}{\vp(b)} l^{-1-\lambda}\zeta_q(1+\lambda).\notag
\end{align}
\end{proposition}




\begin{proof}
Note that a replacement of $q$ with $q^*$ does not affect the functional equation, as well as the two residues. Thus, it is sufficient to prove the proposition for square-free $q$, which always holds in our following proof.
We first express the coprime conditions in the sum of the $\mathcal{D}_q$-function in terms of multiplicative characters with \eqref{eqcoprimechi}, then
\begin{align}\label{eqDsharpD}
\mathcal{D}_{q}\left(s,\lambda,\frac h{l},r\right)=\sum_{b\mid q}\frac{\mu(b)}{\vp(b)} \sum_{\chi(\bmod b)}\ol{\chi}(-r)D\left(s,\lambda,\frac {hq}{lq},\chi_0\chi\right),
\end{align}
where $\chi_0$ is the principal character modulo $q$.
We apply the functional equation \eqref{eqfeqD} to the generalized Estermann $D$-function. Then it follows
\begin{align}\label{eqD_q+1}
\mathcal{D}_{q}\left(s,\lambda,\frac h{l},r\right)=&2(2\pi)^{-2-\lambda+2s}(lq)^{\lambda-2s} \Gamma(1-s)\Gamma(1+\lambda-s)\\
&\times\left[-\cos\left(\frac{\pi}2(2s-\lambda)J_1(1-s) +\cos\left(\frac{\pi\lambda}{2}\right)J_2(1-s)\right)\right]\notag
\end{align}
with
\begin{align}
J_i(s)=\sum_{b\mid q}\frac{\mu(b)}{\vp(b)}\sum_{\chi (\bmod b)}\ol{\chi}(-r)\mathcal{A}_i\left(s,\lambda,\frac {hq}{lq},\chi_0\chi\right)\notag
\end{align}
for $i=1,2$.
Applying the expression \eqref{eqA1-2} and then executing the sums over $b$ and $\chi$ with \eqref{eqcoprimechi}, we have
\begin{align}\label{eqchiJ1}
J_1(s)=\sum_{\substack{1\le u,v\le lq\\ (uv,q)=1\\ (uv+r,q)=1}}e\left(\frac{uvh}{l}\right) F\left(s,\frac{u}{lq}\right) F\left(s+\lambda,\frac{v}{lq}\right).
\end{align}

The coprime conditions on $uv$ are obstructions to calculate the sum, and we should remove them again.
With $q=(q,r)q_{_r}$, it is obvious that
\begin{align}
\begin{array}{c}
(uv,q)=1\\
(uv+r,q)=1
\end{array}\Longleftrightarrow
\begin{array}{c}
(uv,q)=1\\
(uv+r,q_r)=1
\end{array}.\notag
\end{align}
Applying this and removing the coprime conditions by M\"{o}bius inversion, we have
\begin{align}
\sum_{\substack{1\le u,v\le lq\\ (uv,q)=1\\ (uv+r,q)=1}}=\sum_{a\mid q}\mu(a)\sum_{b\mid q_r}\mu(b)\sum_{\substack{1\le u,v\le lq\\ a\mid uv\\ b\mid uv+r}}=\sum_{a\mid q}\sum_{b\mid q_r}\mu(ab)\sum_{\substack{1\le u,v\le lq\\ a\mid uv\\ b\mid uv+r}}\notag
\end{align}
by assuming $(a,b)=1$ in the sum at no cost.

Besides, we apply \eqref{eqdivisor} to write the conditions $a\mid uv$ and $b\mid uv+r$ in terms of the Ramanujan sum, and then it follows
\begin{align}
\sum_{\substack{1\le u,v\le lq\\ (uv,q)=1\\ (uv+r,q)=1}}=\sum_{a\mid q}\sum_{b\mid q_r}\frac{\mu(ab)}{ab}\sum_{a_1\mid a}\sum_{b_1\mid b}\sum_{1\le u,v\le lq}c_{a_1}(uv)c_{b_1}(uv+r).\notag
\end{align}
Inserting this into \eqref{eqchiJ1} with the exponential sum formula for the Ramanujan sum, we find
\begin{align}
J_1(s)=&\sum_{a\mid q}\sum_{b\mid q_r}\frac{\mu(ab)}{ab}\sum_{a_1\mid a}\sum_{b_1\mid b} \mathop{\sum\nolimits^*}_{i(\bmod a_1)}\mathop{\sum\nolimits^*}_{j(\bmod b_1)}
e\left(\frac{jr}{b_1}\right)\notag\\
&\times\sum_{1\le u,v\le lq} e\left(\frac{uvh}{l}\right)e\left(\frac{iuv}{a_1}\right) e\left(\frac{juv}{b_1}\right) F\left(s,\frac{u}{lq}\right) F\left(s+\lambda,\frac{v}{lq}\right).\notag
\end{align}
Then a simple arrangement shows
\begin{align}\label{eqJ1lq}
J_1(s)=&\sum_{a\mid q}\sum_{b\mid q_r}\frac{\mu(ab)}{ab}\sum_{a_1\mid a}\sum_{b_1\mid b} \mathop{\sum\nolimits^*}_{i(\bmod a_1)}\mathop{\sum\nolimits^*}_{j(\bmod b_1)}
e\left(\frac{jr}{b_1}\right)\\
&\times\sum_{1\le u,v\le lq} e\left(\frac{uv(ha_1b_1+jla_1+ilb_1)}{la_1b_1}\right) \sum_{m}\frac1{m^{s}}e\left(\frac{mu}{lq}\right) \sum_{n}\frac1{n^{s+\lambda}}e\left(\frac{nv}{lq}\right).\notag
\end{align}

In the expression \eqref{eqJ1lq}, we should note that there are many terms in the sum over $u,v$, which make no essential contribution to $J_1$. In particular, all terms with $(u,q)>1$ or $(v,q)>1$ are in the case, which certainly vanish if we sum over $a,b,a_1,b_1,i,j$ first. What's more, we note that many terms in the sum over $m,n$ vanish if we sum over other variables first. As we will see, these nonessential terms do not contribute to $\mathcal{A}_1$, but lead to obstructions to bound the error when we apply the functional equation. Thus, we should kick them out now. More precisely, we observe that all terms in the sum over $m,n$ satisfying at least one of the following three conditions
\begin{align}
(mn,a_1b_1)>1,\ \ \ q_{a_1b_1}\nmid m,\ \ \ q_{a_1b_1}\nmid n\notag
\end{align}
can be removed at no cost. Let us see this in detail.
For notational convenience, we apply the following notation
\begin{align}
h_{i,j}=ha_1b_1+jla_1+ilb_1.\notag
\end{align}
Then it is easy to see that $(h_{i,j},la_1b_1)=1$
and
\begin{align}
\left\{
\begin{aligned}
&h_{i,j}\equiv ilb_1  \pmod{a_1},\\
&h_{i,j}\equiv jla_1 \pmod{b_1},\\
&h_{i,j}\equiv ha_1b_1 \pmod{l}.
\end{aligned}
\right.\notag
\end{align}
We show in full detail for the sum over $m$, and the case for $n$ is identical.
More precisely, we first kick out all nonessential terms with $(v,q)>1$ at no cost, and then we sum over $u$ to find
\begin{align}
&\sum_{\substack{1\le u,v\le lq\\ (v,q)=1}}e\left(\frac{uvh_{i,j}}{la_1b_1}\right)\sum_{m}\frac1{m^{s}}e\left(\frac{mu}{lq}\right) \sum_{n}\frac1{n^{s+\lambda}}e\left(\frac{nv}{lq}\right)\\
=&\sum_{\substack{1\le v\le lq\\ (v,q)=1}}\sum_n\frac1{n^{s+\lambda}}e\left(\frac{nv}{lq}\right) \sum_{m}\frac1{m^{s}}\sum_{1\le u\le lq}e\left(\frac{vh_{i,j}q_{a_1b_1}+m}{lq}u\right),\notag
\end{align}
where
\begin{align}
\sum_{1\le u\le lq}e\left(\frac{vh_{i,j}q_{a_1b_1}+m}{lq}u\right)=\left\{
\begin{aligned}
&lq \ \ \ \text{for} \ \ \ vh_{i,j}q_{a_1b_1}\equiv -m \pmod{lq},\\
&0\ \ \ \ \text{otherwise}.\notag
\end{aligned}
\right.
\end{align}
This provides that contribution arises just from such terms that
\begin{align}
-m\equiv vh_{i,j}q_{a_1b_1}~ (\bmod lq),\notag
\end{align}
in particular $q_{a_1b_1}\mid m$ and $(m,a_1b_1)=1$ since $(vh_{i,j}q_{a_1b_1},a_1b_1)=1$. Identical discussions show the same result for the sum over $n$.
Thus, we add conditions
\[
(m,a_1b_1)=1,\ \ \ q_{a_1b_1}\mid m, \ \ \ (n,a_1b_1)=1,\ \ \ q_{a_1b_1}\mid n
\]
to the sum over $m,n$ at no cost.

We extract $q_{a_1b_2}$ from the sum to get
\begin{align}
J_1(s)=&\sum_{a\mid q}\sum_{b\mid q_r}\frac{\mu(ab)}{ab}\sum_{a_1\mid a}\sum_{b_1\mid b} (q_{a_1b_1})^{2-2s-\lambda} \mathop{\sum\nolimits^*}_{i(\bmod a_1)}\mathop{\sum\nolimits^*}_{j(\bmod b_1)}
e\left(\frac{jr}{b_1}\right)\notag\\
&\times\sum_{1\le u,v\le la_1a_2}e\left(\frac{uvh_{i,j}}{la_1b_1}\right)\sum_{(m,a_1b_1)=1}\frac1{m^{s}} e\left(\frac{mu}{la_1b_1}\right) \sum_{(n,a_1b_1)=1}\frac1{n^{s+\lambda}}e\left(\frac{nv}{la_1b_1}\right).\notag
\end{align}
Then by the second identity of Lemma \ref{lemiD},
\begin{align}\label{eqJ1}
J_1(s)=&\sum_{a\mid q}\sum_{b\mid q_r}\frac{\mu(ab)}{ab}\\
&\times\sum_{a_1\mid a}\sum_{b_1\mid b} (q_{a_1b_1})^{2-2s-\lambda} la_1b_1\mathop{\sum\nolimits^*}_{i(\bmod a_1)}\mathop{\sum\nolimits^*}_{j(\bmod b_1)}
e\left(\frac{jr}{b_1}\right)D\left(s,-\lambda,-\frac{\ol{h_{i,j}}}{la_1b_1},\chi'_0\right),\notag
\end{align}
with $\chi'_0$ being the principal character modulo $a_1b_1$.
A similar deduction shows
\begin{align}\label{eqJ2}
J_2(s)=&\sum_{a\mid q}\sum_{b\mid q_r}\frac{\mu(ab)}{ab}\sum_{a_1\mid a}\sum_{b_1\mid b} (q_{a_1b_1})^{2-2s-\lambda} la_1b_1 \mathop{\sum\nolimits^*}_{i(\bmod a_1)}\mathop{\sum\nolimits^*}_{j(\bmod b_1)}
e\left(\frac{jr}{b_1}\right)D\left(s,-\lambda,\frac{\ol{h_{i,j}}}{la_1b_1},\chi'_0\right).
\end{align}
Applying \eqref{eqJ1} and \eqref{eqJ2} into \eqref{eqD_q+1}, one obtains the functional equation after making the change of variables $s\rightarrow \frac12+s$.

To calculate the residues, we apply Lemma \ref{lemGED} into \eqref{eqDsharpD}. Then, the residue at $1$ is
\begin{align}
(lq)^{-2+\lambda}\sum_{b\mid q}\frac{\mu(b)}{\vp(b)}\sum_{\chi (\bmod b)}\ol{\chi}(-r)\sum_{1\le u,v\le lq}\chi_0(uv)\chi(uv)e\left(\frac{uvh}{l}\right) \zeta\left(1-\lambda,\frac{v}{lq}\right).\notag
\end{align}
In this formula, applying \eqref{eqcoprimechi} as before gives
\begin{align}
res(1)=(lq)^{-2+\lambda}\sum_{\substack{1\le u,v\le lq\\ (uv,q)=1\\ (uv+r,q)=1}} e\left(\frac{uvh}{l}\right) \zeta\left(1-\lambda,\frac{v}{lq}\right).\notag
\end{align}
Now we apply
\begin{align}
\begin{array}{c}
(uv,q)=1\\
(uv+r,q)=1
\end{array}\Longleftrightarrow
\begin{array}{c}
(uv,q)=1\\
(uv+r,q_r)=1\notag
\end{array}
\end{align}
and then remove the coprime conditions $(uv+r,q_r)=1$ with \eqref{eqcoprime} to find
\begin{align}\label{+eq2.67}
res(1)=&(lq)^{-2+\lambda}\sum_{b\mid q_r}\frac{\mu(b)}{b}\sum_{b_1\mid b}c_{b_1}(uv+r) \sum_{\substack{1\le u,v\le lq\\ (uv,q)=1}}e\left(\frac{uvh}{l}\right) \zeta\left(1-\lambda,\frac{v}{q}\right)\\
=&(lq)^{-2+\lambda}\sum_{b\mid q_r}\frac{\mu(b)}{b}\sum_{b_1\mid b}\mathop{\sum\nolimits^*}_{j(\bmod b_1)} e\left(\frac{jr}{b_1}\right) \sum_{\substack{1\le u,v\le lq\\ (uv,q)=1}}e\left(\frac{uv(hb_1+jl)}{lb_1}\right)  \zeta\left(1-\lambda,\frac{v}{lq}\right).\notag
\end{align}
With $u=iq+k$, the sum over $u$ is equal to
\begin{align}
\mathop{\sum\nolimits^*}_{k(\bmod q)} e\left(\frac{vk(hb_1+jl)}{lb_1}\right) \sum_{i(\bmod l)}e\left(\frac{ivq_{b_1}(hb_1+jl)}{l}\right).\notag
\end{align}
Since $(q_{b_1}(hb_1+jl),l)=1$, the last sum over $j$ vanishes unless $l\mid v$, in which case it is equal to $l$.
Applying this as well as two changes of variables $v\rightarrow lv$ and $vk\rightarrow k$ in \eqref{+eq2.67}, we have
\begin{align}
res(1)=&l(lq)^{-2+\lambda}\sum_{b\mid q_r}\frac{\mu(b)}{b}\sum_{b_1\mid b} \mathop{\sum\nolimits^*}_{v(\bmod q)} \zeta\left(1-\lambda,\frac{v}{q}\right)\mathop{\sum\nolimits^*}_{j(\bmod b_1)} e\left(\frac{jr}{b_1}\right)\mathop{\sum\nolimits^*}_{k(\bmod q)}e\left(\frac{k(hb_1+jl)}{b_1}\right). \notag
\end{align}
Since $(hb_1+jl,b_1)=1$, we apply the Chinese Remainder Theorem and write the sum over $k$ in terms of the Ramanujan sum to find
\begin{align}
\mathop{\sum\nolimits^*}_{k(\bmod q)}e\left(\frac{k(hb_1+jl)}{b_1}\right) &=\vp(q_{b_1})c_{b_1}(hb_1+jl)=\vp(q)\frac{\mu(b_1)}{\vp(b_1)}.\notag
\end{align}
Consequently, the sum over $j$ also turns out to be the Ramanujan sum, equal to $\mu(b_1)$ since $(b_1,r)=1$. Thus, we conclude that
\begin{align}
res(1)=&l(lq)^{-2+\lambda} \mathop{\sum\nolimits^*}_{v(\bmod q)} \zeta\left(1-\lambda,\frac{v}{q}\right)\left(\vp(q)\sum_{b\mid q_r}\frac{\mu(b)}{b}\sum_{b_1\mid b}\frac{\mu^2(b_1)}{\vp(b_1)}\right).\notag
\end{align}
Note that
\[
\mathop{\sum\nolimits^*}_{v(\bmod q)} \zeta\left(1-\lambda,\frac{v}{q}\right)=q^{1-\lambda}\zeta_q(1-\lambda)
\]
and
\[
\vp(q)\sum_{b\mid q_r}\frac{\mu(b)}{b}\sum_{b_1\mid b}\frac{\mu^2(b_1)}{\vp(b_1)}=\vp(q)\sum_{b\mid q_r}\frac{\mu(b)}{\vp(b)}.
\]
Therefore, we have
\begin{align}
res(1)=&\frac{\vp(q)}{q}\sum_{b\mid q_r}\frac{\mu(b)}{\vp(b)} l^{-1+\lambda}\zeta_q(1-\lambda).\notag
\end{align}
The treatment for the residue at $1+\lambda$ follows identically.
\end{proof}







\section{Initial treatment for off-diagonal terms $A_O$}\label{secnext}

\subsection{Some arithmetic sums}
We require the computation of some arithmetic sums.


\begin{lemma}\label{lemus}
For any integers $q,n$ with $(n,q)=1$, $\text{Re}(\alpha)<0$ we have
\begin{align}
\sigma_{\alpha}(n)=\zeta_{q}(1-\alpha)\sum_{(l,q)=1} \frac{c_{l}(n)}{l^{1-\alpha}}.\notag
\end{align}
\end{lemma}
\begin{proof}
This is an analog of \cite[Lemma 5.2]{You11} but with an extra coprime condition in the sum over $l$.
The proof is obvious, which follows immediately from the formula $c_l(n)=\sum_{d\mid(l,n)}d\mu(l/d)$ with a reversal of summations.
\end{proof}

\begin{lemma}\label{lemsumf}
Let $k\mid q$, $\text{Re}(s)>1$, and $\text{Re}(\lambda)>-1$. Then
\begin{align}
\sum_{(r,q)=k} \frac{1}{r^s}\sum_{(l,q)=1}\frac{c_{l}(r)} {l^{2+\lambda}} =\frac{1}{k^s}\prod_{\substack{p\mid k\\ p\nmid q/k}}\left(1-\frac1{p^s}\right)^{-1}\frac{\zeta_q(s)\zeta_q(1+\lambda+s)} {\zeta_q(2+\lambda)}.\notag
\end{align}
\end{lemma}
\begin{proof}
We apply the change of variables $r\rightarrow k^{\sharp}r$, where $k^{\sharp}$ is largest divisor of $r$ owning the same distinct prime factors as $k$. So there are
\[
(r,q)=1,\ \ \ \ (k^{\sharp}/k,q/k)=1,\ \ \ \  (k^{\sharp},l)=1,
\]
by which we have
\begin{align}
c_{l}(k^{\sharp}r)=c_l(r),\ \ \ \ \ \ \ \ \ \ \sum_{k^{\sharp}}\frac1{\mathop{k^{\sharp}}^s}=\frac1{k^s}\prod_{\substack{p\mid k\\ p\nmid q/k}}\left(1-\frac1{p^s}\right)^{-1},\notag
\end{align}
and then
\begin{align}
\sum_{(r,q)=k} \frac{1}{r^s}\sum_{(l,q)=1}\frac{c_{l}(r)} {l^{2+\lambda}} =\frac{1}{k^s}\prod_{\substack{p\mid k\\ p\nmid q/k}}\left(1-\frac1{p^s}\right)^{-1} \sum_{(r,q)=1}\frac{1}{r^s}\sum_{(l,q)=1}\frac{c_l(r)}{l^{2+\lambda}}.\notag
\end{align}
After we execute the sum over $l$ by Lemma \ref{lemus}, it becomes
\begin{align}
\frac{1}{k^s}\prod_{\substack{p\mid k\\ p\nmid q/k}}\left(1-\frac1{p^s}\right)^{-1}\frac1{\zeta_q(2+\lambda)}\sum_{(r,q)=1} \frac{\sigma_{-1-\lambda}(r)}{r^s} =\frac{1}{k^s}\prod_{\substack{p\mid k\\ p\nmid q/k}}\left(1-\frac1{p^s}\right)^{-1}\frac{\zeta_q(s)\zeta_q(1+\lambda+s)} {\zeta_q(2+\lambda)},\notag
\end{align}
which establishes the lemma.
\end{proof}

We require an approximate functional equation for the divisor function, which we will apply to separate variables in $\sigma_\lambda(n+r)$.
\begin{lemma}\label{lemafed}
Let $n$ be any positive integer.
For $\lambda\in \mathbb{C}$, we have
\begin{align}\label{eqafed}
\sigma_\lambda(n)=&\sum_{(l,q)=1}\frac{c_{l}(n)}{l^{1-\lambda}} \varpi _\lambda\left(\frac l{\sqrt{n}}\right)+n^\lambda \sum_{(l,q)=1}\frac{c_{l}(n)}{l^{1+\lambda}}\varpi _{-\lambda}\left(\frac l{\sqrt{n}}\right),
\end{align}
where
\begin{align}
\varpi_\lambda(x)=\int_{(a)}x^{-w}\zeta_q(1-\lambda+w)\frac{G(w)}wdw\notag
\end{align}
with $a>|\text{Re}(\lambda)|$ and $G$ defined in Definition \ref{DefG}.
\end{lemma}


\begin{proof}
This lemma is an analog of \cite[Lemma 5.4]{You11} and can be proved similarly. The difference is due to an extra coprime condition $(l,q)=1$ in the sum, and thus Lemma \ref{lemus} should be applied to replace \cite[Lemma 5.2]{You11} in the proof.
\end{proof}


The following lemma is required in the deduction of the main term.
\begin{lemma}\label{lemsumk}
For any nonzero $z\in\mathbb{C}$, we have
\begin{align}\label{eqsumk}
&\prod_{p\mid q}\frac{p-2}{p-1}\sum_{k\mid q}\frac1 {k^z}\prod_{p\mid k}\frac{p-1}{p-2}\prod_{\substack{p\mid k\\ p\nmid q/k}}\left(1-\frac1{p^{z}}\right)^{-1}\sum_{d\mid k}\vp(d)\mu\left(\frac{q}{d}\right)\\
&=\frac{\vp^*(q)q}{\vp(q)q^z}\prod_{p\mid q}\left(1-\frac1{p^z}\right)^{-1}\left(1-\frac1{p^{1-z}}\right).\notag
\end{align}
\end{lemma}
\begin{proof}
It is an exercise to check that both sides of the identity are multiplicative functions of $q$. Thus, it is sufficient to check the identity with $q=p^m$, $m\ge1$. Note that the term $\mu\left(\frac qd\right)$ vanishes unless $\frac qd$ equal to $1$ or a prime number, which implies that $k$ can only take values at $p^m$ and $p^{m-1}$. We check the identity in two cases according to $m\ge2$ and $m=1$. For $m\ge 2$, a direct calculation gives that the left-hand side in \eqref{eqsumk} is equal to
\begin{align}
&\frac{p-2}{p-1}\left[\frac{p-1}{p-2}\frac1 {p^{mz}}\left(1-\frac1{p^{z}}\right)^{-1}\left(\vp\left(p^m\right)-\vp\left(p^{m-1}\right)\right) -\frac{p-1}{p-2}\frac1 {p^{(m-1)z}}\vp\left(p^{m-1}\right)\right]\notag\\
=&\frac{\vp(p^m)}{p^{mz}}\left[\left(1-\frac1{p^z}\right)^{-1}\left(1-\frac1p\right) -\frac1{p^{1-z}}\right]\notag\\ =&\frac{\vp^*(p^m)p^m}{\vp(p^m)p^{mz}}\left(1-\frac1{p^z}\right)^{-1}\left(1-\frac1{p^{1-z}}\right).\notag
\end{align}
For $m=1$, the left-hand side in \eqref{eqsumk} is
\begin{align}
\frac{p-2}{p-1}\left[\frac{p-1}{p-2}\frac1 {p^{z}}\left(1-\frac1{p^{z}}\right)^{-1}\left(\vp\left(p\right)-1\right)-1\right]
=&\frac1{p^{z}}\left(1-\frac1{p^z}\right)^{-1}\frac{(p-2)(p-p^z)}{p-1}\notag\\
=&\frac{\vp^*(p)p}{\vp(p)p^z}\left(1-\frac1{p^z}\right)^{-1}\left(1-\frac1{p^{1-z}}\right).\notag
\end{align}
Thus, we complete the proof.
\end{proof}

\subsection{Separations of variables}\label{secsv}
By writing $n=m+r$, we have
\begin{align}
B_{M,N}&=\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right) \sum_{r\equiv0~(\bmod d)}\\
&\times\sum_{\substack{(m,q)=1\\ (m+r,q)=1}}\frac{\sigma_{\alpha-\beta}(m) \sigma_{\gamma-\delta}(m+r)} {m^{\frac12+\alpha}(m+r)^{\frac12+\gamma}}
V_{\alpha,\beta,\gamma,\delta} \left(\frac{m(m+r)}{q^2}\right) W\left(\frac mM\right)W\left(\frac{m+r}N\right), \notag
\end{align}
which is much the same as $B_{M,N}$ in \cite{You11} but with extra coprime conditions
\[
(m,q)=1,\ \ \ \ \ (m+r,q)=1
\]
in the sum. These coprime conditions do not affect the separations of variables too much. Two necessary variations are the expansion of $\sigma_{\gamma-\delta}(m+r)$ and the expression of the sum over $m$ in terms of $\mathcal{D}_q$-function.

We apply the approximate functional equation \eqref{eqafed} to expand $\sigma_{\gamma-\delta}(m+r)$, and this splits $B_{M,N}$ into two terms
\begin{align}\label{eqBMN}
B_{M,N}=C_{M,N}+\wt{C}_{M,N},
\end{align}
where $C_{M,N}$ is the contribution from the first term of the approximate functional equation and $\wt{C}_{M,N}$ from the second term. More precisely,
\begin{align}\label{eqCMN}
C_{M,N}&=\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right) \sum_{r\equiv0~(\bmod d)}\sum_{(l,q)=1}\frac1{l^{1-\gamma+\delta}}\\
&\times\sum_{\substack{(m,q)=1\\ (m+r,q)=1}}\frac{\sigma_{\alpha-\beta}(m) c_l(m+r)} {m^{\frac12+\alpha}(m+r)^{\frac12+\gamma}}
V_{\alpha,\beta,\gamma,\delta} \left(\frac{m(m+r)}{q^2}\right) W\left(\frac mM\right)W\left(\frac{m+r}N\right)\varpi _{\gamma-\delta}\left(\frac l{\sqrt{m+r}}\right), \notag
\end{align}
and
\begin{align}\label{eqwtCC}
\wt{C}_{M,N}(\alpha,\beta,\gamma,\delta)=C_{M,N}(\alpha,\beta,\delta,\gamma).
\end{align}

For notational convenience, we apply the following two notations
\begin{align}\label{eqdefH1}
H_1(s,u_1,u_2,w)=\frac{G(s)G(w)}{sw} g_{\alpha,\beta,\gamma,\delta}(s) \wt{W}(u_1)\wt{W}(u_2) \zeta_q(1-\gamma+\delta+w),
\end{align}
\begin{align}\label{eqH}
H(s,u_1,u_2,v,w)=\frac{\Gamma(v) \Gamma(\frac12+\gamma+s+u_2-v-\frac w2)}{\Gamma(\frac12+\gamma+s+u_2-\frac w2)}H_1(s,u_1,u_2,w),
\end{align}
both of which decay rapidly as any one of the variables gets large in the imaginary direction.


\begin{lemma}\label{lemC_{M,N}}
With $c_s=c_v=\frac14+\ve$, $c_w=\ve$, and $c_{u_1}=c_{u_2}=0$, we have
\begin{align}\label{eqC_MN}
C_{M,N}&=\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)
\sum_{{r\equiv0~(\bmod d)}}\frac{1}{r^{\frac12+\gamma}}\sum_{(l,q)=1}\frac1{l^{1-\gamma+\delta}}\mathop{\sum\nolimits^*}_{h(\bmod l)}e\left(\frac{hr}{l}\right) \left(\frac1{2\pi i}\right)^5\\
&\times\int_{(c_s)}\int_{(c_{u_1})}\int_{(c_{u_2})}\int_{(c_w)} \int_{(c_v)}\frac{q^{2s}M^{u_1}N^{u_2}} {r^{s+u_2-v-\frac w2}l^w}\mathcal{D}_{q}\left(\tfrac12+\alpha+s+u_1+v,\alpha-\beta,\frac h{l},r\right)\notag\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \times H(s,u_1,u_2,v,w)dvdwdu_2du_1ds. \notag
\end{align}
\end{lemma}
\begin{proof}
The proof is the same as \cite[Lemma 5.5]{You11} since the coprime conditions
\[
(l,q)=1,\ \ \ \ (m,q)=1,\ \ \ \ (m+r,q)=1
\]
are just retained to the end and do not affect the separations of variables.
\end{proof}

By Proposition \ref{lemgED2}, we move $c_v$ to $\ve$, crossing the two poles of $\mathcal{D}_{q}$ at $\tfrac12+s+u_1+v=1-\alpha$ and $\tfrac12+s+u_1+v=1-\beta$, whose residues contribute main terms to $C_{M,N}$. The integral along the new paths contributes an error, denoted by $E_{M,N}$, in which we move $c_s$ to $\ve$ without encountering any poles. Let $P_{M,N}(\alpha,\beta,\gamma,\delta)$ denote the contribution of the former pole, and $P_{M,N}(\beta,\alpha,\gamma,\delta)$  corresponds to the latter pole by symmetry. We have
\begin{align}\label{eqCPPE}
C_{M,N}=P_{M,N}(\alpha,\beta,\gamma,\delta) +P_{M,N}(\beta,\alpha,\gamma,\delta)+E_{M,N},
\end{align}
where the main terms will be calculated in Section \ref{secmainAO} and the error $E_{M,N}$ will be bounded in Section \ref{secError}.

\section{The main term of $A_O$}\label{secmainAO}
To deduce the main term of $A_O$, we should first deduce an exact expression for $P_{M,N}$, which we present in the following.
\begin{lemma}\label{lemP_{M,N}}
For $c_s=\frac14$, $c_{u_1}=c_{u_2}=0$, and $c_w=\ve$, we have
\begin{align}\label{eqP_{M,N}}
P_{M,N}=&\zeta_q(1-\alpha+\beta)\left(\frac1{2\pi i}\right)^4\int_{(c_s)}\int_{(c_{u_1})}\int_{(c_{u_2})} \int_{(c_w)}q^{-\alpha-\gamma-u_1-u_2+\frac w2}\prod_{p\mid q}\left(1-\frac1{p^{1-(\alpha+\gamma+2s+u_1+u_2-\frac w2)}}\right)\\
&\times M^{u_1}N^{u_2} \frac{\zeta(\alpha+\gamma+2s+u_1+u_2-\tfrac w2)\zeta_q(1+\beta+\delta+2s+u_1+u_2+\tfrac w2)} {\zeta_q(2-\alpha+\beta-\gamma+\delta+w)}\notag\\
&\times H_2(s,u_1,u_2,w) dwdu_2du_1ds,\notag
\end{align}
where
\begin{align}\label{eqH_2}
H_2(s,u_1,u_2,w)=\frac{\Gamma(\frac12-\alpha-u_1-s) \Gamma(\alpha+\gamma+2s+u_1+u_2-\frac w2)}{\Gamma(\frac12+\gamma+s+u_2-\frac w2)}H_1(s,u_1,u_2,w).
\end{align}
\end{lemma}

\begin{proof}
By Proposition \ref{lemgED2}, the residue of $\mathcal{D}_q$-function at $\tfrac12+s+u_1+v=1-\alpha$ is
\begin{align}
\frac{\vp(q)}{q}\sum_{b\mid q_r}\frac{\mu(b)}{\vp(b)} l^{-1+\alpha-\beta}\zeta_q(1-\alpha+\beta).\notag
\end{align}
Thus, we have
\begin{align}\label{eqmathcalP_{M,N}(t)}
P_{M,N}=&\zeta_q(1-\alpha+\beta)\frac{\vp(q)}{q\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)\sum_{b\mid q}\frac{\mu(b)}{\vp(b)}\sum_{\substack{r\equiv0 ~(\bmod d)\\ (r,b)=1}}\sum_{(l,q)=1}\frac{c_l(r)}{l^{2-\alpha+\beta-\gamma+\delta}}\\
&\times\left(\frac1{2\pi i}\right)^4\int_{(c_s)}\int_{(c_{u_1})}\int_{(c_{u_2})}\int_{(c_w)} \frac{q^{2s}M^{u_1}N^{u_2}} {r^{\alpha+\gamma+2s+u_1+u_2-\frac w2}l^w}H_2\left(s,u_1,u_2,w\right)dwdu_2du_1ds,\notag
\end{align}
observing that $c_s=\frac14+\ve$, $c_{u_1}=c_{u_2}=0$, $c_w=\ve$, and
\begin{align}
  H_2(s,u_1,u_2,w)=H\left(s,u_1,u_2,\tfrac12-\alpha-u_1-s,w\right).\notag
\end{align}
The ratio of gamma factors in $H_2$ is
\begin{align}
 \frac{\Gamma(\tfrac12-\alpha-u_1-s) \Gamma(\alpha+\gamma+2s+u_1+u_2-\frac w2)}{\Gamma(\frac12+\gamma+s+u_2-\frac w2)},\notag
\end{align}
so we may move $c_{u_2}$ to $1$ without encountering any poles.
Now both summations over $l$ and $r$ converge absolutely, which we will execute by using Lemma \ref{lemsumf}.

For convenience, we apply here two notations
\[
k=(q,r),\ \ \ \ \ \ \ z=\alpha+\gamma+2s+u_1+u_2-\tfrac w2.
\]
With these notations, we have
\begin{align}
P_{M,N}=&\zeta_q(1-\alpha+\beta)\frac{\vp(q)}{q\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)\sum_{b\mid q}\frac{\mu(b)}{\vp(b)}\sum_{\substack{k\mid q_b\\ d\mid k}}\sum_{(r,q)=k}\sum_{(l,q)=1}\frac{c_l(r)}{l^{2-\alpha+\beta-\gamma+\delta}}\notag\\
&\times \left(\frac1{2\pi i}\right)^4\int\cdots\int  \frac{q^{2s}M^{u_1}N^{u_2} } {r^z l^w}H_2\left(s,u_1,u_2,w\right)dwdu_2du_1ds,\notag
\end{align}
where the sum over $l,r$ is ready to apply Lemma \ref{lemsumf}, and it follows
\begin{align}\label{eqP1}
P_{M,N}=&\zeta_q(1-\alpha+\beta)\frac{\vp(q)}{q\vp^*(q)}\sum_{d\mid q} \vp(d)\mu\left(\frac qd\right)\sum_{b\mid q} \frac{\mu(b)}{\vp(b)}\\
&\times \left(\frac1{2\pi i}\right)^4\int\cdots\int q^{2s}\sum_{\substack{k\mid q_b\\ d\mid k}}\frac{1}{k^{z}}\prod_{\substack{p\mid k\\ p\nmid q/k}}\left(1-\frac1{p^z}\right)^{-1}M^{u_1}N^{u_2}\notag\\
&\times\frac{\zeta_q(z)\zeta_q(1+\beta+\delta+2s+u_1+u_2+\tfrac w2)} {\zeta_q(2-\alpha+\beta-\gamma+\delta+w)}H_2(s,u_1,u_2,w) dwdu_2du_1ds.\notag
\end{align}

We first come to the sum over $b$, $d$, and $k$ in the above formula, which is essential to calculate the following arithmetic sum
\begin{align}
\sum_{\substack{d\mid q\\ b\mid q}}\mu(b)\mu\left(\frac qd\right)\frac{\vp(d)}{\vp(b)}\sum_{\substack{k\mid q_b\\ d\mid k}}\frac{1}{k^z}\prod_{\substack{p\mid k\\ p\nmid q/k}}\left(1-\frac1{p^z}\right)^{-1}&=\sum_{k\mid q}\frac{1}{k^z}\prod_{\substack{p\mid k\\ p\nmid q/k}}\left(1-\frac1{p^z}\right)^{-1}\sum_{d\mid k}\vp(d)\mu\left(\frac qd\right)\sum_{b\mid q_k}\frac{\mu(b)}{\vp(b)}.\notag
\end{align}
It is easy to see that
\[
\sum_{b\mid q_k}\frac{\mu(b)}{\vp(b)}=\prod_{p\mid q}\frac{p-2}{p-1}\prod_{p\mid k}\frac{p-1}{p-2}.
\]
After inserting this formula, one notes that the arithmetic sum is ready to apply Lemma \ref{lemsumk}, which gives
\[
\sum_{\substack{d\mid q\\ b\mid q}}\mu(b)\mu\left(\frac qd\right)\frac{\vp(d)}{\vp(b)}\sum_{\substack{k\mid q_b\\ d\mid k}}\frac{1}{k^z}\prod_{\substack{p\mid k\\ p\nmid q/k}}\left(1-\frac1{p^z}\right)^{-1}=\frac{\vp^*(q)q}{\vp(q)q^z}\prod_{p\mid q}\left(1-\frac1{p^z}\right)^{-1}\left(1-\frac1{p^{1-z}}\right).
\]
Applying this into \eqref{eqP1} gives
\begin{align}
P_{M,N}=&\zeta_q(1-\alpha+\beta)\left(\frac1{2\pi i}\right)^4\int\cdots\int q^{-\alpha-\gamma-u_1-u_2+\frac w2}\prod_{p\mid q}\left(1-\frac1{p^{1-z}}\right)M^{u_1}N^{u_2}\notag\\
&\times\frac{\zeta(z)\zeta_q(1+\beta+\delta+2s+u_1+u_2+\tfrac w2)} {\zeta_q(2-\alpha+\beta-\gamma+\delta+w)}H_2(s,u_1,u_2,w) dwdu_2du_1ds,\notag
\end{align}
where we move $c_{u_2}$ back to $0$ to complete our proof, observing that the pole of the Riemann zeta-function at $z=\alpha+\gamma+2s+u_1+u_2-\frac w2=1$ is canceled by the factor $1-\frac1{p^{1-z}}$.
\end{proof}

\begin{remark}
At the end of the above proof, the pole of $\zeta(\alpha+\gamma+2s+u_1+u_2-\frac w2)$ is canceled in $P_{M,N}$, while it should give a non-negligible contribution in \cite{You11}, especially when $M$ and $N$ are close. The pole is the obstruction to calculate $P_{M,N}$ directly in Young \cite{You11}, while it is actually a by-product of removing coprime conditions $(mn,q)=1$ in the original formula. For prime $q$, we can skirt this obstruction by recovering the partition first, but it does not work for general $q$. Actually, even if we recover the partition first, the residue is still very large, especially when $q$ contains small prime factors.
\end{remark}

The upper bound of the main term mentioned in \eqref{eqbmain} is presented especially in the following lemma.
\begin{lemma}\label{lemubP_{M,N}}
For $M\ll Nq^\ve$, we have
\begin{align}
P_{M,N}(\alpha,\beta,\gamma,\delta)\ll M^{\frac12}N^{-\frac12} q^\ve.\notag
\end{align}
\end{lemma}
\begin{proof}
Observing the location of poles of $H_2$, we may move the contours in \eqref{eqP_{M,N}} to
\[
c_s=\ve,\ \ \ \ c_{u_1}=\tfrac12-2\ve,\ \ \ \ c_{u_2}=-\tfrac12+\ve,\ \ \ \ c_w=\ve,
\]
without encountering any poles. Then the lemma follows by a trivial estimation on the integral along these lines.
\end{proof}

We use $\mathcal{A}_O(\alpha,\beta,\gamma,\delta)$ for the main term of $A_O(\alpha,\beta,\gamma,\delta)$ and $\mathcal{A}^*_O(\alpha,\beta,\gamma,\delta)$ for the main term of $A^*_O(\alpha,\beta,\gamma,\delta)$.
We evaluate $\mathcal{A}_O$ in the following proposition and then switch the shifts to give the expression for $\mathcal{A}^*_O$.
\begin{proposition}\label{promathcalB}
We have
\begin{align}
\mathcal{A}_O(\alpha,\beta,\gamma,\delta)=&\mathcal{P}(\alpha,\beta,\gamma,\delta) +\mathcal{P}(\beta,\alpha,\gamma,\delta)+\mathcal{P}(\alpha,\beta,\delta,\gamma) +\mathcal{P}(\beta,\alpha,\delta,\gamma) +O\left(q^{-\frac12+\ve}\right),\notag
\end{align}
where
\begin{align}\label{eqmathcalP}
\mathcal{P}(\alpha,\beta,\gamma,\delta) &=\frac{\zeta_q(1-\alpha+\beta)\zeta_q(1-\gamma+\delta)} {\zeta_q(2-\alpha+\beta-\gamma+\delta)}q^{-\alpha-\gamma}\\
&\times\frac1{2\pi i} \int_{(\frac 14)}\frac{G(s)}{s}g_{\alpha,\beta,\gamma,\delta}(s)\frac{\Gamma(\frac12-\alpha-s) \Gamma(\alpha+\gamma+2s)}{\Gamma(\frac12+\gamma+s)} \notag\\
&\times \prod_{p\mid q}\left(1-\frac1{p^{1-(\alpha+\gamma+2s)}}\right)\zeta(\alpha+\gamma+2s) \zeta_q(1+\beta+\delta+2s) ds.\notag
\end{align}
\end{proposition}
\begin{proof}
We sum up from \eqref{eqAOBMN}, \eqref{eqBMN}, \eqref{eqwtCC}, and \eqref{eqCPPE} that
\begin{align}\label{eqmathcalA}
\mathcal{A}_O(\alpha,\beta,\gamma,\delta)
=P(\alpha,\beta,\gamma,\delta) +P(\beta,\alpha,\gamma,\delta)+P(\alpha,\beta,\delta,\gamma) +P(\beta,\alpha,\delta,\gamma),
\end{align}
where
\begin{align}
P(\alpha,\beta,\gamma,\delta)=\sum_{M\ll Nq^\ve}P_{M,N}(\alpha,\beta,\gamma,\delta).\notag
\end{align}
Observing the location of poles of $H_2$, we can easily take
\[
c_{u_1}=-A, \ \ \ \ \ \ \ \ c_{u_2}=A
\]
for large $A$ in \eqref{eqP_{M,N}} to see that $P_{M,N}\ll q^{-2020}$ with $M\gg Nq^\ve$.
That is to say, we can extend $P$ to the sum over all $M, N$ with a negligible error.
Then we restore the dyadic partition of unity by applying a formula due to Mellin transform
\begin{align}\label{eTqransform}
\sum_{M,N}\left(\frac 1{2\pi i}\right)^2\int_{(c_{u_1})}\int_{(c_{u_2})}F(u_1,u_2)M^{u_1}N^{u_2} \wt{W}(u_1)\wt{W}(u_2)du_1du_2=F(0,0)
\end{align}
for a `nice' function $F$ (see also the proof of \cite[Lemma 6.2]{You11}). To see the nice function, we apply two changes of variables $s\rightarrow s-u_1$ and $w\rightarrow w-2u_1+2u_2$ in \eqref{eqP_{M,N}}, then
\begin{align}
P_{M,N}=&\zeta_q(1-\alpha+\beta)\left(\frac1{2\pi i}\right)^4\int_{(c_s)}\int_{(c_{u_1})}\int_{(c_{u_2})} \int_{(c_w)}q^{-\alpha-\gamma-2u_1+\frac w2}\prod_{p\mid q}\left(1-\frac1{p^{1-(\alpha+\gamma+2s-\frac w2)}}\right)\notag\\
&\times M^{u_1}N^{u_2} \frac{\zeta(\alpha+\gamma+2s-\tfrac w2)\zeta_q(1+\beta+\delta+2s-2u_1+2u_2+\tfrac w2)} {\zeta_q(2-\alpha+\beta-\gamma+\delta-2u_1+2u_2+w)}\notag\\
&\times H_2(s-u_1,u_1,u_2,w-2u_1+2u_2) dwdu_2du_1ds,\notag
\end{align}
where
\begin{align}
H_2(s-u_1,u_1,&u_2,w-2u_1+2u_2)=\frac{\Gamma(\frac12-\alpha-s) \Gamma(\alpha+\gamma+2s-\frac w2)}{\Gamma(\frac12+\gamma+s-\frac w2)}g_{\alpha,\beta,\gamma,\delta}(s-u_1)\notag\\
&\times\wt{W}(u_1)\wt{W}(u_2)\frac{G(s-u_1)G(w-2u_1+2u_2)}{(s-u_1)(w-2u_1+2u_2)}\zeta_q(1-\gamma+\delta-2u_1+2u_2+w).\notag
\end{align}
So after summing over $M$ and $N$ in \eqref{eqP_{M,N}}, we get
\begin{align}\label{eqP}
P=&\zeta_q(1-\alpha+\beta)\left(\frac1{2\pi i}\right)^2\int_{(\frac14)} \int_{(\ve)}q^{-\alpha-\gamma+\frac w2}\prod_{p\mid q}\left(1-\frac1{p^{1-(\alpha+\gamma+2s-\frac w2)}}\right)\\
&\times\frac{\zeta(\alpha+\gamma+2s-\tfrac w2)\zeta_q(1+\beta+\delta+2s+\tfrac w2)} {\zeta_q(2-\alpha+\beta-\gamma+\delta+w)}H_2(s,0,0,w) dwds,\notag
\end{align}
where
\begin{align}
H_2(s,0,0,w)=&\frac{\Gamma(\frac12-\alpha-s) \Gamma(\alpha+\gamma+2s-\frac w2)}{\Gamma(\frac12+\gamma+s-\frac w2)}\frac{G(s)G(w)}{sw} g_{\alpha,\beta,\gamma,\delta}(s) \zeta_q(1-\gamma+\delta+w).\notag
\end{align}

To eliminate the $w$-line integration in \eqref{eqP}, we move it  to $\text{Re}(w)=-1+\ve$. Observing that the pole of $\zeta(1-\gamma+\delta+w)$ is canceled by $G(\gamma-\delta)=0$, we have only crossed the pole at $w=0$, which contributes the main term $\mathcal{P}(\alpha,\beta,\gamma,\delta)$.
A direct estimate shows that the integral along the new lines is bounded by $\ll q^{-\frac12+\ve}$, which establishes the lemma.
\end{proof}

Due to the relationship between $A_O^*$ and $A_O$ mentioned in \eqref{A*OAO}, we switch the parameters to get
\begin{proposition}\label{promathcalB1}
We have
\begin{align}
\mathcal{A}_O^*(\alpha,\beta,\gamma,\delta)=&\mathcal{P}(\gamma,\delta,\alpha,\beta) +\mathcal{P}(\delta,\gamma,\alpha,\beta)+\mathcal{P}(\gamma,\delta,\beta,\alpha) +\mathcal{P}(\delta,\gamma,\beta,\alpha) +O\left(q^{-\frac12+\ve}\right),\notag
\end{align}
where $\mathcal{P}(\alpha,\beta,\gamma,\delta)$ is given by \eqref{eqmathcalP}.
\end{proposition}







\section{Dual off-diagonal terms $A_{\ol{O}}$}\label{secAolO}
We focus on $B_{\ol{M,N}}$ with $M\le N$, and the case $M>N$ follows by swapping the positions of $M, N$.
By writing $r=m+n$, we have
\begin{align}
B_{\ol{M,N}}&=\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right) \sum_{r\equiv0~(\bmod d)}\notag\\
&\times\sum_{\substack{(m,q)=1\\ (r-m,q)=1}}\frac{\sigma_{\alpha-\beta}(m) \sigma_{\gamma-\delta}(r-m)} {m^{\frac12+\alpha}(r-m)^{\frac12+\gamma}}
V_{\alpha,\beta,\gamma,\delta} \left(\frac{m(r-m)}{q^2}\right) W\left(\frac mM\right)W\left(\frac{r-m}N\right),\notag
\end{align}
whose variables we may separate as before. At first, we apply the approximate functional equation \eqref{eqafed} to $\sigma_{\gamma-\delta}(r-m)$ to write
\begin{align}\label{eqolBCC}
B_{\ol{M,N}}=C_{\ol{M,N}}(\alpha,\beta,\gamma,\delta)+C_{\ol{M,N}}(\alpha,\beta,\delta,\gamma),
\end{align}
where
\begin{align}\label{eqColMN}
C_{\ol{M,N}}&=\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right) \sum_{r\equiv0~(\bmod d)}\sum_{(l,q)=1}\frac1{l^{1-\gamma+\delta}}\\
&\times\sum_{\substack{(m,q)=1\\ (r-m,q)=1}}\frac{\sigma_{\alpha-\beta}(m) c_l(r-m)} {m^{\frac12+\alpha}(r-m)^{\frac12+\gamma}}
V_{\alpha,\beta,\gamma,\delta} \left(\frac{m(r-m)}{q^2}\right) W\left(\frac mM\right)W\left(\frac{r-m}N\right)\varpi _{\gamma-\delta}\left(\frac l{\sqrt{r-m}}\right). \notag
\end{align}

We apply the notation
\begin{align}
\ol{H}(s,u_1,u_2,v,w)=\frac{\Gamma(v)\Gamma(\frac12-\gamma-s-u_2+\frac w2)}{\Gamma(\frac12-\gamma-s-u_2+v+\frac w2)}H_1(s,u_1,u_2,w),
\end{align}
which decays rapidly as any one of $\text{Im}(s), ~\text{Im}(u_1), ~\text{Im}(u_2), ~\text{Im}(w)$ gets large. In addition, $\ol{H}(s,u_1,u_2,v,w)$ has at most a polynomial growth on $\text{Im}(v)$, to be specific, it follows from Stirling's approximation that
\[
\frac{\Gamma(v)\Gamma(\frac12-\gamma-s-u_2+\frac w2)}{\Gamma(\frac12-\gamma-s-u_2+v+\frac w2)}\ll (1+|\text{Im}(v)|)^{\text{Re}(\gamma+s+u_2-\frac w2)-\frac12}.
\]
In particular, if $\text{Re}(\gamma+s+u_2-\frac w2)-\frac12<0$, $\ol{H}(s,u_1,u_2,v,w)$ has a polynomial decay on $\text{Im}(v)$.

Then one can separate variables as in Section 7 of \cite{You11} to obtain the following analog of Lemma \ref{lemC_{M,N}}.
\begin{lemma}
With $c_{s}=c_v=\frac14+\ve$, $c_{u_1}=c_{u_2}=0$, and $c_w=2$, we have
\begin{align}\label{eqC_{M,N}}
C_{\ol{M,N}}&=\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)
\sum_{{r\equiv0~(\bmod d)}}\frac{1}{r^{\frac12+\gamma}}\sum_{(l,q)=1}\frac1{l^{1-\gamma+\delta}}\mathop{\sum\nolimits^*}_{h(\bmod l)}e\left(\frac{-hr}{l}\right) \left(\frac1{2\pi i}\right)^5\\
&\times\int_{(c_s)}\int_{(c_{u_1})}\int_{(c_{u_2})}\int_{(c_w)} \int_{(c_v)}\frac{q^{2s}M^{u_1}N^{u_2}} {r^{s+u_2-v-\frac w2}l^w}\mathcal{D}_{q}\left(\tfrac12+\alpha+s+u_1+v,\alpha-\beta,\frac h{l},-r\right)\notag\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \times \ol{H}(s,u_1,u_2,v,w)dvdwdu_2du_1ds. \notag
\end{align}
\end{lemma}

One notes that $C_{\ol{M,N}}$ is almost the same as $C_{M,N}$ except different gamma factors in $\ol{H}$, as well as a different path for the $w$-integral, which is applied to keep the absolute convergence for the $v$-integral. We move the $v$-line of integral to $c_v=\ve$, crossing the two poles of $\mathcal{D}_q$-function at $\frac12+s+u_1+v=1-\alpha$ and  $\frac12+s+u_1+v=1-\beta$ with respective residues $P_{\ol{M,N}}(\alpha,\beta,\gamma,\delta)$ and $P_{\ol{M,N}}(\beta,\alpha,\gamma,\delta)$. For the remained integral, denoted by $E_{\ol{M,N}}$, we move $c_s$ to $\ve$ without encountering any poles. In conclusion, we have
\begin{align}\label{eqolCPPE}
C_{\ol{M,N}}=P_{\ol{M,N}}(\alpha,\beta,\gamma,\delta) +P_{\ol{M,N}}(\beta,\alpha,\gamma,\delta)+E_{\ol{M,N}},
\end{align}
where $P_{\ol{M,N}}$ will be calculated in the next section and $E_{\ol{M,N}}$ will be bounded together with $E_{{M,N}}$ in Section \ref{secError}.

\subsection{The main terms}
We first deduce an exact expression for  $P_{\ol{M,N}}$ in the following lemma.
\begin{lemma}\label{lemolPMN}
For $M\le N$, $c_s=\frac14$, $c_{u_1}=c_{u_2}=0$, and $c_w=\ve$, we have
\begin{align}\label{eqolPMN}
P_{\ol{M,N}}=&\zeta_q(1-\alpha+\beta)\left(\frac1{2\pi i}\right)^4\int_{(c_s)}\int_{(c_{u_1})}\int_{(c_{u_2})} \int_{(c_w)}q^{-\alpha-\gamma-u_1-u_2+\frac w2}\prod_{p\mid q}\left(1-\frac1{p^{1-(\alpha+\gamma+2s+u_1+u_2-\frac w2)}}\right)\\
&\times M^{u_1}N^{u_2}\frac{\zeta(\alpha+\gamma+2s+u_1+u_2-\tfrac w2)\zeta_q(1+\beta+\delta+2s+u_1+u_2+\tfrac w2)} {\zeta_q(2-\alpha+\beta-\gamma+\delta+w)}\notag\\
&\times \ol{H}_2(s,u_1,u_2,w) dwdu_2du_1ds,\notag
\end{align}
where
\begin{align}\label{eqolH_2}
\ol{H}_2(s,u_1,u_2,w)=\frac{\Gamma(\tfrac12-\alpha-u_1-s)\Gamma(\frac12-\gamma-s-u_2+\frac w2)}{\Gamma(1-\alpha-\gamma-2s-u_1-u_2+\frac w2)}H_1(s,u_1,u_2,w).
\end{align}
\end{lemma}


\begin{proof}
After calculating the residue, we have
\begin{align}
P_{\ol{M,N}}=&\zeta_q(1-\alpha+\beta)\frac{\vp(q)}{q\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)\sum_{b\mid q}\frac{\mu(b)}{\vp(b)}\sum_{\substack{r\equiv0 ~(\bmod d)\\ (r,b)=1}}\sum_{(l,q)=1}\frac{c_l(r)}{l^{2-\alpha+\beta-\gamma+\delta}}\notag\\
&\times\left(\frac1{2\pi i}\right)^4\int_{(c_s)}\int_{(c_{u_1})}\int_{(c_{u_2})}\int_{(c_w)} \frac{q^{2s}M^{u_1}N^{u_2}} {r^{\alpha+\gamma+2s+u_1+u_2-\frac w2}l^w}\ol{H}_2\left(s,u_1,u_2,w\right)dwdu_2du_1ds,\notag
\end{align}
which is the same expression as \eqref{eqmathcalP_{M,N}(t)} but with $H_2$ replaced by $\ol{H}_2$ only. The difference between $H_2$ and $\ol{H}_2$ is that the ratio of  gamma factors
\[
 \frac{\Gamma(\tfrac12-\alpha-u_1-s) \Gamma(\alpha+\gamma+2s+u_1+u_2-\frac w2)}{\Gamma(\frac12+\gamma+s+u_2-\frac w2)}
\]
appearing in $H_2$ is replaced by
\[
\frac{\Gamma(\tfrac12-\alpha-u_1-s)\Gamma(\frac12-\gamma-s-u_2+\frac w2)}{\Gamma(1-\alpha-\gamma-2s-u_1-u_2+\frac w2)}
\]
in $\ol{H}_2$.

We initially move the contours of the integration such that
\[
c_s=\ve, \ \ \ \ c_{u_1}=c_{u_2}=0,\ \ \ \ c_w=3\ve,
\]
passing no poles. Then we move the $s$-line of integral to $c_s=\frac12+\ve$, crossing a pole at  $\frac12-\alpha-u_1-s=0$.

We consider the residue first. With $k=(q,r)$, the pole at $\frac12-\alpha-u_1-s=0$ gives
\begin{align}\label{eqdualresidue}
&\zeta_q(1-\alpha+\beta)\frac{\vp(q)}{q\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)\sum_{b\mid q}\frac{\mu(b)}{\vp(b)}\sum_{\substack{k\mid q_b\\ d\mid k}}\sum_{(r,q)=k}\sum_{(l,q)=1}\frac{c_l(r)}{l^{2-\alpha+\beta-\gamma+\delta}}\\
&\times\left(\frac1{2\pi i}\right)^3\int_{(c_{u_1})}\int_{(c_{u_2})}\int_{(c_w)} \frac{q^{1-2\alpha-2u_1}M^{u_1}N^{u_2}} {r^{1-\alpha+\gamma-u_1+u_2-\frac w2}l^w}H_1\left(\tfrac12-\alpha-u_1,u_1,u_2,w\right)dwdu_2du_1,\notag
\end{align}
in which we may move $c_{u_2}$ to $2\ve$  to make both summations over $l$ and $r$ converge absolutely.
So one may execute the sum over $l,r$ with Lemma \ref{lemsumf}, and then execute the sum over $d,b,k$ with Lemma \ref{lemsumk} as before to reduce \eqref{eqdualresidue} to
\begin{align}\label{eqdualresidue1}
\zeta_q(1&-\alpha+\beta)\left(\frac1{2\pi i}\right)^3\int_{(c_{u_1})}\int_{(c_{u_2})}\int_{(c_w)} q^{-\alpha-\gamma-u_1-u_2+\frac w2} \prod_{p\mid q}\left(1-\frac1{p^{\alpha-\gamma+u_1-u_2+\frac w2}}\right)M^{u_1}N^{u_2}\\
&\times \frac{\zeta\left(1-\alpha+\gamma-u_1+u_2-\frac w2\right)\zeta_q\left(2-2\alpha+\beta+\delta-u_1+u_2+\frac w2\right)}{\zeta_q(2-\alpha+\beta-\gamma+\delta+w)}\notag\\
&\times H_1\left(\tfrac12-\alpha-u_1,u_1,u_2,w\right)dwdu_2du_1.\notag
\end{align}

For the integral along the new lines, the sum over $l,r$ converges absolutely, and then similar calculations as for \eqref{eqP_{M,N}} show that it has the same expression as the right-hand side of \eqref{eqolPMN} except $c_s=\frac12+\ve$. We move $c_s$ back to $\frac14$ to obtain the expression for $P_{\ol{M,N}}$ and cross the pole at $\frac12-\alpha-u_1-s=0$ again, observing that the pole of the zeta-function is canceled. An easy calculation shows that this residue is equal to \eqref{eqdualresidue1}, so the two residues cancel, and we establish the lemma.
\end{proof}

The same to Lemma \ref{lemolPMN}, it is easy to see the following upper bound for $P_{\ol{M,N}}$ from \eqref{eqolPMN}.
\begin{lemma}
For $M\le N$,
we have
\begin{align}\label{eqboundPolMN}
P_{\ol{M,N}}\ll M^{-\frac12}N^{\frac12}q^\ve.
\end{align}
\end{lemma}

For the case $M>N$, we switch $M, N$ and the shifts in \eqref{eqolPMN} to get
\begin{align}\label{eqPolNM}
P_{\ol{M,N}}(\alpha,\beta,\gamma,\delta)=P_{\ol{N,M}}(\gamma,\delta,\alpha,\beta)
\end{align}
with $P_{\ol{N,M}}(\gamma,\delta,\alpha,\beta)$ given by \eqref{eqolPMN}.

Now we assemble $P_{\ol{M,N}}$ to evaluate $A_{\ol{O}}$. Let $\mathcal{A}_{\ol{O}}(\alpha,\beta,\gamma,\delta)$ denote the main term of $A_{\ol{O}}(\alpha,\beta,\gamma,\delta)$, and we deduce it in the following proposition.
\begin{proposition}\label{proolmathcalA}
We have
\begin{align}
\mathcal{A}_{\ol{O}}(\alpha,\beta,\gamma,\delta)=&\ol{\mathcal{P}}(\alpha,\beta,\gamma,\delta) +\ol{\mathcal{P}}(\beta,\alpha,\gamma,\delta)+\ol{\mathcal{P}}(\alpha,\beta,\delta,\gamma) +\ol{\mathcal{P}}(\beta,\alpha,\delta,\gamma) +O\left(q^{-\frac13+\ve}\right),\notag
\end{align}
where
\begin{align}\label{eqmathcalP}
\ol{\mathcal{P}}(\alpha,\beta,\gamma,\delta) &=\frac{\zeta_q(1-\alpha+\beta)\zeta_q(1-\gamma+\delta)} {\zeta_q(2-\alpha+\beta-\gamma+\delta)}q^{-\alpha-\gamma}\\
&\times\frac1{2\pi i} \int_{(\frac 14)}\frac{G(s)}{s}g_{\alpha,\beta,\gamma,\delta}(s)\frac{\Gamma(\frac12-\alpha-s) \Gamma(\frac12-\gamma-s)}{\Gamma(1-\alpha-\gamma-2s)} \notag\\
&\times \prod_{p\mid q}\left(1-\frac1{p^{1-(\alpha+\gamma+2s)}}\right)\zeta(\alpha+\gamma+2s) \zeta_q(1+\beta+\delta+2s) ds.\notag
\end{align}
\end{proposition}
\begin{proof}
It follows from \eqref{eqAolOBMN}, \eqref{eqolBCC}, and \eqref{eqolCPPE} that
\begin{align}
\mathcal{A}_{\ol{O}}(\alpha,\beta,\gamma,\delta)=\ol{P}(\alpha,\beta,\gamma,\delta)+\ol{P}(\beta,\alpha,\gamma,\delta)+\ol{P}(\alpha,\beta,\delta,\gamma)+\ol{P}(\beta,\alpha,\delta,\gamma),\notag
\end{align}
where
\[
\ol{P}(\alpha,\beta,\gamma,\delta)=\sum_{M,N}P_{\ol{M,N}}(\alpha,\beta,\gamma,\delta).
\]

Note from \eqref{eqolPMN} and \eqref{eqPolNM} that the expressions for the two cases $M\le N$ and $M>N$ are not identical, so we should not recover the partition with these expressions directly. The difference is due to the $w$-line of integral. We eliminate the $w$-line of integral in \eqref{eqolPMN} by moving it to $c_w=-\frac14+\ve$. Since the pole of $\zeta(1-\gamma+\delta+w)$ is cancelled by $G(\gamma-\delta)=0$, we only cross the pole at $w=0$. For the remained integral, we move $c_s$ to $\frac 16$ followed by $c_w\rightarrow -\frac23+\ve$, obtaining an upper bound $O\left(q^{-\frac13+\ve}\right)$ by a trivial estimate. Thus, we have
\begin{align}\label{eqPolMN}
P_{\ol{M,N}}&=\zeta_q(1-\alpha+\beta)\zeta(1-\gamma+\delta)\left(\frac1{2\pi i}\right)^3\int_{(c_s)}\int_{(c_{u_1})}\int_{(c_{u_2})} q^{-\alpha-\gamma-u_1-u_2} M^{u_1}N^{u_2}\\
&\times\prod_{p\mid q}\left(1-\frac1{p^{1-(\alpha+\gamma+2s+u_1+u_2)}}\right)\frac{\zeta(\alpha+\gamma+2s+u_1+u_2)\zeta_q(1+\beta+\delta+2s+u_1+u_2)} {\zeta_q(2-\alpha+\beta-\gamma+\delta)}\notag\\
&\times \frac{\Gamma(\tfrac12-\alpha-u_1-s)\Gamma(\frac12-\gamma-s-u_2)}{\Gamma(1-\alpha-\gamma-2s-u_1-u_2)}\frac{G(s)}{s} g(s) \wt{W}(u_1)\wt{W}(u_2)du_2du_1ds+O\left(q^{-\frac13+\ve}\right).\notag
\end{align}
With the expression \eqref{eqPolMN}, it is easy to check that
\[
P_{\ol{M,N}}(\alpha,\beta,\gamma,\delta)=P_{\ol{N,M}}(\gamma,\delta,\alpha,\beta)+O\left(q^{-\frac13+\ve}\right),
\]
so we observe from \eqref{eqPolNM} that the expression \eqref{eqPolMN} holds for the $M<N$ case too.
After the variable change $s\rightarrow s-\frac{u_1+u_2}2$ in \eqref{eqPolMN}, the zeta functions do not rely on $u_1$ and $u_2$ any more, and the product of gamma factors and $G$ becomes
\[
\frac{\Gamma\left(\tfrac12-\alpha-s-\frac{u_1-u_2}2\right)\Gamma(\frac12-\gamma-s+\frac{u_1-u_2}2)}{\Gamma(1-\alpha-\gamma-2s)}G\left(s-\frac{u_1+u_2}2\right),
\]
which decays rapidly in the imaginary direction of $u_1$ and $u_2$.
Thus, we sum over $M, N$ and recover the dyadic partition with \eqref{eTqransform} to get
\[
\ol{P}(\alpha,\beta,\gamma,\delta)=\ol{\mathcal{P}}(\alpha,\beta,\gamma,\delta)+O\left(q^{-\frac13+\ve}\right),
\]
with $\ol{\mathcal{P}}$ given by \eqref{eqmathcalP}, which establishes the lemma.
\end{proof}






\section{The main term of the Theorem}\label{secmain}

In this section, we assemble all main terms arising from both off-diagonal terms and diagonal terms to deduce the main term of Theorem \ref{thmmains}, which proceeds much the same as \cite{You11}, and we sketch it here.

At first, we combine $A_O$, $A^*_O$, and $A_{\ol{O}}$ to deduce the main contribution of off-diagonal terms in $A_1$, which follows easily by \eqref{eqA_1} with Propositions \ref{promathcalB}, \ref{promathcalB1}, and \ref{proolmathcalA} that
\begin{align}\label{eqmoA_1}
Q(\alpha,\beta,\gamma,\delta)+Q(\beta,\alpha,\gamma,\delta)+Q(\alpha,\beta,\delta,\gamma)+Q(\beta,\alpha,\delta,\gamma)+O(q^{-\frac13+\ve}),
\end{align}
where
\[
Q(\alpha,\beta,\gamma,\delta)=\mathcal{P}(\alpha,\beta,\gamma,\delta)+\mathcal{P}(\gamma,\delta,\alpha,\beta)+\ol{\mathcal{P}}(\alpha,\beta,\gamma,\delta).
\]

\begin{lemma}
We have
\begin{align}
Q(\alpha,\beta,\gamma,\delta)=&\frac{\zeta_q(1-\alpha+\beta)\zeta_q(1-\gamma+\delta)}{\zeta_q(2-\alpha+\beta-\gamma+\delta)}\left(\frac q\pi\right)^{-\alpha-\gamma}\frac1{2\pi i}\int_{(\frac14)}\frac{G(s)}{s}\pi^{2s}g_{\alpha,\beta,\gamma,\delta}(s)\notag\\
&\times\zeta_q(1-\alpha-\gamma-2s)\zeta_q(1+\beta+\delta+2s)\frac{\Gamma\left(\frac{\frac12-\alpha-s}{2}\right)}{\Gamma\left(\frac{\frac12+\alpha+s}{2}\right)} \frac{\Gamma\left(\frac{\frac12-\gamma-s}{2}\right)}{\Gamma\left(\frac{\frac12+\gamma+s}{2}\right)}ds.\notag
\end{align}
\end{lemma}

\begin{proof}
The lemma can be proved identically to \cite[Lemma 8.1]{You11}, and the only difference is to apply the functional equation
\[
\prod_{p\mid q}\left(1-\frac1{p^{1-(\alpha+\gamma+2s)}}\right)\zeta(\alpha+\gamma+2s)=\frac1{\pi^{\frac12-\alpha-\gamma-2s}} \frac{\Gamma\left(\frac{1-\alpha-\gamma-2s}{2}\right)}{\Gamma\left(\frac{\alpha+\gamma+2s}{2}\right)}\zeta_q(1-\alpha-\gamma-2s)
\]
in place of the functional equation of $\zeta(\alpha+\gamma+2s)$ in the original proof.
\end{proof}

Secondary, by the relationship \eqref{eqrA_1A_{-1}}, one can easily deduce the main contribution of off-diagonal terms in $A_{-1}$ from \eqref{eqmoA_1},  by switching the signs of the shifts and multiplying by $X_{\alpha,\beta,\gamma,\delta}$, which is
\[
Q_{-}(\alpha,\beta,\gamma,\delta)+Q_{-}(\beta,\alpha,\gamma,\delta)+Q_{-}(\alpha,\beta,\delta,\gamma)+Q_{-}(\beta,\alpha,\delta,\gamma)+O(q^{-\frac13+\ve})
\]
with
\[
Q_{-}(\alpha,\beta,\gamma,\delta)=X_{\alpha,\beta,\gamma,\delta}Q(-\gamma,-\delta,-\alpha,-\beta).
\]

Now we combine main terms of off-diagonal terms in $A_1$ and $A_{-1}$. We couple $Q(\alpha,\beta,\gamma,\delta)$ with $Q_{-}(\beta,\alpha,\delta,\gamma)$ and the other three pairs by switching the shifts, which follow easily from the following lemma.
\begin{lemma}
We have
\begin{align}
Q(\alpha,\beta,\gamma,\delta)+Q_{-}(\beta,\alpha,\delta,\gamma)=X_{\alpha,\gamma}\frac{\zeta_q(1-\alpha+\beta)\zeta_q(1-\alpha-\gamma)\zeta_q(1+\beta+\delta)\zeta_q(1-\gamma+\delta)} {\zeta_q(2-\alpha+\beta-\gamma+\delta)}.
\end{align}
\end{lemma}

\begin{proof}
The lemma can be proved identically to \cite[Lemma 8.3]{You11} with $\zeta$ replaced by $\zeta_q$ only.
\end{proof}

By combining these  main terms of off-diagonal terms with the diagonal terms in Lemma \ref{lemD}, we get the main term of $M$ that
\begin{align}
&Z_{q}(\alpha,\beta,\gamma,\delta) +X_{\alpha,\beta,\gamma,\delta}Z_{q}(-\gamma,-\delta,-\alpha,-\beta) +X_{\alpha,\gamma}Z_{q}(\beta,-\gamma,\delta,-\alpha) \notag\\
&+X_{\beta,\gamma}Z_{q}(\alpha,-\gamma,\delta,-\beta)\notag+X_{\alpha,\delta}Z_{q}(\beta,-\delta,\gamma,-\alpha)  +X_{\beta,\delta}Z_{q}(\alpha,-\delta,\gamma,-\beta)+O(q^{-\frac13+\ve}),\notag
\end{align}
which gives the main term of Theorem \ref{thmmains}.


\section{Error terms with $M$ and $N$ close}\label{secError}
This section is devoted to bounding the error terms $E_{M,N}$ and $E_{\ol{M,N}}$. For ease presentation, we set all the shifts $\alpha, \beta$, etc., equal to $0$, and one can easily generalize the arguments below to handle small enough nonzero parameters.

\subsection{The error term $E_{M,N}$}\label{secEMN}
We first recall that $E_{M,N}$ is given by the right-hand side of \eqref{eqC_MN}, but with the contours of the integration at
\begin{align}\label{eqcontour}
c_s=c_w=c_v=\ve,\ \ \ \ c_{u_1}=c_{u_2}=0.
\end{align}
After applying the functional equation to $\mathcal{D}_{q}$, we replace the sum over $a$ and $b$ with a factor $q^\ve$ and rename $a_1,b_1$ as $a, b$ for notational convenience, subsequently, $\chi'_0$ denotes the principal character modulo $ab$. Then we have
\begin{align}
E_{M,N}\ll q^\ve\mathop{\sum\sum\sum}_{\substack{a,b,d\mid q\\ (b,d)=1}}\mu^2(ab)(|E_+|+|E_-|),\notag
\end{align}
where\begin{align}
E_\pm&=\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}\sum_{(l,q)=1}\frac1{lab}\mathop{\sum\nolimits^*}_{i(\bmod{a})} \mathop{\sum\nolimits^*}_{j(\bmod{b})} \mathop{\sum\nolimits^*}_{h(\bmod{l})} e\left(\frac{jr}{b}\right)e\left(\frac{hr}{l}\right)\notag\\
&\times\left(\frac1{2\pi i}\right)^5
\int_{(c_s)}\int_{(c_w)}\int_{(c_{u_1})}\int_{(c_{u_2})} \int_{(c_v)}\frac{q^{2s}M^{u_1}N^{u_2}} {r^{\frac12+s+u_2-v-\frac w2}l^{w}\left(lab\right)^{2s+2u_1+2v}}\notag\\
&\times D\left(\tfrac12-s-u_1-v,0,\pm\frac{\ol{h_{i,j}}}{lab},\chi'_0\right) H_\pm(s,u_1,u_2,v,w)dvdu_2du_1dwds, \notag
\end{align}
and where
\begin{align}\label{eqHpm}
H_\pm(s,u_1,u_2,v,w)=2(2\pi)^{-1+2s+2u_1+2v}\Gamma(\tfrac12-s-u_1-v)^2 H(s,u_1,u_2,v,w)S_\pm
\end{align}
with
\begin{align}
S_+=1,\ \ \ \ \ S_-=\sin(\pi(s+u_1+v)).\notag
\end{align}
We will treat the sums over $r, j, h$ carefully by writing them in terms of the Kloosterman sum and the Ramanujan sum. Here $H_\pm$ have rapid decay on all variables since $H$ does and that the exponential decay of the gamma factors cancels the exponential growth of $S_-$.

\begin{notation}
Note that we should treat $E_{\pm}$ with such $d$ and square-free $a,b$ that $(a,b)=1$ and $(b,d)=1$, especially, there are $ab\le q$ and $bd\le q$. These conditions should be kept in mind in the following.
\end{notation}

We now move $c_{u_1}$ to $-\frac12-3\ve$ and expand $D$ into absolutely convergent Dirichlet series, then
\begin{align}\label{eqE_pm}
E_\pm=&\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}\sum_{(l,q)=1}\frac1{lab}\sum_{(m,ab)=1}\frac{d(m)\mathcal{S}_\pm(m,r)}{r^{\frac12}m^{\frac12}} \\
&\times\left(\frac1{2\pi i}\right)^5
\int_{(c_s)}\int_{(c_w)}\int_{(c_{u_1})}\int_{(c_{u_2})} \int_{(c_v)}\frac{q^{2s}M^{u_1}N^{u_2}m^{s+u_1+v}} {r^{s+u_2-v-\frac w2}l^w\left(lab\right)^{2s+2u_1+2v}}\notag\\
&\times H_\pm(s,u_1,u_2,v,w)dvdu_2du_1dwds \notag
\end{align}
with the exponential sum
\begin{align}\label{eqS_pm}
\mathcal{S}_\pm(m, r)=\mathop{\sum\nolimits^*}_{i(\bmod a)} \mathop{\sum\nolimits^*}_{j(\bmod b)} \mathop{\sum\nolimits^*}_{h(\bmod l)} e\left(\frac{hr}{l}\right) e\left(\frac{jr}{b}\right)e\left(\pm\frac{m\ol{h_{i,j}}}{lab}\right).
\end{align}
With a negligible error, we may restrict the $m$-sum in $E_{\pm}$ to $m\le M^{-1}Na^2b^2q^\ve$, which we assume in the rest of this section.
To see this, we move the paths of the integration in $E_{\pm}$ to
\[c_s=2,\ \ \ \ c_{u_1}=-A,\ \ \ \ c_{u_2}=A,\ \ \ \ c_w=2A,\ \ \ \ c_v=\ve
\]
for large $A$, crossing no poles. Then trivial estimates show that the contribution of all terms with $m\gg M^{-1}Na^2b^2q^\ve$ is negligible.

We simplify the exponential sum $\mathcal{S}_\pm(m, r)$ in the following lemma.

\begin{lemma}\label{lemes}
Let $a, b$, and $l$ be three integers coprime with each other. We have
\begin{align}
\mathcal{S}_\pm(m, r)=c_{a}(m)S(r\ol{a}^2,\pm m;lb).\notag
\end{align}
\end{lemma}
\begin{proof}
We recall from Proposition \ref{lemgED2} that
\begin{align}
\left\{
\begin{aligned}
&h_{i,j}\equiv ilb  \pmod{a},\\
&h_{i,j}\equiv jla \pmod{b},\\
&h_{i,j}\equiv hab \pmod{l}.\notag
\end{aligned}
\right.
\end{align}
Applying Lemma \ref{lemmaspe(x)}, we split the last exponential function as
\begin{align}
e\left(\pm\frac{m\ol{h_{i,j}}}{lab}\right) = e\left(\pm\frac{m\ol{i}(\ol{lb})^2}{a}\right) e\left(\pm\frac{m\ol{(hb+jl)}\ol{a}^2}{lb}\right).\notag
\end{align}
Inserting this into \eqref{eqS_pm}, we sum over $i$ first, which turns out to be the Ramanujan sum $c_{a}(m)$. Besides, we apply the Chinese Remainder Theorem to write the sums over $j$ and $h$ in terms of Kloosterman sums to get
\[
\mathop{\sum\nolimits^*}_{j(\bmod {b})}\mathop{\sum\nolimits^*}_{h(\bmod{l})} e\left(\frac{jr}{b}\right)e\left(\frac{hr}{l}\right)e\left(\pm\frac{m\ol{(hb+jl)}\ol{a}^2}{lb}\right)=S(r\ol{a}^2,\pm m;lb),
\]
and the lemma follows.
\end{proof}


Note that we may change the order of the summations in \eqref{eqE_pm} as we want. More precisely, we may move the lines of integration in \eqref{eqE_pm} without encountering any poles so that all integrations and summations converge absolutely, taking
\[
c_s=c_v=\ve, \ \ \ c_{u_1}=-\tfrac34,\ \ \ c_{u_2}=c_w=\tfrac32
\]
as an example, and similarly, we can move the contours back to \eqref{eqcontour} after changing the order.
By applying Lemma \ref{lemes} and changing the order of the summations, we have
\begin{align}
E_\pm=&\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}\sum_{\substack{m\le M^{-1}Na^2b^2q^\ve\\ (m,ab)=1}} \frac{d(m)c_{a}(m)}{r^{\frac12}m^{\frac12}} \notag\\
&\times\sum_{(l,q)=1}\frac{S(r\ol{a}^2,\pm m;lb)}{lab}\left(\frac1{2\pi i}\right)^5
\int_{(c_s)}\int_{(c_w)}\int_{(c_{u_1})}\int_{(c_{u_2})} \int_{(c_v)}\frac{q^{2s}M^{u_1}N^{u_2}m^{s+u_1+v}} {r^{s+u_2-v-\frac w2}l^w\left(lab\right)^{2s+2u_1+2v}}\notag\\
&\times H_\pm(s,u_1,u_2,v,w)dvdu_2du_1dwds. \notag
\end{align}
Note that
\[
c_{a}(m)=\mu(a)
\]
for $(m,a)=1$.\footnote{The coprime condition $(m,a)=1$ benefits from our kicking out nonessential terms in the deduction of the functional equation for $\mathcal{D}_q$. If not the case, an extra factor from $c_a(m)$ will lead to an obstruction.}
We split the coprime condition $(l,q)=1$ into $(l,a)=1$ and $(l,q_{a})=1$ and remove $(l,q_{a})=1$ by M\"{o}bius inversion
\[
\sum_{c\mid (l,q_{a})}\mu(c).
\]
Then it follows
\begin{align}
E_\pm\ll \sum_{c\mid q_a}|\mathcal{E}_{\pm}|,\notag
\end{align}
where
\begin{align}
\mathcal{E}_{\pm}=&\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{r\equiv0~(\bmod d)\\ (r,b)=1}} \sum_{\substack{m\le M^{-1}Na^2b^2q^\ve\\ (m,ab)=1}} \frac{d(m)}{r^{\frac12}m^{\frac12}} \sum_{(l,a)=1}\frac{S(r\ol{a}^2,\pm m;lbc)}{labc}\notag\\
&\times\left(\frac1{2\pi i}\right)^5
\int_{(c_s)}\int_{(c_w)}\int_{(c_{u_1})}\int_{(c_{u_2})} \int_{(c_v)}\frac{q^{2s}M^{u_1}N^{u_2}m^{s+u_1+v}} {r^{s+u_2-v-\frac w2}(lc)^w\left(labc\right)^{2s+2u_1+2v}}\notag\\
&\times H_\pm(s,u_1,u_2,v,w)dvdu_2du_1dwds.\notag
\end{align}
We should keep in mind that $a,b,c$ are square-free divisors of $q$, satisfying $(a,bc)=1$.

Let $\mathcal{E}_{\pm}(R,M^*)$ be the same expression as $\mathcal{E}_{\pm}$ but with $r$ and $m$ restricted to respective dyadic segments $R\le r\le 2R$ and $M^*\le m\le 2M^*$.
Now we are left to estimate $\mathcal{E}_{\pm}(R,M^*)$ with $M^*\ll M^{-1}Na^2b^2q^\ve$.

We apply \eqref{eqSS} with $Q=a^2bc$, $\tau=a^2$, and $s=bc$, which gives
\[
S(r\ol{a}^2,\pm m;lbc)=e\left(\mp m\frac{\ol{bc}}{a^2}\right)S_{\infty, 1/bc}(r, \pm m; \gammaup).
\]
Note that $S_{\infty, 1/bc}(r,\pm m; \gammaup)$ is defined if and only if, $\gammaup=labc$ with integer $l$ coprime with $a$, that is to say, the sum over $l$ in $\mathcal{E}_{\pm}$ is equal to the sum over $\gammaup$.

Let
\begin{align}\label{eqphi}
\phi_{\pm}(x)=\left(\frac1{2\pi i}\right)^5
\int\cdots \int\frac{a^wb^wq^{2s}M^{u_1}N^{u_2}} {r^{2s+u_1+u_2}m^{\frac w2}} \left(\frac{x}{4\pi}\right)^{2s+2u_1+2v+w}H_\pm(s,u_1,u_2,v,w)dvdu_2du_1dwds
\end{align}
with the contours of the integration given as in \eqref{eqcontour}. By moving the integration to suitable paths without encountering any poles, one can easily check that $\phi_{\pm}(0)=0$ as well as $\phi_{\pm}^{(j)}(x)\ll (1+x)^{-2-\ve}$, $j=0,1,2$.
Then replacing the sum over $l$ by the sum over $\gammaup$ gives
\begin{align}\label{eqERM}
\mathcal{E}_{\pm}(R,M^*)=&\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{R\le r\le 2R\\ {\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}}}\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}} \frac{e\left(\mp m\frac{\ol{bc}}{a^2}\right)d(m)}{r^{\frac12}m^{\frac12}} \Sigma_{\pm}(r,m),
\end{align}
where
\[
\Sigma_{\pm}(r,m)=\sum_{\gammaup}^{\Gamma}\frac{S_{\infty, 1/bc}(r, \pm m; \gammaup)}{\gammaup}\phi_{\pm}\left(\frac{4\pi\sqrt{mr}}{\gammaup}\right).
\]
We apply the Kuznetsov formula with $\mathfrak{a}=\infty$ and $\mathfrak{b}=1/bc$ to get
\begin{align}
\Sigma_{+}(r,m)=&\sum_{\substack{k\ge2\\ k~ \text{even}}}\sum_{f\in\mathcal{B}_k(Q)}\tilde{\phi}_+(k)\Gamma(k)\sqrt{rm}\ol{\rho}_{f}(r)\rho_f(\mathfrak{b},m)\notag\\
&+\sum_{f\in\mathcal{B}(Q)}\hat{\phi}_+(\kappa_f)\frac{\sqrt{rm}}{\cosh(\pi \kappa_f)}\ol{\rho}_f(r)\rho_f(\mathfrak{b},m)\notag\\ &+\frac1{4\pi}\sum_{\mathfrak{c}}\int_{-\infty}^{\infty}\hat{\phi}_+(\kappa)\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\rho_{\mathfrak{b},\mathfrak{c}}(m,\kappa)d\kappa\notag
\end{align}
and
\begin{align}
\Sigma_{-}(r,m)=&\sum_{f\in\mathcal{B}(Q)}\breve{\phi}_-(\kappa_f)\frac{\sqrt{rm}}{\cosh(\pi \kappa_f)}\ol{\rho}_f(r)\rho_f(\mathfrak{b},-m)\notag\\ &+\frac1{4\pi}\sum_{\mathfrak{c}}\int_{-\infty}^{\infty}\breve{\phi}_-(\kappa)\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\rho_{\mathfrak{b},\mathfrak{c}}(-m,\kappa)d\kappa.\notag
\end{align}
Applying these into \eqref{eqERM}, we rewrite
\[
\mathcal{E}_{\pm}(R,M^*)=\mathcal{E}_{h\pm}(R,M^*)+\mathcal{E}_{m\pm}(R,M^*)+\mathcal{E}_{c\pm}(R,M^*)
\]
to correspond to the holomorphic forms, the Maass forms, and the Eisenstein series. Obviously, there is $\mathcal{E}_{h-}(R,M^*)=0$.


By integral transforms of Bassel functions as in \cite[Section 9.3]{You11}, we have
\begin{align}
\tilde{\phi}_+(k)=\left(\frac1{2\pi i}\right)^5\int\cdots\int\frac{a^wb^wq^{2s}M^{u_1}N^{u_2}}{r^{2s+u_1+u_2} m^{\frac w2}}\tilde{H}_+(s,u_1,u_2,v,w,k)dvdu_1du_2dwds,\notag
\end{align}
where
\begin{align}
\tilde{H}_+(s,u_1,u_2,v,w,k)=&\frac{\Gamma(\frac{k-1}2+s+u_1+v+\tfrac w2)}{\Gamma(\frac{k+1}2-s-u_1-v-\tfrac w2)}\Gamma(\tfrac12-s-u_1-v)^2\notag\\
&\times \frac{\Gamma(v)\Gamma(\frac12+s+u_2-v-\frac w2)}{\Gamma(\frac12+s+u_2-\frac w2)} \frac{G(s)G(w)}{sw}g(s)\wt{W}(u_1)\wt{W}(u_2)\zeta(1+w)c^*\notag
\end{align}
with $c^*$ being some bounded factors like powers of $2$, $\pi$, etc., which do not affect the convergence of the integrals. Also, there exist
\begin{align}
\hat{\phi}_+(\kappa)=\left(\frac1{2\pi i}\right)^5\int\cdots\int\frac{a^wb^wq^{2s}M^{u_1}N^{u_2}}{r^{2s+u_1+u_2} m^{\frac w2}}\hat{H}_+(s,u_1,u_2,v,w,\kappa)dvdu_1du_2dwdsp\notag
\end{align}
with
\begin{align}
\hat{H}_+(s,u_1,u_2,v,w,\kappa)=&\cos(\pi(s+u_1+v+\tfrac w2))\Gamma(s+u_1+v+\tfrac w2+i\kappa)\Gamma(s+u_1+v+\tfrac w2-i\kappa)\notag\\
&\times \Gamma(\tfrac12-s-u_1-v)^2\frac{\Gamma(v)\Gamma(\frac12+s+u_2-v-\frac w2)}{\Gamma(\frac12+s+u_2-\frac w2)}\notag\\
&\times\frac{G(s)G(w)}{sw}g(s)\wt{W}(u_1)\wt{W}(u_2)\zeta(1+w)c^*,\notag
\end{align}
and
\begin{align}\label{eqphi-}
\breve{\phi}_-(\kappa)=\left(\frac1{2\pi i}\right)^5\int\cdots\int\frac{a^wb^wq^{2s}M^{u_1}N^{u_2}}{r^{2s+u_1+u_2} m^{\frac w2}}\cosh(\pi\kappa)\breve{H}_-(s,u_1,u_2,v,w,\kappa)dvdu_1du_2dwds
\end{align}
with
\begin{align}
\breve{H}_-(s,u_1,u_2,v,w,\kappa)=&\Gamma(s+u_1+v+\tfrac w2+i\kappa)\Gamma(s+u_1+v+\tfrac w2-i\kappa)\notag\\
&\times \sin(\pi(s+u_1+v))\Gamma(\tfrac12-s-u_1-v)^2\frac{\Gamma(v)\Gamma(\frac12+s+u_2-v-\frac w2)}{\Gamma(\frac12+s+u_2-\frac w2)}\notag\\
&\times\frac{G(s)G(w)}{sw}g(s)\wt{W}(u_1)\wt{W}(u_2)\zeta(1+w)c^*.\notag
\end{align}

All treatments of $\mathcal{E}_{h+}(R,M^*)$, $\mathcal{E}_{m\pm}(R,M^*)$, and $\mathcal{E}_{c\pm}(R,M^*)$ follow in the same way by spectral large sieve inequalities given in Lemma \ref{lemsls} with slight variations, which are because of their different Fourier coefficients and Gamma factors. We show in detail for $\mathcal{E}_{c\pm}(R,M^*)$ and point out necessary variations for other cases.

\begin{proposition}\label{proE}
We have
\begin{align}
\mathcal{E}_{c\pm}(R,M^*)\ll q^{-\frac12+\ve}\left(\frac NM\right)^{\frac12}R^{-\ve}.
\end{align}
\end{proposition}

We only show in detail for $\mathcal{E}_{c-}$ since the case of $\mathcal{E}_{c+}$ is similar.
\begin{proof}
We recall that
\begin{align}
\mathcal{E}_{c-}(R,M^*)=&\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{R\le r\le 2R\\ {\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}}}\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}} \frac{e\left(m\frac{\ol{bc}}{a^2}\right)d(m)}{r^{\frac12}m^{\frac12}} \notag\\
&\times\frac1{4\pi}\sum_{\mathfrak{c}}\int_{-\infty}^{\infty}\breve{\phi}_-(\kappa)\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\rho_{\mathfrak{b},\mathfrak{c}}(-m,\kappa)d\kappa.\notag
\end{align}
Let $\mathcal{E}_{c-}(R,M^*;K)$ be the same expression as $\mathcal{E}_{c-}(R,M^*)$ but with the integral over $\kappa$ restricted to $[K,2K]$,
which after our applying \eqref{eqphi-} evolves into
\begin{align}\label{eqformulaEc-}
\mathcal{E}_{c-}&(R,M^*;K)\ll\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{R\le r\le 2R\\ {\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}}}\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}} \frac{e\left(m\frac{\ol{bc}}{a^2}\right)d(m)}{r^{\frac12}m^{\frac12}} \sum_{\mathfrak{c}}\int_{K}^{2K}\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\rho_{\mathfrak{b},\mathfrak{c}}(-m,\kappa)\\
&\ \ \ \times \left(\frac1{2\pi i}\right)^5\int\cdots\int\frac{a^w b^w q^{2s}M^{u_1}N^{u_2}}{r^{2s+u_1+u_2} m^{\frac w2}}\cosh(\pi\kappa)\breve{H}_-(s,u_1,u_2,v,w,\kappa)dvdu_1du_2dwds d\kappa.\notag
\end{align}
The rapid decay of $\breve{H}_-$ allows for truncation of all integrals except the one over $\kappa$ at height $(qK)^\ve$ with a negligible error.
After this truncation, we may move the integration in $\mathcal{E}_{c-}(R,M^*;K)$ without encountering any poles to any contours with
\begin{align}\label{eqcontourE1}
 c_v=c_w=\ve,\ \ \ \ \
 \left\{
\begin{aligned}
&c_s\ge\ve,\\
&-\ve\le c_s+c_{u_1}\le \tfrac 12-2\ve,\\
&c_s+c_{u_2}\ge -\tfrac12+2\ve
\end{aligned}
\ \ \ \ \ \text{for any}~ K
\right.
\end{align}
and
\begin{align}\label{eqcontourE2}
c_v=c_w=\ve,\ \ \ \ \
\left\{
\begin{aligned}
&c_s\ge\ve,\\
&c_s+c_{u_1}\le \tfrac 12-2\ve,\\
&c_s+c_{u_2}\ge -\tfrac12+2\ve
\end{aligned}
\ \ \ \ \ \text{especially for large}~ K\gg q^\ve.
\right.
\end{align}

We apply Stirling's approximation to see that
\begin{align}\label{eqboundGamma}
\cosh(\pi\kappa)\Gamma(s+u_1+v+\tfrac w2+i\kappa)\Gamma(s+u_1+v+\tfrac w2-i\kappa)\ll q^\ve K^{-1+2c_s+2c_{u_1}+2c_v+c_w}.
\end{align}
By Cauchy-Schwarz inequality,
\begin{align}\label{eqspectralsum}
&\sum_{\substack{R\le r\le 2R\\ {\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}}}\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}} \frac{e\left(m\frac{\ol{bc}}{a^2}\right)d(m)}{r^{\frac12+2s+u_1+u_2}m^{\frac12+\frac w2}} \sum_{\mathfrak{c}}\int_{K}^{2K}\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\rho_{\mathfrak{b},\mathfrak{c}}(-m,\kappa)d\kappa\\
&\ll\left(\sum_{\mathfrak{c}}\int_{K}^{2K}\frac{1}{\cosh(\pi \kappa)}\bigg|\sum_{\substack{R\le r\le 2R\\ {\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}}}\frac1{r^{\frac12+2s+u_1+u_2}}\sqrt{r}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\bigg|^2d\kappa\right)^{\frac12}\notag\\
&\ \ \ \ \times\left(\sum_{\mathfrak{c}}\int_{K}^{2K}\frac{1}{\cosh(\pi \kappa)}\bigg|\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}} \frac{e\left(m\frac{\ol{bc}}{a^2}\right)d(m)}{m^{\frac12+\frac w2}}\sqrt{m}\rho_{\mathfrak{b},\mathfrak{c}}(-m,\kappa)\bigg|^2\right)^{\frac12}.\notag
\end{align}
We bound the second factor first, which follows directly by the large sieve inequality of Lemma \ref{lemsls} that
\begin{align}
\ll \left(K^2+\frac{M^*}{Q}\right)^{\frac12}\left(\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}}\frac1{m^{1+\ve}}\right)^{\frac12}\ll \left(K+\frac {N^{\frac12}}{M^{\frac12}}\right)b^{\frac12}q^\ve\notag
\end{align}
since $\muup(\mathfrak{b})=Q^{-1}$, $M^*\ll M^{-1}Na^2b^2q^\ve$, and $Q=a^2bc$. For the first factor, we apply the change of variables $r\rightarrow d^{\sharp}r$, where $d^{\sharp}$ is largest divisor of $r$ owning the same distinct prime factors as $d$. Then we have
\[
\sum_{\substack{R\le r\le 2R\\ {\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}}}=\sum_{d^{\sharp}}\sum_{\substack{R\le r\le 2R\\ {\substack{r=r'd^{\sharp}\\ (r',bd^{\sharp})=1}}}}.
\]
Note that the number of $d^{\sharp}$ is not more than $R^\ve$ since $d^{\sharp}$ must have the same distinct prime factors as $d$.
Then by \eqref{eqErho}, the first factor is at most
\begin{align}\label{eqsumr}
&\sum_{d^{\sharp}}\left(\sum_{\mathfrak{c}}\int_{K}^{2K}\frac{1}{\cosh(\pi \kappa)}\bigg|\sum_{\substack{R\le r\le 2R\\ {\substack{r=r'd^\sharp\\ (r',bd^\sharp)=1}}}}\frac1{r^{\frac12+2s+u_1+u_2}}\sqrt{r}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\bigg|^2d\kappa\right)^{\frac12}\\
&\ll q^\ve\sum_{d^{\sharp}}\sum_{d_1\mid (d^\sharp,Q)}\left(\sum_{\mathfrak{c}}\int_{K}^{2K}\frac{1}{\cosh(\pi \kappa)}\bigg|\sum_{\substack{R/d^\sharp\le r'\le 2R/d^\sharp\\ (r',bd^\sharp)=1}}\frac1{(r'd^\sharp)^{\frac12+2s+u_1+u_2}}\sqrt{d_1r'}\ol{\rho}_{\mathfrak{c}}(d_1r',\kappa)\bigg|^2d\kappa\right)^{\frac12}.\notag
\end{align}
 After applying the large sieve inequality of Lemma \ref{lemsls}, we find that it is bounded by
\begin{align}\label{eqrddsharp}
&\ll q^\ve\sum_{d^{\sharp}}\sum_{d_1\mid (d^\sharp,Q)}\left(K^2+\frac{d_1R}{d^\sharp Q}\right)^{\frac12}\left(\sum_{\substack{R/d^\sharp\le r'\le 2R/d^\sharp\\ (r',bd^\sharp)=1}}\frac1{(r'd^\sharp)^{1+4c_s+2c_{u_1}+2c_{u_2}}}\right)^{\frac12}\\
&\ll q^\ve\sum_{d^{\sharp}}\sum_{d_1\mid (d^\sharp,Q)}\left(K^2+\frac{d_1R}{d^\sharp Q}\right)^{\frac12}\left(\frac{1}{d^\sharp R^{4c_s+2c_{u_1}+2c_{u_2}}}\right)^{\frac12}\notag\\
&\ll \frac{q^\ve}{R^{2c_s+c_{u_1}+c_{u_2}-\ve}}\left(\frac{K}{d^\frac12}+\frac{R^{\frac12}}{dQ^{\frac12}}\sum_{d_1\mid (d,Q)}d_1^{\frac12}\right).\notag
\end{align}
In \eqref{eqrddsharp}, we have replaced $d^\sharp$ by $d$ in the final step at no cost, observing that $d^\sharp$ has the same prime factors as $d$ and satisfies $d\le d^\sharp$.
In conclusion, we bound the spectral sum in \eqref{eqspectralsum} as
\begin{align}\label{eqboundss}
&\sum_{\substack{R\le r\le 2R\\ {\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}}}\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}} \frac{e\left(m\frac{\ol{bc}}{a^2}\right)d(m)}{r^{\frac12+2s+u_1+u_2}m^{\frac12+\frac w2}} \sum_{\mathfrak{c}}\int_{K}^{2K}\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\rho_{\mathfrak{b},\mathfrak{c}}(-m,\kappa)d\kappa\\
&\ll \frac{q^\ve}{R^{2c_s+c_{u_1}+c_{u_2}-\ve}}\left(K+\frac {N^{\frac12}}{M^{\frac12}}\right)\left(\frac{Kb^{\frac12}}{d^\frac12}+\frac{R^{\frac12}}{d Q^{\frac12}}\sum_{d_1\mid (d,Q)}d_1^{\frac12}b^{\frac12}\right)\notag\\
&\ll \frac{q^{\frac12+\ve}}{dR^{2c_s+c_{u_1}+c_{u_2}-\ve}}\left(K+\frac {N^{\frac12}}{M^{\frac12}}\right) \left(K+q^{-\frac12}R^{\frac12}\right),\notag
\end{align}
observing that $d_1 b\le Q$ and $db\le q$, two results deduced directly from
\[
d_1\mid Q,\ \ \ \ \ b\mid Q,\ \ \ \ \ (d_1,b)=1
\]
and
\[
d\mid q,\ \ \ \ \ b\mid q,\ \ \ \ \ (d,b)=1
\]
respectively.

Applying \eqref{eqboundGamma} and \eqref{eqboundss} into \eqref{eqformulaEc-} gives
\begin{align}\label{eqER}
\mathcal{E}_{c-}(R,M^*;K)
\ll q^{-\frac12+2c_s+\ve}M^{c_{u_1}}N^{c_{u_2}}K^{-1+2c_s+2c_{u_1}+3\ve}\left(K+\frac {N^{\frac12}}{M^{\frac12}}\right)\left(K+q^{-\frac12}R^{\frac12}\right)R^{-2c_s-c_{u_1}-c_{u_2}+\ve}.
\end{align}
If $R\le K^2q$, it follows that
\begin{align}\label{eqEsmallR}
\mathcal{E}_{c-}(R,M^*;K)\ll q^{-\frac12+2c_s+\ve}M^{c_{u_1}}N^{c_{u_2}}K^{2c_s+2c_{u_1}+3\ve}\left(K+\frac {N^{\frac12}}{M^{\frac12}}\right) R^{-2c_s-c_{u_1}-c_{u_2}+\ve}.
\end{align}
When $K\ll M^{-\frac12}N^{\frac12}q^\ve$, we take the contours  as in \eqref{eqcontourE1} with
\[
c_s=\ve,\ \ \ c_{u_1}=\tfrac 12-3\ve,\ \ \ c_{u_2}=-\tfrac12+3\ve;
\]
When $K\gg M^{-\frac12}N^{\frac12}q^\ve$, we take the contours  as in \eqref{eqcontourE2} with
\[
c_s=\ve,\ \ \ c_{u_1}=-\tfrac 12-3\ve,\ \ \ c_{u_2}=\tfrac12+3\ve.
\]
With these in \eqref{eqEsmallR}, a direct calculation gives
\begin{align}
\mathcal{E}_{c-}(R,M^*;K)\ll q^{-\frac12+\ve}\left(\frac NM\right)^{\frac12}R^{-\ve}K^{-\ve}.\notag
\end{align}
If $R\ge K^2q$, then \eqref{eqER} evolves into
\begin{align}\label{eqElargeR}
\mathcal{E}_{c-}(R,M^*;K)\ll q^{-1+2c_s+\ve}M^{c_{u_1}}N^{c_{u_2}}K^{-1+2c_s+2c_{u_1}+3\ve} \left(K+\frac {N^{\frac12}}{M^{\frac12}}\right)R^{\frac12-2c_s-c_{u_1}-c_{u_2}+\ve}.
\end{align}
When $K\ll M^{-\frac12}N^{\frac12}q^\ve$, we take the contours  as in \eqref{eqcontourE1} with
\[
c_s=\tfrac 14+\ve,\ \ \ c_{u_1}=\tfrac 14-3\ve,\ \ \  c_{u_2}=-\tfrac 14+3\ve;
\]
When $K\gg M^{-\frac12}N^{\frac12}q^\ve$, we take the contours  as in \eqref{eqcontourE2} with
\[
c_s=\tfrac14+\ve,\ \ \ c_{u_1}=-\tfrac 14-3\ve,\ \ \ c_{u_2}=\tfrac14+3\ve.
\]
With these in \eqref{eqElargeR}, a direct calculation gives
\begin{align}
\mathcal{E}_{c-}(R,M^*;K)\ll q^{-\frac12+\ve}\left(\frac NM\right)^{\frac14}R^{-\ve}K^{-\ve}.\notag
\end{align}
Thus, we establish the lemma by summing over $K$.
\end{proof}


We now come to $\mathcal{E}_{h+}(R,M^*;K)$. Note that
\begin{align}\label{eqHG}
\frac{\Gamma(\frac{k-1}2+s+u_1+v+\tfrac w2)}{\Gamma(\frac{k+1}2-s-u_1-v-\tfrac w2)}\ll (K+|v|)^{-1+2c_s+2c_{u_1}+2c_v+c_w},
\end{align}
so the same treatment as the Eisenstein spectrum case gives the same bound for the holomorphic spectrum.

We bound the Maass spectrum in the following.
\begin{proposition}\label{proM}
We have
\begin{align}
\mathcal{E}_{m\pm}(R,M^*)\ll q^{-\frac12+\theta+\ve}\left(\frac NM\right)^{\frac12}R^{-\ve}.
\end{align}
\end{proposition}
\begin{proof}
The proposition follows by the same argument used for the Eisenstein spectrum. Following the proof of Proposition \ref{proE}, we only need a replacement of \eqref{eqErho} by \eqref{eqMrho} to separate $d^\sharp$ from the Fourier coefficients. We can bound the extra factor $\lambda_{f^*}$ arising from the replacement with the known bound due to  Kim and Sarnak \cite{Kim03} that
\[
\lambda_{f^*}(n)\ll n^{\theta+\ve}
\]
with $\theta=7/64$. More precisely, the following estimate due to \eqref{eqMrho} is applied in the analog of \eqref{eqsumr} to replace \eqref{eqErho}
\[
\sum_{\substack{R\le r\le 2R\\ r=r'd^\sharp\\ (r',bd^\sharp)=1}}\frac1{r^{\frac12+2s+u_1+u_2}}\sqrt{r}\ol{\rho}_f(r)\ll {d^{\sharp}}^{\theta+\ve}\sum_{d_1\mid (d^{\sharp},Q)}\Bigg|\sum_{\substack{R/d^\sharp\le r'\le 2R/d^\sharp\\ (r',bd^\sharp)=1}}\frac1{(r'd^\sharp)^{\frac12+2s+u_1+u_2}}\sqrt{d_1r'}\rho_f(d_1r')\Bigg|.
\]
Note that the extra factor ${d^\sharp}^{\theta+\ve}$ contributes at most a factor $q^{\theta+\ve}$ to $\mathcal{E}_{m\pm}(R,M^*)$ after our replacing $d^\sharp$ with $d$ as in \eqref{eqrddsharp}. Then the rest of the proof follows identically to the Eisenstein spectrum.
\end{proof}


At last, we sum up from the holomorphic spectrum, the Maass spectrum, and the Eisenstein spectrum to obtain that
\[
\mathcal{E}_{\pm}(R,M^*)\ll q^{-\frac12+\theta+\ve}\left(\frac NM\right)^{\frac12}R^{-\ve},
\]
which establishes Theorem \ref{thmEMN} by a summation over $R$ and $M^*$.


\subsection{The dual error term $E_{\ol{M,N}}$}
We apply the functional equation to $\mathcal{D}_q$, and the same treatment as in Section \ref{secEMN} gives
\begin{align}
E_{\ol{M,N}}\ll q^\ve\mathop{\sum\sum\sum}_{\substack{a,b,d\mid q \\ (b,d)=1}}\mu^2(ab)(|\ol{E}_+|+|\ol{E}_-|),\notag
\end{align}
where\begin{align}\label{eqdefolE}
\ol{E}_\pm=&\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}\sum_{(l,q)=1}\frac1{lab}\sum_{(m,ab)=1}\frac{d(m)\mathcal{S}_\pm(m,-r)}{r^{\frac12}m^{\frac12}}\\
&\times\left(\frac1{2\pi i}\right)^5
\int_{(c_s)}\int_{(c_w)}\int_{(c_{u_1})}\int_{(c_{u_2})} \int_{(c_v)}\frac{q^{2s}M^{u_1}N^{u_2}m^{s+u_1+v}} {r^{s+u_2-v-\frac w2}l^w\left(lab\right)^{2s+2u_1+2v}}\notag\\
&\times \ol{H}_\pm(s,u_1,u_2,v,w)dvdu_2du_1dwds \notag
\end{align}
with $\mathcal{S}_\pm(m,-r)$ defined by \eqref{eqS_pm} and  $\ol{H}_{\pm}$ given by an expression identical to \eqref{eqHpm} except that $H$ is replaced by $\ol{H}$.

Recalling that the ratio of the Gamma factors in $\ol{H}$ is
\[
\frac{\Gamma(v)\Gamma(\frac12-s-u_2+\frac w2)}{\Gamma(\frac12-s-u_2+v+\frac w2)},
\]
we can move the path of the integration so that both summations over $l$ and $m$ converge absolutely, and thus we can change their order. However, we can not make the $r$-sum convergent absolutely without crossing any poles, \footnote{Note that the sum over $r$ in \eqref{eqdefolE} is in the order we require, which we have prepared before separating variables in $C_{\ol{M,N}}$.} and we can not truncate the sum over $m$ to $m\le M^{-1}Na^2b^2q^\ve$ as we have done for $E_{M,N}$.

By the same computations as in Section \ref{secEMN}, we may reduce the problem to bounding
\begin{align}\label{eqolERM}
\ol{\mathcal{E}}_{\pm}(R,M^*)=&\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{R\le r\le 2R\\ r\equiv0~(\bmod d)\\ (r,b)=1}}\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}} \frac{e\left(\pm m\frac{\ol{bc}}{a^2}\right)d(m)}{r^{\frac12}m^{\frac12}} \ol{\Sigma}_{\pm}(r,m),
\end{align}
where
\[
\ol{\Sigma}_{\pm}(r,m)=\sum_{\gammaup}\frac{S_{\infty, 1/bc}(r, \mp m; \gammaup)}{\gammaup}\Theta _{\pm}\left(\frac{4\pi\sqrt{mr}}{\gammaup}\right)
\]
with $\Theta_\pm$ given by an expression identical to \eqref{eqphi} except that $H_{\pm}$ is replaced by $\ol{H}_\pm$. Note that we have lost the bound $M^*\ll M^{-1}Na^2b^2q^\ve$ here.

We apply the Kuznetsov formula as before, then
\begin{align}
\ol{\Sigma}_{-}(r,m)=&\sum_{\substack{k\ge2\\ k~ \text{even}}}\sum_{f\in\mathcal{B}_k(Q)}\tilde{\Theta}_-(k)\Gamma(k)\sqrt{rm}\ol{\rho}_{f}(r)\rho_f(\mathfrak{b},m)\notag\\
&+\sum_{f\in\mathcal{B}(Q)}\hat{\Theta}_-(\kappa_f)\frac{\sqrt{rm}}{\cosh(\pi \kappa_f)}\ol{\rho}_f(r)\rho_f(\mathfrak{b},m) \notag\\ &+\frac1{4\pi}\sum_{\mathfrak{c}}\int_{-\infty}^{\infty}\hat{\Theta}_-(\kappa)\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\rho_{\mathfrak{b},\mathfrak{c}}(m,\kappa)d\kappa\notag
\end{align}
and
\begin{align}
\ol{\Sigma}_{+}(r,m)=&\sum_{f\in\mathcal{B}(Q)}\breve{\Theta}_+(\kappa_f)\frac{\sqrt{rm}}{\cosh(\pi \kappa_f)}\ol{\rho}_f(r)\rho_f(\mathfrak{b},-m) \notag\\
&+\frac1{4\pi}\sum_{\mathfrak{c}}\int_{-\infty}^{\infty}\breve{\Theta}_+(\kappa)\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\rho_{\mathfrak{b},\mathfrak{c}}(-m,\kappa)d\kappa.\notag
\end{align}
Applying these into \eqref{eqolERM}, we rewrite
\[
\ol{\mathcal{E}}_{\pm}(R,M^*)=\ol{\mathcal{E}}_{h\pm}(R,M^*)+\ol{\mathcal{E}}_{m\pm}(R,M^*)+\ol{\mathcal{E}}_{c\pm}(R,M^*)
\]
to correspond to the holomorphic forms, the Maass forms and the Eisenstein series. Also, there is $\ol{\mathcal{E}}_{h+}(R,M^*)=0$.


Integral transforms as before give
\begin{align}
\tilde{\Theta}_-(k)=\left(\frac1{2\pi i}\right)^5\int\cdots\int\frac{a^wb^wq^{2s}M^{u_1}N^{u_2}}{r^{2s+u_1+u_2} m^{\frac w2}}\tilde{\ol{H}}_-(s,u_1,u_2,v,w,k)dvdu_1du_2dwds\notag
\end{align}
with
\begin{align}
\tilde{\ol{H}}_-(s,u_1,u_2,v,w,k)=&\frac{\Gamma(\frac{k-1}2+s+u_1+v+\tfrac w2)}{\Gamma(\frac{k+1}2-s-u_1-v-\tfrac w2)}\sin(\pi(s+u_1+v))\Gamma(\tfrac12-s-u_1-v)^2\notag\\
&\times \frac{\Gamma(v) \Gamma(\frac12-s-u_2+\frac w2)}{\Gamma(\frac12-s-u_2+v+\frac w2)} \frac{G(s)G(w)}{sw}g(s)\wt{W}(u_1)\wt{W}(u_2)\zeta(1+w)c^*.\notag
\end{align}
In addition,
\begin{align}
\hat{\Theta}_-(\kappa)=\left(\frac1{2\pi i}\right)^5\int\cdots\int\frac{a^wb^wq^{2s}M^{u_1}N^{u_2}}{r^{2s+u_1+u_2} m^{\frac w2}}\hat{\ol{H}}_-(s,u_1,u_2,v,w,\kappa)dvdu_1du_2dwds\notag
\end{align}
with
\begin{align}
\hat{\ol{H}}_-(s,u_1,u_2,v,w,\kappa)=&\cos(\pi(s+u_1+v+\tfrac w2))\Gamma(s+u_1+v+\tfrac w2+i\kappa)\Gamma(s+u_1+v+\tfrac w2-i\kappa)\notag\\
&\times \sin(\pi(s+u_1+v))\Gamma(\tfrac12-s-u_1-v)^2\frac{\Gamma(v) \Gamma(\frac12-s-u_2+\frac w2)}{\Gamma(\frac12-s-u_2+v+\frac w2)}\notag\\
&\times\frac{G(s)G(w)}{sw}g(s)\wt{W}(u_1)\wt{W}(u_2)\zeta(1+w)c^*,\notag
\end{align}
and
\begin{align}
\breve{\Theta}_+(\kappa)=\left(\frac1{2\pi i}\right)^5\int\cdots\int\frac{a^wb^wq^{2s}M^{u_1}N^{u_2}}{r^{2s+u_1+u_2} m^{\frac w2}}\cosh(\pi\kappa)\breve{\ol{H}}_+(s,u_1,u_2,v,w,\kappa)dvdu_1du_2dwds\notag
\end{align}
with
\begin{align}
\breve{\ol{H}}_+(s,u_1,u_2,v,w,\kappa)=&\Gamma(s+u_1+v+\tfrac w2+i\kappa)\Gamma(s+u_1+v+\tfrac w2-i\kappa)\notag\\
&\times \Gamma(\tfrac12-s-u_1-v)^2\frac{\Gamma(v)\Gamma(\frac12-s-u_2+\frac w2)}{\Gamma(\frac12-s-u_2+v+\frac w2)}\notag\\
&\times\frac{G(s)G(w)}{sw}g(s)\wt{W}(u_1)\wt{W}(u_2)\zeta(1+w)c^*.\notag
\end{align}

\begin{proposition}\label{prodcE}
We have
\begin{align}\label{eqdcE}
\ol{\mathcal{E}}_{c\pm}(R,M^*)\ll q^{-\frac12+\ve}\left(\frac NM\right)^{\frac12}R^{-\ve}{M^*}^{-\ve}.
\end{align}
\end{proposition}

\begin{proof}
We come to $\ol{\mathcal{E}}_{c+}(R,M^*)$ first, which is matched with $\breve{\Theta}_+(\kappa)$. We write
\[
\ol{\mathcal{E}}_{c+}(R,M^*)=\sum_K\ol{\mathcal{E}}_{c+}(R,M^*;K),
\]
where
\begin{align}
\ol{\mathcal{E}}_{c+}(&R,M^*;K)\ll\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{R\le r\le 2R\\ {\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}}}\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}} \frac{e\left(m\frac{\ol{bc}}{a^2}\right)d(m)}{r^{\frac12}m^{\frac12}} \sum_{\mathfrak{c}}\int_{K}^{2K}\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\notag\\
&\times \rho_{\mathfrak{b},\mathfrak{c}}(-m,\kappa)\left(\frac1{2\pi i}\right)^5\int\cdots\int\frac{a^wb^wq^{2s}M^{u_1}N^{u_2}}{r^{2s+u_1+u_2} m^{\frac w2}}\breve{\ol{H}}_+(s,u_1,u_2,v,w,\kappa)dvdu_1du_2dwds d\kappa.\notag
\end{align}
Note in $\breve{\ol{H}}_+$ that the gamma factor $\Gamma(\tfrac12-s-u_1-v)^2$ assures the repaid decay on $v$, so we may truncate all the integrals except $\kappa$ at height $(qK)^\ve$ with a negligible error. The treatment is much the same as the $\mathcal{E}_{c-}(R,M^*;K)$ case, but there are some differences. Especially, due to the new ratio of gamma factors
\[
\frac{\Gamma(v)\Gamma(\frac12-s-u_2+\frac w2)}{\Gamma(\frac12-s-u_2+v+\frac w2)},
\]
it does not seem that we can bound $\ol{\mathcal{E}}_{c+}(R,M^*;K)$ well just moving the paths with no poles crossed.
We move $c_v$ to $-2\ve$, crossing the pole at $v=0$ with the residue denoted by $res_{c+}(R,M^*;K)$, which has the same expression as $\ol{\mathcal{E}}_{c+}(R,M^*;K)$ but with the $v$-integral disappeared and $\breve{\ol{H}}_+$ replaced by
\begin{align}
\Gamma(s+u_1+\tfrac w2+i\kappa)\Gamma(s+u_1+\tfrac w2-i\kappa)\Gamma(\tfrac12-s-u_1)^2\frac{G(s)G(w)}{sw}g(s)\wt{W}(u_1)\wt{W}(u_2)\zeta(1+w)c^*.\notag
\end{align}
Note that we may move the $u_2$-line of integral freely in $res_{c+}(R,M^*;K)$.
Now
\[
\ol{\mathcal{E}}_{c+}(R,M^*;K)=\ol{\mathcal{E}}'_{c+}(R,M^*;K)+res_{c+}(R,M^*;K),
\]
where $\ol{\mathcal{E}}'_{c+}(R,M^*;K)$ is the same integral as $\ol{\mathcal{E}}_{c+}(R,M^*;K)$ but with $c_v=-2\ve$.


With $c_w=2\ve$ for $M^*\le K^2Q$ and $c_w=1+\ve$ for $M^*> K^2Q$ in the spectral sum
\[
\sum_{\substack{R\le r\le 2R\\ {\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}}}\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}} \frac{e\left(m\frac{\ol{bc}}{a^2}\right)d(m)}{r^{\frac12+2s+u_1+u_2}m^{\frac12+\frac w2}} \sum_{\mathfrak{c}}\int_{K}^{2K}\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\rho_{\mathfrak{b},\mathfrak{c}}(m,\kappa)d\kappa,
\]
we have the following upper bound for it
\begin{align}\label{eqboundss2}
\ll \left\{
\begin{aligned}
&\frac{q^{\frac12+\ve}K}{dR^{2c_s+c_{u_1}+c_{u_2}-\ve}}\left(K+q^{-\frac12}R^{\frac12}\right){M^*}^{-\ve},\ \ \ \ \ \ \ \ \ \text{for}\ \ M^*\le K^2Q,\\
&\frac{q^{\frac12+\ve}}{abdR^{2c_s+c_{u_1}+c_{u_2}-\ve}}\left(K+q^{-\frac12}R^{\frac12}\right){M^*}^{-\ve},\ \ \ \ \ \text{for} \ \ M^*> K^2Q.
\end{aligned}
\right.
\end{align}
This upper bound is deduced much the same as \eqref{eqboundss}, but the following bound
\begin{align}
\ll\left\{
\begin{aligned}
&K{M^*}^{-\ve}, \ \ \ \ \ \ \ \ \ \ \ \ \text{for} \ \ M^*\le K^2Q,\\
&a^{-1}b^{-\frac12}{M^*}^{-\ve},\ \ \ \ \text{for} \ \ M^*> K^2Q\notag
\end{aligned}
\right.
\end{align}
is applied to bound the second factor in \eqref{eqspectralsum} instead, which is a direct result of the large sieve inequality of Lemma \ref{lemsls}.

To treat $\ol{\mathcal{E}}'_{c+}(R,M^*;K)$, we take
\begin{align}
 c_v=-2\ve,\ \ \ \ \ \ \ \ c_w=\left\{
\begin{aligned}
&2\ve, \ \ \ \ \ \ \ \text{for} \ \ M^*\le K^2Q,\\
&1+\ve,\ \ \ \text{for} \ \ M^*> K^2Q.\notag
\end{aligned}
\right.
\end{align}
Applying \eqref{eqboundGamma} and \eqref{eqboundss2} into the expression of $\ol{\mathcal{E}}'_{c+}(R,M^*;K)$ gives \footnote{With $c_v>0$, there is an analog of \eqref{eq+olE} for $\ol{\mathcal{E}}'_{c+}(R,M^*;K)$, where we can not make both exponents of $K$ and $R$ negative simultaneously just moving the paths crossing no poles.}
\begin{align}\label{eq+olE}
\ol{\mathcal{E}}'_{c+}(R,M^*;K)\ll q^{-\frac12+2c_s+\ve}M^{c_{u_1}}N^{c_{u_2}}K^{2c_s+2c_{u_1}-2\ve}\left(K+q^{-\frac12}R^{\frac12}\right)R^{-2c_s-c_{u_1}-c_{u_2}+\ve}{M^*}^{-\ve},
\end{align}
by which we will deduce the upper bound for $\ol{\mathcal{E}}'_{c+}(R,M^*;K)$ with appropriate $c_s, c_{u_1}, c_{u_2}$.

For small $K\ll q^\ve$, we need not care for the exponent of $K$, but we should avoid poles of $\Gamma(s+u_1+v+\tfrac w2\pm i\kappa)$ in moving the paths.
For such $K$, we take
\begin{align}
 c_{u_1}=c_{u_2}=\ve,\ \ \ \ \ \ \ \ c_s=\left\{
\begin{aligned}
&\ve, \ \ \ \ \ \ \ \ \ \ \text{for} \ \ R\le K^2q,\\
&\tfrac14+\ve,\ \ \ \text{for} \ \ R> K^2q.\notag
\end{aligned}
\right.
\end{align}
to get
\[
\ol{\mathcal{E}}'_{c+}(R,M^*;K)\ll q^{-\frac12+\ve}R^{-\ve}K^{-\ve}{M^*}^{-\ve}.
\]


When $K\gg q^\ve$ is large, we will not cross the poles of $\Gamma(s+u_1+v+\tfrac w2\pm i\kappa)$, but should keep the exponent of $K$ negative.
If $K\gg q^\ve$ and $R\le K^2q$, we move the contours so that
\[
c_s=\tfrac34\ve, \ \ c_{u_1}=-\tfrac12, \ \ c_{u_2}=\tfrac12,
\]
with which it follows
\[
\ol{\mathcal{E}}'_{c+}(R,M^*;K)\ll q^{-\frac12+\ve}\left(\frac NM\right)^{\frac12}R^{-\ve}K^{-\ve}{M^*}^{-\ve}.
\]
If $K\gg q^\ve$ and $R> K^2q$, we move the contours so that
\[
c_s=\tfrac14+\tfrac34\ve, \ \ c_{u_1}=-\tfrac14, \ \ c_{u_2}=\tfrac14,
\]
and then a direct calculate gives
\[
\ol{\mathcal{E}}'_{c+}(R,M^*;K)\ll q^{-\frac12+\ve}\left(\frac NM\right)^{\frac14}R^{-\ve}K^{-\ve}{M^*}^{-\ve}.
\]


The treatment for $res(R,M^*;K)$ is similar and easier since $c_{u_2}$ is free. The analog to \eqref{eq+olE} for the residue is
\[
res_{c+}(R,M^*;K)\ll q^{-\frac12+2c_s+\ve}M^{c_{u_1}}N^{c_{u_2}}K^{2c_s+2c_{u_1}+2\ve}\left(K+q^{-\frac12}R^{\frac12}\right)R^{-2c_s-c_{u_1}-c_{u_2}+\ve}{M^*}^{-\ve}.
\]
Then the same discussions, with $c_{u_1}$ reducing several $\ve$ and $c_{u_2}$ adding several $\ve$, show the same bound for the residue. Thus, summing these bounds over $K$ gives the bound \eqref{eqdcE} for $\ol{\mathcal{E}}_{c+}(R,M^*)$.

Similarly, we write
\[
\ol{\mathcal{E}}_{c-}(R,M^*)=\sum_K\ol{\mathcal{E}}_{c-}(R,M^*;K)
\]
with
\begin{align}
\ol{\mathcal{E}}_{c-}(&R,M^*;K)\ll\frac{\vp(d)}{\vp^*(q)}
\sum_{\substack{R\le r\le 2R\\ {\substack{r\equiv0~(\bmod d)\\ (r,b)=1}}}}\sum_{\substack{M^*\le m\le2M^*\\ (m,ab)=1}} \frac{e\left(- m\frac{\ol{bc}}{a^2}\right)d(m)}{r^{\frac12}m^{\frac12}} \sum_{\mathfrak{c}}\int_{K}^{2K}\frac{\sqrt{rm}}{\cosh(\pi \kappa)}\ol{\rho}_{\mathfrak{c}}(r,\kappa)\notag\\
&\times \rho_{\mathfrak{b},\mathfrak{c}}(-m,\kappa)\left(\frac1{2\pi i}\right)^5\int\cdots\int\frac{a^wb^wq^{2s}M^{u_1}N^{u_2}}{r^{2s+u_1+u_2} m^{\frac w2}}\hat{\ol{H}}_-(s,u_1,u_2,v,w,\kappa)dvdu_1du_2dwds d\kappa.\notag
\end{align}

Due to the ratio of gamma factors
\[
\frac{\Gamma(v)\Gamma(\frac12-s-u_2+\frac w2)}{\Gamma(\frac12-s-u_2+v+\frac w2)},
\]
Observing the gamma factors in $\hat{\ol{H}}_-$, we note that it has at best polynomial decay in $v$-aspect. So we should truncate all the integrals except the two over $v$ and $\kappa$ at height $(qK)^\ve$ with a negligible error. We split $\ol{\mathcal{E}}_{c-}(R,M^*;K)$ into two parts by writing the $v$-line of integral as
\[
\int_{(c_v)}=\int_{c_v-3i(qK)^\ve}^{c_v+3i(qK)^\ve}+\int_{c_v\pm3i(qK)^\ve}^{c_v\pm i\infty}.
\]
Since $v$ is small in the first part, a treatment identical to the $\ol{\mathcal{E}}_{c+}(R,M^*;K)$ case gives the same bound
\[
\ll q^{-\frac12+\ve}\left(\frac NM\right)^{\frac12}R^{-\ve}K^{-\ve}{M^*}^{-\ve}.
\]
In the treatment of the first part, the horizontal integral arising in our moving the $v$-line can be bounded identically to the second part case, which we show in detail next.

Note in the second part that moving $c_{u_1}$ to the right will not cross any poles. So we initially take the contours so that
\begin{align}\label{eqcontourDE1}
c_s=\ve,\ \ \ c_{u_1}=-c_{u_2}=c>1,\ \ \ c_v=\ve,\ \ \
 \left\{
\begin{aligned}
&c_w=2\ve, \ \ \ \ \ \ \ \text{for} \ \ M^*\le K^2Q,\\
&c_w=1+\ve,\ \ \ \text{for} \ \ M^*> K^2Q.
\end{aligned}
\right.
\end{align}
Then applying Stirling's approximation, we have
\begin{align}
\cos(\pi(s&+u_1+v+\tfrac w2))\Gamma(s+u_1+v+\tfrac w2+i\kappa)\Gamma(s+u_1+v+\tfrac w2-i\kappa)\notag\\
&\times \sin(\pi(s+u_1+v))\Gamma(\tfrac12-s-u_1-v)^2\frac{\Gamma(v)\Gamma(\frac12-s-u_2+\frac w2)}{\Gamma(\frac12-s-u_2+v+\frac w2)}\notag\\
\ll&(qK)^\ve e^{\pi|v'|}e^{-\frac\pi 2|v'-\kappa|}e^{-\frac\pi 2|v'+\kappa|}(1+|v|)^{-\frac32-c+c_s+\frac12 c_{w}}\notag
\end{align}
with $v'-v=\text{Im}(s+u_1+\frac w2)\ll (qK)^\ve$. A trivial estimate gives
\begin{align}\label{eqv-integral}
(qK)^\ve\int_0^\infty e^{\pi|v'|}e^{-\frac\pi 2|v'-\kappa|}e^{-\frac\pi 2|v'+\kappa|}(1+|v|)^{-\frac32-c+c_s+\frac12 c_{w}}dv\ll K^{-\frac12-c+c_s+\frac12 c_{w}}.
\end{align}


Applying \eqref{eqv-integral} and \eqref{eqboundss2} into the expression of $\ol{\mathcal{E}}_{c-}(R,M^*;K)$, we get the second part is bounded by
\begin{align}\label{eqolER}
\ll q^{-\frac12+2c_s+\ve}M^{c}N^{-c}K^{\frac12-c+c_s+\ve}\left(K+q^{-\frac12}R^{\frac12}\right)R^{-2c_s+\ve}{M^*}^{-\ve}.
\end{align}
With $c$ large, $c_s=\ve$ if $R\le K^2q$ and $c_s=\frac14+\ve$ if $R\le K^2q$, it is easy to see that \eqref{eqolER} is bounded by $q^{-\frac12+\ve}$ since $M\le N$.

Thus, we conclude that
\[
\ol{\mathcal{E}}_{c-}(R,M^*;K)\ll q^{-\frac12+\ve}\left(\frac NM\right)^{\frac12}R^{-\ve}K^{-\ve}{M^*}^{-\ve},
\]
which after summing over $K$ gives the bound  $\ol{\mathcal{E}}_{c-}(R,M^*)$ in \eqref{eqdcE} and then establishes the proposition.
\end{proof}



Similar calculations but with some minor variations give the same bound for the holomorphic spectrum $\ol{\mathcal{E}}_{h-}(R,M^*)$. The variations are due to different gamma factors, and we outline their treatments here. Applying Stirling's approximation to the gamma factors in $\tilde{\ol{H}}_-$, we have
\begin{align}\label{eqgammaH-}
&\frac{\Gamma(\frac{k-1}2+s+u_1+v+\tfrac w2)}{\Gamma(\frac{k+1}2-s-u_1-v-\tfrac w2)}\sin(\pi(s+u_1+v))\Gamma(\tfrac12-s-u_1-v)^2\frac{\Gamma(v) \Gamma(\frac12-s-u_2+\frac w2)}{\Gamma(\frac12-s-u_2+v+\frac w2)}\\
&\ \ \ \ \ \ll (K+|v|)^{-1+2c_s+2c_{u_1}+2c_v+c_w}(1+|v|)^{-\frac12-c_s-2c_{u_1}+c_{u_2}-2c_{v}-\frac12c_w}.\notag
\end{align}
We split $\ol{\mathcal{E}}_{h-}(R,M^*;K)$ into three parts by writing the $v$-line of integral as
\[
\int_{(c_v)}=\int_{c_v-3i(qK)^\ve}^{c_v+3i(qK)^\ve}+\int_{c_v\pm3i(qK)^\ve}^{c_v\pm iK^{1-\ve}}+\int_{c_v\pm iK^{1-\ve}}^{c_v\pm i\infty}.
\]
In the first part, $v$ is small, and we can bound it identically to $\ol{\mathcal{E}}_{c+}(R,M^*;K)$. For the third part, $|v|$ is large, and integrating \eqref{eqgammaH-} over $v$ provides the same bound as \eqref{eqv-integral}, so it follows identically to the second part of $\ol{\mathcal{E}}_{c-}(R,M^*;K)$. In the second part, we take $c_v=-A$ with $A>0$ large crossing no poles, and then \eqref{eqgammaH-} is not exceeding $O\left(K^{-2020}q^\ve\right)$. At last, the two horizontal integrals arising in our moving the $v$-line can be bounded identically to the first part and the third part respectively.

Also, the calculations identical to Proposition \ref{prodcE} but with some minor variations made in Proposition \ref{proM} give the following estimate for the Maass spectrum.
\begin{proposition}
We have
\begin{align}
\ol{\mathcal{E}}_{m\pm}(R,M^*)\ll q^{-\frac12+\theta+\ve}\left(\frac NM\right)^{\frac12}R^{-\ve}{M^*}^{-\ve}.
\end{align}
\end{proposition}

Thus, we conclude from the holomorphic spectrum, the Maass spectrum, and the Eisenstein spectrum that
\[
\ol{\mathcal{E}}_{\pm}(R,M^*)\ll q^{-\frac12+\theta+\ve}\left(\frac NM\right)^{\frac12}R^{-\ve}{M^*}^{-\ve},
\]
and we sum over $R$ and $M^*$ to get
\[
E_{\ol{M,N}}\ll q^{-\frac12+\theta+\ve}M^{-\frac12}N^{\frac12}q^\ve.
\]


\section{Off-diagonal terms with $M$ and $N$ far away}\label{secfar}
This section is devoted to bounding these $B_{M,N}$ as well as $B_{\ol{M,N}}$ with $M$ and $N$ far away from each other. We put our focus on $B_{M,N}$, and the treatment for $B_{\ol{M,N}}$ proceeds much the same after we classify them according to $M>N$ and $M<N$.
For notational convenience, we set all shifts to zero as before, and the arguments extend easily to nonzero parameter values. Also, we define the exponents $\mu,\nu$ as
\[
M=q^{\mu},\ \ \ \ N=q^{\nu}.
\]

Let $\eta=\frac1{80}-\frac 1{40}\theta$. Applying \eqref{eqEMN} as well as the trivial bound \eqref{eqtrivalbound}, we may reduce the proof of Theorem \ref{thmmains} to proving
\begin{align}
B_{M,N}\ll q^{-\eta+\ve}\notag
\end{align}
for
\begin{align}\label{eq1munu}
2-2\eta\le\mu+\nu\le 2,\ \ \ \ \ \  1-2\theta-2\eta\le\nu-\mu.
\end{align}

We recall that
\begin{align}
B_{M,N}=\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)\mathop{\sum\sum}_{\substack{m\equiv n(\bmod d)\\ (mn,q)=1}}\frac{d(m)d(n)}{m^{\frac12}n^{\frac12}}V\left(\frac{mn}{q^2}\right) W\left(\frac mM\right)W\left(\frac nN\right),\notag
\end{align}
which with the definition of $V$ becomes
\begin{align}
B_{M,N}=&\frac1{2\pi i}\int_{(\ve)}\left(\frac{q^2}{MN}\right)^{s}\frac{G(s)}{s}g(s)\notag\\
&\times\frac1{\vp^*(q)\sqrt{MN}}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)\mathop{\sum\sum}_{\substack{m\equiv n(\bmod d)\\ (mn,q)=1}}d(m)d(n) W_s\left(\frac mM\right)W_s\left(\frac nN\right),\notag
\end{align}
where $W_s(x)=x^{\frac12+s}W(x)$. Since $G(s)g(s)$ decays rapidly in the imaginary direction, we can omit the effect of $s$ to get
\begin{align}
B_{M,N}\ll\frac1{\vp^*(q)\sqrt{MN}}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)\mathop{\sum\sum}_{\substack{m\equiv n(\bmod d)\\ (mn,q)=1}}d(m)d(n) W\left(\frac mM\right)W\left(\frac nN\right).\notag
\end{align}

We write $d(n)=\sum\limits_{n_1n_2=n}1$ and apply the dyadic partition of unity to both $n_1$ and $n_2$ with $n_1\as N_1$, $n_2\as N_2$, and $N_1N_2\as N$. Without loss of generality, we assume $N_2\ge N_1$. Now the problem reduces to bounding $B_{M,N_1,N_2}$, defined as
\begin{align}
B_{M,N_1,N_2}=&\frac1{\vp^*(q)\sqrt{MN}}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right)\notag\\
&\times\mathop{\sum\sum}_{\substack{m\equiv n_1n_2(\bmod d)\\ (mn_1n_2,q)=1}}d(m) W\left(\frac mM\right)W\left(\frac {n_1}{N_1}\right)W\left(\frac {n_2}{N_2}\right)W\left(\frac {n_1n_2}{N_1N_2}\right).\notag
\end{align}
We define two exponents $\nu_1,\nu_2$ as
\[
N_1=q^{\nu_1},\ \ \ \ \ \ N_2=q^{\nu_2}
\]
with
\begin{align}\label{eqnu1nu2u}
\nu=\nu_1+\nu_2,\ \ \ \ \ \nu_1\le \nu_2.
\end{align}


After applying Mellin transform to $W\left(\frac{n_1n_2}{N_1N_2}\right)$, we can remove it as $V$, then
\begin{align}
B_{M,N_1,N_2}\ll\frac1{\sqrt{MN}}\sum_{(m,q)=1}d(m)\Delta_{N_1,N_2}(m)W\left(\frac mM\right)\notag
\end{align}
with
\begin{align}
\Delta_{N_1,N_2}(m)=\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right) \sum_{(n_1,q)=1}\sum_{\substack{n_2\equiv m\ol{n}_1(\bmod d)\\ (n_2,q)=1}}W\left(\frac {n_1}{N_1}\right) W\left(\frac {n_2}{N_2}\right).\notag
\end{align}
Since $(mn_1,q)=1$, we rewrite the condition in the sum over $n_2$ as
\begin{align}
\begin{array}{c}
n_2\equiv m\ol{n}_1 \pmod d\\
(n_2,q)=1
\end{array}\Longleftrightarrow
\begin{array}{c}
n_2\equiv m\ol{n}_1 \pmod d\\
(n_2,q_d)=1
\end{array},\notag
\end{align}
and remove the coprime condition $(n_2,q_d)=1$ by M\"{o}bius inversion, then
\begin{align}
\Delta_{N_1,N_2}(m)=&\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right) \sum_{a\mid q_d}\mu(a)\sum_{(n_1,q)=1} \sum_{\substack{n_2\equiv m\ol{n}_1(\bmod d)\\ n_2\equiv0(\bmod a)}}  W\left(\frac {n_1}{N_1}\right)W\left(\frac {n_2}{N_2}\right).\notag
\end{align}
By the Chinese Remainder Theorem, we apply Poisson summation formula to $n_2$, then
\begin{align}
\Delta_{N_1,N_2}(m)=&\frac1{\vp^*(q)}\sum_{d\mid q}\vp(d)\mu\left(\frac qd\right) \sum_{a\mid q_d}\mu(a)\sum_{(n_1,q)=1}\frac{N_2}{ad}\sum_{h} e\left(\frac{hx}{ad}\right) W\left(\frac {n_1}{N_1}\right)\wh{W}\left(\frac {h}{ad/N_2}\right),\notag
\end{align}
where $\wh{W}$ is the Fourier transform of $W$ and $x\pmod{ad}$ is given by
\begin{align}
\left\{
\begin{aligned}
&x\equiv 0 \pmod {a},\\
&x\equiv m\ol{n}_1 \pmod d.\notag
\end{aligned}
\right.
\end{align}
By Lemma \ref{lemmaspe(x)}, there is
\begin{align}
e\left(\frac{hx}{ad}\right)=e\left(\frac{hm\ol{a}\ol{n}_1}{d}\right),\notag
\end{align}
and thus
\begin{align}
\Delta_{N_1,N_2}(m)=&\frac{N_2}{\vp^*(q)}\sum_{d\mid q}\frac{\vp(d)}{d}\mu\left(\frac qd\right) \sum_{a\mid q_d}\frac{\mu(a)}{a}\sum_{(n_1,q)=1}\sum_{h} e\left(\frac{hm\ol{a}\ol{n}_1}{d}\right) W\left(\frac {n_1}{N_1}\right)\wh{W}\left(\frac {h}{ad/N_2}\right).\notag
\end{align}


To treat $\Delta_{N_1,N_2}(m)$, we first come to the $h=0$ term, which is equal to
\begin{align}
\frac{N_2}{\vp^*(q)}\wh{W}(0)\sum_{(n_1,q)=1}W\left(\frac{n_1}{N_1}\right)\sum_{d\mid q}\frac{\vp(d)}{d}\mu\left(\frac qd\right)\sum_{a\mid q_d}\frac{\mu(a)}{a},\notag
\end{align}
where the sum over $d$ and $a$ cancels adequately.
To see this, we alter the order of the summation, then
\begin{align}
\sum_{d\mid q}\frac{\vp(d)}{d}\mu\left(\frac qd\right)\sum_{a\mid q_d}\frac{\mu(a)}{a}=\sum_{a\mid q}\frac{\mu(a)}{a}\sum_{d\mid q_a}\frac{\vp(d)}{d}\mu\left(\frac qd\right).\notag
\end{align}
To the right-hand side, if there is a prime $p$ with $p\mid a$ and $p^2\mid q$, then $\mu(q/d)=0$ and it vanishes. If not the case, there is $q_a=\frac qa$ and it follows by Lemma \ref{lemsumd} that
\[
\sum_{d\mid q_a}\frac{\vp(d)}{d}\mu\left(\frac qd\right)=\mu(a)\sum_{d\mid \frac q{a}}\frac{\vp(d)}{d}\mu\left(\frac q{ad}\right)=a\frac{\mu(q)}q.
\]
Applying this into the sum over $a$, we find that it vanishes too.


We are remained to treat the terms with $h\neq0$. To be specific, we have
\begin{align}
B_{M,N_1,N_2}\ll\sum_{d\mid q}\frac{\vp(d)}d\mu^2\left(\frac qd\right) \sum_{a\mid q_d}\mu^2(a)|R(d,a)|,\notag
\end{align}
where
\begin{align}\label{eqR}
R(d,a)=\frac{N_2}{a\vp^*(q)\sqrt{MN}}\sum_{(m,q)=1}d(m)W\left(\frac mM\right)\sum_{(n_1,q)=1}\sum_{h\neq0} e\left(\frac{hm\ol{a}\ol{n}_1}{d}\right) W\left(\frac {n_1}{N_1}\right)\wh{W}\left(\frac {h}{H}\right)
\end{align}
with
\[
H=adN_2^{-1}\ll qN_2^{-1}.
\]
Due to the rapid decay of the Fourier transform $\wh{W}$, contribution from such terms with $|h|\ge Hq^{\ve}$ is negligible. So we may assume that the $h$-sum is over $h\le Hq^\ve$ with $Hq^\ve=adq^\ve N_2^{-1}\gg 1$, which gives particularly
\begin{align}\label{eqnu2}
\nu_2\le 1+\ve.
\end{align}
\begin{notation}
In the rest of this section, we should keep in mind that $a\mid q$ is square-free and $d\mid q$ satisfies $(a,d)=1$, and these also hold for $a_0,a_1,\dots$ and $d_0,d_1,\dots$, which always denote divisors of $a$ and $d$ respectively.
\end{notation}

The expression of $R(d,a)$ in \eqref{eqR} is not ready for directly processing right now, due to the coprime condition $(m,q)=1$ in the $m$-sum. The condition is perfect if we bound the $n$-sum by Weil bound. However, this bound is not sufficient for some cases, especially when $N_1$ is close to $q^{\frac12}$, for which we should also obtain savings from the $m$-sum.
By the reciprocity law, one may reduce the exponential function to $e\left(-\frac{\ol{d}hm}{an_1}\right)$, with which we may regard the $m$-sum as the generalized Estermann $D$-function. However, the coprime condition $(m,q)=1$ turns out to be an obstruction in applying Proposition \ref{progED1}. As we will see, removing the coprime condition, especially $(m,d)=1$ by M\"{o}bius inversion
\[
\sum_{(m,d)=1}a(m)=\sum_{d_0\mid d}\mu(d_0)\sum_{d_0\mid m}a(m),
\]
will bring an extra factor $d_0^{\frac12}$ to the final bound, which makes the bound invalid unless $d_0$ is small. To eliminate this influence, we should not pin our hope on the cancelation of $\mu(d_0)$. Instead, we observe that the common divisor $d_0=(m,d)$  could lower the bound for the $n$-sum. So the problem will be solved if we also remove the coprime condition $(m,d)=1$ when bounding the $n$-sum.

More precisely, we split the coprime condition in the sum over $m$ as
\begin{align}
(m,q)=1\Longleftrightarrow\left\{
\begin{aligned}
&(m,d)=1,\\
&(m,q_d)=1,\notag
\end{aligned}
\right.
\end{align}
and remove $(m,d)=1$ by M\"{o}bius inversion to get
\begin{align}
R(d,a)=\sum_{d_0\mid d}\mu(d_0)R_{d_0}(d,a),
\end{align}
where
\[
R_{d_0}(d,a)=\frac{N_2}{a\vp^*(q)\sqrt{MN}}\mathop{\sum_{(m,q_d)=1}}_{d_0\mid m}d(m)W\left(\frac mM\right)\sum_{(n_1,q)=1}\sum_{h\neq0} e\left(\frac{hm\ol{a}\ol{n}_1}{d}\right) W\left(\frac {n_1}{N_1}\right)\wh{W}\left(\frac {h}{H}\right).
\]
Now the problem reduces to proving
\begin{align}\label{eqmainRd0}
R_{d_0}(d,a)\ll q^{-\eta+\ve}.
\end{align}
Note that the $\ve$ in this bound allows our neglect of $\ve$ in following calculations of the exponents.

Before further evolution of $R_{d_0}(d,a)$, we first recall the region of exponents here. We sum up from \eqref{eq1munu}, \eqref{eqnu1nu2u}, and \eqref{eqnu2} that the region to treat is
\begin{align}\label{eqregion}
&2-2\eta\le \nu+\mu\le 2,\ \ \ \ \ \ \ \ \ \ \ \ 1-2\theta-2\eta\le\nu-\mu,\\
&\nu=\nu_1+\nu_2,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \tfrac12-\theta-2\eta\le\nu_1\le \nu_2\le 1,\notag
\end{align}
where $\frac12-\theta-2\eta$ in the last inequality is deduced from other inequalities.
Let $\rho$ be the exponent defined as $d_0=q^\rho$. We will bound $R_{d_0}(d,a)$ separately by partitioning this region into several parts, according to
\begin{enumerate}
  \item $1+2\eta-3\rho\le\nu-\mu$;
  \item $1-2\theta-2\eta\le\nu-\mu<1+2\eta-3\rho$:
  \begin{itemize}
           \item $\frac12+2\eta-\frac52\rho\le \nu_1$;
           \item $\frac12-\theta-2\eta< \nu_1< \frac12+2\eta-\frac52\rho$.
         \end{itemize}
\end{enumerate}
Since $\theta$ and $\eta$ are small, we should note that this partition works only when $\rho$ is very small. If not the case, the second inequality in \eqref{eqregion} dominates the bound in part (1), and thus part (2) disappears.


\subsection{Estimate for $\nu-\mu$ large}
In this section, we prove \eqref{eqmainRd0} with exponents in part (1), which turns out to be a direct result of Weil bound for Kloosterman sums. With
\[
d_1=d/d_0,
\]
it follows
\[
R_{d_0}(d,a)=\frac{N_2}{a\vp^*(q)\sqrt{MN}}\sum_{(m,q_d)=1}d(d_0m)W\left(\frac {d_0m}M\right)\sum_{(n_1,q)=1}\sum_{h\neq0} e\left(\frac{hm\ol{a}\ol{n}_1}{d_1}\right) W\left(\frac {n_1}{N_1}\right)\wh{W}\left(\frac {h}{H}\right).
\]
Furthermore, we may assume $(hm,d_1)=1$ at no cost when
apply Weil bound to the sum over $n_1$ and bound other sums trivially. Then we have
\begin{align}
R_{d_0}(d,a)\ll \frac{N_2q^{\ve}}{aq\sqrt{MN}}\frac{MH}{d_0}\left(d_1^{\frac12+\ve}+N_1d_1^{-1}\right),\notag
\end{align}
which after a simple calculation with $H=adN_2^{-1}$ becomes
\begin{align}
R_{d_0}(d,a)\ll d_0^{-\frac32}q^{\frac12+\ve}\left(\frac NM\right)^{-\frac12}+\left(\frac NM\right)^{-\frac12}\ll q^{-\eta+\ve}\notag
\end{align}
since $N_1,ad\le q^{1+\ve}$.





\subsection{Estimate for $\nu-\mu$ near $1$ and $\nu_1$ large}
This section is devoted to proving \eqref{eqmainRd0} with exponents in the first case of part (2), where $\frac NM$ is near $q$, and $N_1$ can not go far away from  $[q^{\frac12},q^{\frac34}]$. More precisely, an easy calculation shows
\[
\tfrac12+2\eta-\tfrac52\rho\le \nu_1< \tfrac34+\tfrac12\eta-\tfrac34\rho.
\]

We combine $m$ with $h$ to one longer variable $l=mh$. Then a trivial estimate shows
\begin{align}\label{eqRd0-DS}
R_{d_0}(d,a)\ll\frac{N_2q^\ve}{aq\sqrt{MN}}\sum_{l\le L}\left|\sum_{(n_1,q)=1} e\left(\frac{\ol{a}\ol{n}_1l}{d_1}\right) W\left(\frac {n_1}{N_1}\right)\right|
\end{align}
with
\[
d_1=d/d_0,\ \ \ \ \ \ \ \ L=MHq^\ve/d_0.
\]
This double sum is of independent interest for its many applications, and we will give a uniform upper bound for it in Proposition \ref{proDS} in the last section. As a direct consequence of this uniform bound, we have the following Lemma.

\begin{lemma}\label{lemDS}
Let $(\alpha_k)$ be a sequence of complex numbers satisfying $\alpha_k\ll k^\ve$. For any positive integers $L,K$, and $g$ satisfying $(g,d_1)=1$, we have
\begin{align}
\sum_{l\le L}\left|\sum_{\substack{k\le K\\ (k,d_1)=1}} \alpha_k e\left(\frac{gl\ol{k}}{d_1}\right) \right|\ll \left(LKd_1^{-\frac12}+LK^{\frac12}+L^{\frac12}K+L^{\frac12}d_1^{\frac34}\right)(LKd_1)^\ve.\notag
\end{align}
\end{lemma}

We apply Lemma \ref{lemDS} into \eqref{eqRd0-DS}  to get
\[
R_{d_0}(d,a)\ll\frac{N_2q^\ve}{aq\sqrt{MN}}\left(LN_1d_1^{-\frac12}+LN_1^{\frac12}+L^{\frac12}N_1+L^{\frac12}d_1^{\frac34}\right).
\]
An easy calculation with
\[
L\ll \frac{adM}{d_0N_2}q^\ve,\ \ \ \ ad\le q,\ \ \ \ d_1=\frac d{d_0}\ \ \ \  N_1N_2=N
\]
gives
\begin{align}\label{eqRd03terms}
R_{d_0}(d,a)\ll q^{-\frac12+\ve}\left(\frac MN\right)^{\frac12}N_1+q^\ve\left(\frac MN\right)^{\frac12}N_1^{\frac12}+q^{-\frac12+\ve}N_1^{\frac12}+d_0^{-\frac54}q^{\frac14+\ve}N_1^{-\frac12}.
\end{align}
For $\nu-\mu$ close to $1$ and
$\nu_1< \frac34+\frac12\eta-\frac34\rho$, all terms in \eqref{eqRd03terms} except the last one are obviously negligible, then
\begin{align}
R_{d_0}(d,a)&\ll q^{\frac14-\frac54\rho-\frac12\nu_1+\ve}.\notag
\end{align}
With
$\nu_1\ge\frac12+2\eta-\frac52\rho$,
we have
\[
\tfrac14-\tfrac54\rho-\tfrac12\nu_1\le \tfrac14-\tfrac54\rho-\tfrac12\times\left(\tfrac12+2\eta-\tfrac52\rho\right)= -\eta,
\]
proving the bound \eqref{eqmainRd0}.

\subsection{Estimate for $\nu-\mu$ near $1$ and $\nu_1$ near $\frac12$}
In this section, we prove \eqref{eqmainRd0} for the last case in part (2), where the exponents satisfy
\[
1-2\theta-2\eta\le\nu-\mu<1+2\eta-3\rho,\ \ \ \ \ \ \  \tfrac12-\theta-2\eta< \nu_1< \tfrac12+2\eta-\tfrac52\rho,
\]
with which $\frac NM$ and $N_1$ are close to $q$ and $q^{\frac 12}$ respectively.

We begin by recalling that
\begin{align}
R_{d_0}(d,a)=\frac{N_2}{a\vp^*(q)\sqrt{MN}}\sum_{\substack{(m,q_d)=1\\ d_0\mid m}}d(m)W\left(\frac mM\right)\sum_{(n_1,q)=1}\sum_{h\neq0} e\left(\frac{hm\ol{a}\ol{n}_1}{d}\right) W\left(\frac {n_1}{N_1}\right)\wh{W}\left(\frac {h}{H}\right).\notag
\end{align}
By Lemma \ref{lemmaspe(x)}, there is
\begin{align}
e\left(\frac{hm\ol{a}\ol{n}_1}{d}\right)=e\left(-\frac{\ol{d}hm}{an_1}\right) e\left(\frac{hm}{adn_1}\right).\notag
\end{align}
Since $h\ll ad N_2^{-1}q^\ve $, it is easy to see that
\[
\frac{hm}{adn_1}\ll\frac MNq^\ve\ll q^{-1+2\theta+2\eta+\ve},
\]
and thus
\[
e\left(\frac{hm}{adn_1}\right)=1+O\left(q^{-1+2\theta+2\eta+\ve}\right).
\]
A trivial estimate shows that the contribution to $R_{d_0}(d,a)$ from this $O$-term is negligible, then
\begin{align}
R_{d_0}(d,a)\ll&\frac{N_2q^\ve}{aq\sqrt{MN}}\sum_{{h\neq0}} \wh{W}\left(\frac {h}{H}\right)\sum_{(n_1,q)=1} W\left(\frac {n_1}{N_1}\right)\sum_{\substack{(m,q_d)=1\\ d_0\mid m}}d(m)e\left(-\frac{\ol{d}hm}{an_1}\right) W\left(\frac mM\right).\notag
\end{align}
Since the sum over $h$ will be bounded trivially with absolute values, we may assume $(h,a)=1$ in the above formula at no cost.

We will treat the sum over $m$ by writing it in terms of the generalized Estermann $D$-function, for which we should remove the extra coprime condition. After splitting the coprime condition into
\begin{align}
(m,q_d)=1\Longleftrightarrow\left\{
\begin{aligned}
&(m,a)=1,\\
&(m,q_{ad})=1,\notag
\end{aligned}
\right.
\end{align}
we may remove both conditions $d_0\mid m$ and $(m,q_{ad})=1$ by \eqref{eqdivisor} and \eqref{eqcoprime} respectively. Then we have
\begin{align}
R_{d_0}(d,a)\ll&\frac{N_2q^\ve}{aq\sqrt{MN}}\frac1{d_0}\sum_{d_{01}\mid d_0}\sum_{bc\mid q_{ad}}\frac{\mu(bc)}{bc}\notag\\
&\sum_{\substack{h\neq0\\ (h,a)=1}} \wh{W}\left(\frac {h}{H}\right)\sum_{(n_1,q)=1} W\left(\frac {n_1}{N_1}\right) \sum_{m}c_{bd_{01}}(m)d(m)\chi_0(m)e\left(-\frac{\ol{d}hm}{an_1}\right) W\left(\frac mM\right)\notag
\end{align}
with $\chi_0$ being the principal character modulo $a$.
Applying the exponential sum formula for $c_{bd_{01}}(m)$ and writing $g=(h,n_1)$, we get
\begin{align}\label{eqRd0R}
R_{d_0}(d,a)\ll \frac{N_2q^\ve}{aq\sqrt{MN}}\sum_{d_{01}\mid d_0}\sum_{b\mid q_{ad}}\frac{\mu^2(b)}{bd_0}\sum_{\substack{0\neq gh\le Hq^\ve\\ (h,a)=1}} \left|\sum_{(n_1,qh)=1} W\left(\frac {n_1g}{N_1}\right)U(h,n_1)\right|,
\end{align}
where
\begin{align}
U(h,n_1)=\mathop{\sum\nolimits^*}_{j(\bmod bd_{01})}\sum_{m}d(m)\chi_0(m) e\left(-\frac{\ol{d}hm}{an_1}\right) e\left(\frac{jm}{bd_{01}}\right)W\left(\frac mM\right).\notag
\end{align}



The treatment of $U(h,n_1)$ relies on properties of the generalized Estermann $D$-function.
Since $(an_1,bd_{01})=1$, we combine the two exponential functions by Lemma \ref{lemmaspe(x)} to get
\begin{align}
e\left(-\frac{\ol{d}hm}{an_1}\right) e\left(\frac{jm}{bd_{01}}\right)=e\left(\frac{my_j}{an_1bd_{01}}\right)\notag
\end{align}
with
\begin{align}
\left\{
\begin{aligned}
&y_j\equiv jan_1 \ \  \pmod {bd_{01}},\\
&y_j\equiv -b\ol{d_2}h \ \ \pmod {an_1},
\end{aligned}
\right.
\ \ \ \ \ \ \ \ \ \text{and}\ \ \ \ \ \ \ \ d_2=d/d_{01}.\notag
\end{align}
It is clear that $(y_j,an_1bd_{01})=1$, and thus after using the Mellin transform of $W$, we may write the sum over $m$ in terms of the generalized Estermann $D$-function, that is
\begin{align}
U(h,n_1)=\frac1{2\pi i}\mathop{\sum\nolimits^*}_{j(\bmod bd_{01})}\int_{(2)}M^sD\left(s,0,\frac{y_j}{an_1bd_{01}},\chi_0\right)\wt{W}(s)ds.\notag
\end{align}

Since $a$ is a square-free number satisfying $(a,n_1bd_{01})=1$ and $(y_j,an_1bd_{01})=1$, we may treat $D\left(s,0,\frac{y_j}{an_1bd_{01}},\chi_0\right)$ by using Proposition \ref{progED1}.
We move the integration to $-1$, and cross a double pole at $s=1$ with a residue not more than
\begin{align}
\ll\mathop{\sum\nolimits^*}_{j(\bmod bd_{01})}\frac{M}{an_1bd_{01}}q^\ve\ll MgN_1^{-1}q^\ve.\notag
\end{align}
So the contribution of the residue to $R_{d_0}(d,a)$ is a negligible error
\begin{align}
\ll \left(\frac{N}{M}\right)^{-\frac12}q^\ve\ll q^{-\frac12+\eta+\theta+\ve}.\notag
\end{align}
Then applying the functional equation of Proposition \ref{progED1} to the generalized Estermann $D$-function with $q=a$, $q_1=a_3$, $q_2=a_2$, $q_3=a_1$ gives
\begin{align}
U(h,n_1)=U_{+}(h,n_1)+U_{-}(h,n_1)+O(MgN_1^{-1}q^\ve),\notag
\end{align}
where
\begin{align}
U_{\pm}(h,n_1)=&\sum_{a_1a_2a_3=a}\frac{\mu(a_1a_2)}{a_1a_2}\notag\\
&\times\mathop{\sum\nolimits^*}_{j(\bmod bd_{01})}\mathop{\sum\nolimits^*}_{r(\bmod {a_2})}\frac1{2\pi i}\int_{(-1)}M^s (n_1\Omega)^{1-2s} D\left(1-s,0,\frac{\pm \ol{x}_{j,r}}{n_1\Omega}\right)\wt{W}_{\pm}(s)ds\notag
\end{align}
with
\[
\Omega=a_2a_3bd_{01}
\]
and $\ol{x}_{j,r}\pmod{n_1\Omega}$ given by
\begin{align}
\left\{
\begin{aligned}
&\ol{x}_{j,r}\equiv \ol{r} \ \ \ \  \ \ \ \ \ \ \ \ \ \ \ \pmod {a_2},\\
&\ol{x}_{j,r}\equiv -a_1\ol{b}d_2\ol{h} \ \ \ \pmod {a_3n_1},\\
&\ol{x}_{j,r}\equiv \ol{ja_2a_3n_1} \ \ \ \ \pmod {bd_{01}}.\notag
\end{aligned}
\right.
\end{align}
Here
\[
\wt{W}_{\pm}(s)=\wt{W}(s)2(2\pi)^{2s-2}\Gamma(1-s)^2C_{\pm}(\pi s),
\]
and $C_{+}(x)=1, C_{-}(x)=-\cos x$,


We expand $D$ into Dirichlet series with the variables change $s\rightarrow 1-s$, then
\begin{align}\label{eqU}
U_{\pm}(h,n_1)=&\sum_{a_1a_2a_3=a} \frac{\mu(a_1a_2)}{a_1a_2} \frac M{n_1\Omega}\mathop{\sum\nolimits^*}_{j(\bmod bd_{01})}\mathop{\sum\nolimits^*}_{r(\bmod {a_2})}\sum_{k}d(k)e\left(\frac{\pm \ol{x}_{j,r}k}{n_1\Omega}\right) W_0\left(\frac{kM}{n_1^2\Omega^2}\right).\notag
\end{align}
Here $W_0$ is defined as
\begin{align}
W_0(x)=\frac1{2\pi i}\int_{(\ve)}x^{-s}\wt{W}_{\pm}(1-s)ds,\notag
\end{align}
and by Stirling's approximation one can easily check its following properties: $x^jW_0^{(j)}(x)\ll_A x^{-A}$ for any $A>0$, and for $x$ small, there are $W_0(x)=c_1\log x+c_2+O(x^{1-\ve})$ and $W_0^{(j)}(x)=c_j+O(x^{1-\ve})$ for $j=1,2,\dots$, see also \cite[P.23]{You11} for details. Thus, the $k$-sum is actually over $k\ll \frac{N_1^2\Omega^2}{g^2M}q^\ve$.
We will only treat $U_{+}$ in detail just because the treatment of $U_{-}$ is identical.

We first execute both sums over $j$ and $r$ in $U_{+}(h,n_1)$ by writing them in terms of the Ramanujan sum. More precisely,  applying Lemma \ref{lemmaspe(x)} twice, we split the exponential function into three factors as
\begin{align}
e\left(\frac{\ol{x}_{j,r}k}{a_2a_3n_1bd_{01}}\right)= e\left(-\frac{\ol{x}y}{a_3n_1}\right) e\left(\frac{z_1k\ol{j}}{bd_{01}}\right)e\left(\frac{z_2k\ol{r}}{a_2}\right),\notag
\end{align}
where
\[
 x=a_2b^2d_{01}h,\ \ \ \  y=a_1d_2k, \ \ \ \ (z_1,bd_{01})=1,\ \ \ \ (z_2,a_2)=1.
\]
Exact values of $z_1$ and $z_2$ are inessential, since the last two exponents contribute a factor $c_{a_2bd_{01}}(k)$ after our summing over $j$ and $r$. Thus,
\begin{align}
U_{+}(h,n_1)=&\sum_{a_1a_2a_3=a} \frac{\mu(a_1a_2)}{a_1a_2} \frac M{n_1\Omega} \sum_{k}d(k)c_{a_2bd_{01}}(k) e\left(-\frac{\ol{x}y}{a_3n_1}\right) W_0\left(\frac{kM}{n_1^2\Omega^2}\right).\notag
\end{align}
Actually, we will bound the sum over $k$ trivially after summing over $n_1$, and the obvious bound $c_{a_2bd_{01}}(k)\ll (a_2bd_{01},k)^{1+\ve}$ implies that the factor $c_{a_2bd_{01}}(k)$ makes no essential contribution to the sum.

Now we conclude that
\begin{align}\label{eqRR+R-}
R_{d_0}(d,a)\ll R_{d_0}^{+}(d,a)+R_{d_0}^{-}(d,a)+\left( q^{-\frac12+\eta+\theta+\ve}\right),
\end{align}
where
\begin{align}\label{eqRlast}
R_{d_0}^{\pm}(d,a)=&\frac{MN_2q^\ve}{a qN_1\sqrt{MN}}\sum_{d_{01}\mid d_0}\sum_{b\mid q_{ad}}\frac{\mu^2(b)}{bd_0}\sum_{a_1a_2a_3=a} \frac{1}{a_1a_2\Omega}\\
&\times\sum_{\substack{0\neq gh\ll Hq^\ve\\ (h,a)=1}}g\sum_{(n_1,qh)=1} W\left(\frac {n_1g}{N_1}\right)\sum_{k}d(k)c_{a_2bd_{01}}(k) e\left(\mp\frac{\ol{x}y}{a_3n_1}\right) W_0\left(\frac{kM}{n_1^2\Omega^2}\right).\notag
\end{align}
Here we have denoted $\frac{N_1}{n_1g}W\left(\frac {n_1g}{N_1}\right)$ by $W\left(\frac {n_1g}{N_1}\right)$.

Let $V(h,k)$ be defined as
\begin{align}
V(h,k)=\sum_{(n_1,qh)=1} e\left(\frac{\ol{x}y}{a_3n_1}\right)  W\left(\frac {n_1g}{N_1}\right) W_0\left(\frac{kM}{n_1^2\Omega^2}\right).\notag
\end{align}
We will save from the sum over $n_1$ by using the theory of exponential sums and bound all other sums trivially. By the partial summation, it suffices to bound the following exponential sum
\begin{align}
\mathfrak{S}(x,y)=\sum_{\substack{(n_1,qh)=1\\ n_1\asymp N_1/g}} e\left(\frac{\ol{x}y}{a_3n_1}\right).\notag
\end{align}
By Lemma \ref{lemmaspe(x)},
\begin{align}
e\left(\frac{\ol{x}y}{a_3n_1}\right) =e\left(-\frac{\ol{a_3n_1}y}{x}\right) e\left(\frac{y}{a_3n_1x}\right).\notag
\end{align}
Recalling that $x=a_2b^2d_{01}h$, we note that any prime $p$ dividing $x$ also divides $qh$. So we split the coprime condition in $\mathfrak{S}(x,y)$ as
\begin{align}
(n_1,qh)=1\Longleftrightarrow\left\{
\begin{aligned}
&(n_1,x)=1,\\
&(n_1,q_x)=1,\notag
\end{aligned}
\right.
\end{align}
and we remove $(n_1,q_x)=1$ by M\"{o}bius inversion to get
\begin{align}
\mathfrak{S}(x,y)=\sum_{l\mid q_x}\mu(l)\mathop{\sum_{(n_1,x)=1,~l\mid n_1}}_{n_1\asymp N_1/g} e\left(-\frac{\ol{a_3n_1}y}{x}\right) e\left(\frac{y}{a_3n_1x}\right).\notag
\end{align}
Since $(l,x)=1$, we may write $n_1=l(j+xr)$ to split $n_1$ into residue classes$\pmod{x}$. Then we eliminate the first exponential function  to get
\begin{align}
\mathfrak{S}(x,y)\ll&q^\ve x\max_{l,j} \left|\sum_{\frac{N_1}{glx}<r\le \frac{N_1}{glx}+C} e\left(\frac{y}{a_3xl(j+xr)}\right)\right|\notag
\end{align}
with $C\ll \frac{N_1}{glx}$.

We bound the above exponential sum by the following lemma, which is Corollary 8.5 in \cite{IK04}.
\begin{lemma}[Corollary 8.5 of \cite{IK04}]\label{lemIK04}
Suppose that $f(u)$ satisfies
\begin{align}
\frac FA\le \frac{u^k}{k!}\left|f^{(k)}(u)\right|\le F\notag
\end{align}
for $k=2,3$ in the segment $[Q,2Q]$. Then for $1\le Q'\le Q\le F$ we have
\begin{align}
\sum_{Q<m\le Q+Q'}e(f(m))\ll AF^{\frac16}Q^{\frac12}\log 3Q,\notag
\end{align}
where the implied constant is absolute.
\end{lemma}

To apply the lemma, we let
\begin{align}
f(u)=\frac{y}{a_3xl(j+xu)},\ \ \ F\asymp\frac{yg}{a_3N_1x},\ \ \ Q=\frac{N_1}{glx},\notag
\end{align}
and $A$ be absolute. When $k$ is large enough to keep $F\ge Q$, we apply the lemma
to get
\begin{align}\label{eqbV}
V(h,k)\ll x\left(\frac{yg}{a_3N_1x}\right)^{\frac16} \left(\frac{N_1}{gx}\right)^{\frac12}q^\ve.
\end{align}
Since
\[x= a_2b^2d_{01}h,\ \ \ y=a_1d_2k,
\]
the bound in \eqref{eqbV} increases as $h$ and $k$ get large.
Thus, with
\[
h=\frac{H}{g}q^\ve= \frac{ad}{gN_2}q^\ve,\ \ \ \ \ \ \ k=\frac{N_1^2\Omega^2}{g^2M}q^\ve=\frac{(a_2a_3bd_{01}N_1)^2}{g^2M}q^\ve
\]
the largest possible values,
 a direct calculation gives
\begin{align}
F\asymp \frac{yg}{a_3N_1x}\ll \frac NM\frac{a_1d_2(a_2a_3bd_{01})^2}{a_3a_2b^2d_{01}ad}q^\ve=\frac NM q^\ve\notag
\end{align}
since
\[
a=a_1a_2a_3,\  \ \ \ \ \ \ \ d=d_{01}d_2.
\]
So the bound \eqref{eqbV} evolves into
\begin{align}
V(h,k)\ll \left(\frac NM\right)^{\frac16}N_1^{\frac12}(a_2b^2d_{01}H)^{\frac12}q^\ve.\notag
\end{align}
After summing over $g$, $h$, $k$ trivially, we have
\begin{align}
R_{d_0}^{\pm}(d,a)\ll&\frac{MN_2}{a qN_1\sqrt{MN}}\sum_{d_{01}\mid d_0}\sum_{b\mid q_{ad}}\frac{1}{bd_0}\sum_{a_1a_2a_3=a} \frac{1}{a_1a_2\Omega}H\frac{N_1^2\Omega^2}{M}\left(\frac NM\right)^{\frac16}N_1^{\frac12}(a_2b^2d_{01}H)^{\frac12}q^\ve.\notag
\end{align}
With $H= \frac{ad}{N_2}$, a direct calculation gives
\begin{align}
R_{d_0}^{\pm}(d,a)\ll&\frac{N_1^2}{aqN^{\frac56}M^{\frac23}}\sum_{d_{01}\mid d_0}\sum_{b\mid q_{ad}}\sum_{a_1a_2a_3=a} \frac {(a_2b^2d_{01})^{\frac12}(ad)^{\frac32}\Omega}{a_1a_2bd_0} q^\ve.\notag
\end{align}
Recalling $\Omega=a_2a_3bd_{01}$, we have
\begin{align}
R_{d_0}^{\pm}(d,a)\ll&\frac{N_1^2}{aqN^{\frac56}M^{\frac23}}\sum_{d_{01}\mid d_0}\sum_{b\mid q_{ad}}\sum_{a_1a_2a_3=a} \frac {a_2^{\frac12}a_3d_{01}^{\frac12}b(ad)^{\frac32}}{a_1} q^\ve\ll d_0^{\frac12}\frac{N_1^2q^{\frac12+\ve}} {N^{\frac56}M^{\frac23}}\notag
\end{align}
since $abd\le q$.



For small $k$, $F<Q$ is possible, in which case Lemma \ref{lemIK04} is not available. However, we can treat this special case by the trivial bound. The inequality $F<Q$ gives
\begin{align}
k<\frac{N_1^2a_3}{a_1d_2g^2}.\notag
\end{align}
Applying this into \eqref{eqRlast} and calculating all the sums trivially, we get
\begin{align}
R_{d_0}^{\pm}(d,a)\ll \frac{N_1^2}q\sqrt{\frac MN}q^\ve,\notag
\end{align}
which is certainly a negligible error with $\frac NM$ near $q$ and $N_1$ near $q^{\frac12}$.

In conclusion, we have that
\begin{align}
R_{d_0}^{\pm}(d,a)&\ll d_0^{\frac12}\frac{N_1^2q^{\frac12+\ve}} {N^{\frac56}M^{\frac23}}=d_0^{\frac12}q^{\frac12+\ve}(MN)^{-\frac34}\left(\frac MN\right)^{\frac1{12}}N_1^2=q^{\frac12+\frac12\rho-\frac34(\mu+\nu)-\frac1{12}(\nu-\mu)+2\nu_1+\ve}.\notag
\end{align}
With
\[
\mu+\nu\ge 2-2\eta,\ \ \ \ \ \ \ \nu-\mu>1-2\theta-2\eta,\ \ \ \ \ \ \ \nu_1< \tfrac12+2\eta-\tfrac52\rho,
\]
a simple calculation shows
\[
\tfrac12+\tfrac12\rho-\tfrac34(\mu+\nu)-\tfrac1{12}(\nu-\mu)+2\nu_1\le -\tfrac1{12}+\tfrac{17}3\eta+\tfrac16\theta-\tfrac92\rho\le -\eta
\]
as $\eta\le \frac1{80}-\frac 1{40}\theta$, which in combination with \eqref{eqRR+R-} proves \eqref{eqmainRd0}.

Thus, we complete the proof of Theorem \ref{thmmains}.




\section{A double sum in Kloosterman sums}\label{secdoublesum}
Let $g, q$ be positive integers satisfying $(g,q)=1$, and $(\alpha_k)$ be any sequence of complex numbers satisfying $\alpha_k\ll k^\ve$. We often encounter double sums of type
\[
\mathfrak{D}(L,K;q)=\sum_{l\le L}\left|\sum_{\substack{k\le K \\ (k,q)=1}} \alpha_k e\left(\frac{gl\ol{k}}{q}\right)\right|.
\]

For the special case with $\alpha_k=1$, applying Weil bound gives at once
\[
\mathfrak{D}(L,K;q)\ll \left(Lq^{\frac12}+LKq^{-1}\right)(LKq)^\ve.
\]
With $\alpha_k=1$ and $q$ prime, Friedlander and Iwaniec \cite[P. 337]{FI85} state a bound
\begin{align}\label{eqDFI85}
\mathfrak{D}(L,K;q)\ll q^{\frac18+\ve}(LK)^{\frac34}+q^{\frac14+\ve}L^{\frac54}K^{\frac14},
\end{align}
valid for $1\le L< K,~LK<q$. They also mention that this is also true with slight modification for composite $q$. However, the bound in \eqref{eqDFI85} is nontrivial for $L$ and $K$ in rather short segments (e.g., $L=q^{\frac18}$ and $K=q^{\frac38+\ve}$).

Also, with $\alpha_k=1$ and $q$ prime, Young \cite[Proposition 4.3]{You11} shows
\[
\mathfrak{D}(L,K;q)\ll L^{\frac12}q^{\frac34+\ve}+LK^{\frac12}q^\ve,
\]
valid for $L, K\ll q^{1+\ve}$.

Motivated by \cite{You11}, we prove a uniform bound for general $q$. For composite $q$, the treatment proceeds with a careful case analysis, based on exact expressions for Kloosterman sums to prime power moduli. These expressions have also been applied in the deduction of a large sieve type result for Kloosterman sums by Blomer and Mili\'cevi\'c, see also \cite[Theorem 5]{BM15}.
\begin{proposition}\label{proDS}
For any positive integers $q, L, K$, and any sequence of complex numbers $(\alpha_k)$ satisfying $a_k\ll k^\ve$, we have
\begin{align}
\mathfrak{D}(L,K;q)\ll
\left\{
\begin{aligned}
&\left(LKq^{-\frac12}+LK^{\frac12}+L^{\frac12}K+L^{\frac12}q^{\frac34}\right)(LKq)^\ve, \ \ \ \ \ \ \ K\gg q^{\frac12},\\
&\left(LK^{\frac12}+L^{\frac12}K^{\frac12}q^{\frac12}\right)(LKq)^\ve,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ K\ll q^{\frac12}.
\end{aligned}
\right.
\end{align}
\end{proposition}

The estimation of the double sums is related to counting integer solutions of a congruence equation
\[
\ol{k}_2-\ol{k}_1\equiv z \pmod q
\]
for any given integer $z$, which is of independent interest. In the special case with $q=p$ prime and $(z,p)=1$, it is equivalent to solving the equation
\[
k_1k_2\equiv z \pmod q.
\]
This special case has been widely studied, see also Chan and Shparlinski \cite{CS10}, Bourgain and Garaev \cite{BG14}, as well as Cilleruelo and Garaev \cite{CG11}. We obtain a uniform upper bound for the general case with any $q$ and $z$ in the following proposition.
\begin{proposition}\label{prok1k2}
For any integer $z$ and positive integers $K_1, K_2$, we have
\begin{align}\label{eqk1k2}
\mathop{\sum\nolimits^*}_{\ol{k}_2-\ol{k}_1\equiv z (\bmod q)} W\left(\frac {k_1}{K_1}\right)W\left(\frac {k_2}{K_2}\right)\ll \left(\frac{K_1K_2}{q}+\frac{K_1+K_2}{q^{\frac12}}+\frac{\max\{K_1,K_2\}(z,q)^{\frac12}}{q^{\frac12}}+q^{\frac12}\right)(K_1K_2q)^\ve.
\end{align}
\end{proposition}

\begin{remark}
Proposition \ref{prok1k2} is especially effective for $K_1,K_2\gg q^{\frac12}$.
If $\min\{K_1,K_2\}\ll q^{\frac12}$ and $\max\{K_1,K_2\}\ll q$, one should apply the trivial bound
\begin{align}\label{eqk1k2trivial}
\mathop{\sum\nolimits^*}_{\ol{k}_2-\ol{k}_1\equiv z (\bmod q)} W\left(\frac {k_1}{K_1}\right)W\left(\frac {k_2}{K_2}\right)\ll \min\{K_1,K_2\}.
\end{align}
\end{remark}
We first deduce Proposition \ref{proDS} from Proposition \ref{prok1k2}, and the proof of Proposition \ref{prok1k2} will be given in the next section.

\begin{proof}[Proof of Proposition \ref{proDS}] We apply the dyadic partition of unity to the sum over $l$ and use Cauchy inequality to get
\begin{align}
\mathfrak{D}(L,K;q)^2\ll L\mathop{\sum\nolimits^*}_{k_1,k_2\le K}\alpha_{k_1}\ol{\alpha}_{k_2}\sum_{l} e\left(\frac{gl(\ol{k}_1-\ol{k}_2)}{q}\right)W\left(\frac lL\right).
\end{align}
We split $l$ into residue classes$\pmod q$ and apply Poisson summation formula to get
\begin{align}
\mathfrak{D}(L,K;q)^2\ll&\frac{L^2}{q}\mathop{\sum\nolimits^*}_{k_1,k_2\le K}\alpha_{k_1}\ol{\alpha}_{k_2} \sum_z \sum_{c(\bmod q)} e\left(\frac{cg(\ol{k}_1-\ol{k}_2)}{q}\right)e\left(\frac{cz}{q}\right)\wh{W}\left(\frac {z}{q/L}\right).\notag
\end{align}
Note that the sum over $c$ is 0 unless
\begin{align}
g(\ol{k}_1-\ol{k}_2)+z\equiv0 \pmod {q},\notag
\end{align}
in which case it equals $q$. That is,
\begin{align}\label{eqD^2}
\mathfrak{D}(L,K;q)^2\ll L^2K^\ve\sum_z \wh{W}\left(\frac {z}{q/L}\right)\mathop{\sum\nolimits^*}_{\substack{k_1,k_2\le K\\ \ol{k}_2-\ol{k}_1\equiv \ol{g}z (\bmod q)}} 1.
\end{align}
Now the problem is reduced to counting integer solutions for the congruence equation $\ol{k}_2-\ol{k}_1=\ol{g}z \pmod q$, considered in Proposition \ref{prok1k2}.

For $K\gg q^{\frac12}$, we bound \eqref{eqD^2} by applying \eqref{eqk1k2}. We sum over $z$ with $0\le z\le q^{1+\ve}L^{-1}$ if $q^{1+\ve}L^{-1}(z,q)^{-1}\ge1$ and with the only term $z=0$ otherwise. More precisely,
\begin{align}
\mathfrak{D}(L,K;q)^2&\ll L^2\sum_z \wh{W}\left(\frac {z}{q/L}\right)\left(\frac{K^2}{q}+\frac{K}{q^{\frac12}}+\frac{K(z,q)^{\frac12}}{q^{\frac12}}+q^{\frac12}\right)(LKq)^\ve\notag\\
&\ll \left(\frac{L^2K^2}{q}+L^2K+LK^2+Lq^{\frac32}\right)(LKq)^\ve.\notag
\end{align}
For $K\ll q^{\frac12}$, the same calculation but with \eqref{eqk1k2trivial} in place of \eqref{eqk1k2} gives
\begin{align}
\mathfrak{D}(L,K;q)^2&\ll L^2K^{1+\ve}\sum_z \wh{W}\left(\frac {z}{q/L}\right)\ll \left(L^2K+LKq\right)(LKq)^\ve.\notag
\end{align}
Combining these two bounds proves the proposition.
\end{proof}

\subsection{Proof of Proposition \ref{prok1k2}}
The proof of Proposition \ref{prok1k2} relies on the estimate of a special exponential sum in the following.
\begin{lemma}\label{lemT}
For any integer $q\ge1$, and $x,y,z\in\mathbb{F}_q$,
we define $\mathcal{T}(x,y,z,q)$ by
\begin{align}
\mathcal{T}(x,y,z,q):=\mathop{\mathop{\sum\nolimits^*}_{a(\bmod q)}\mathop{\sum\nolimits^*}_{b(\bmod q)}}_{\ol{a}-\ol{b}\equiv z(\bmod q)} e\left(\frac{ax+by}{q}\right).\notag
\end{align}
Then we have
\begin{align}\label{eqT}
\mathcal{T}(x,y,z,q)\ll\left\{
\begin{aligned}
&(x,y,q)^{\frac12}q^{\frac12+\ve},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x=0~\text{or}~y=0,\\
&(x,y,q)(x+y,z,q/(x,y,q))^{\frac12}q^{\frac12+\ve},\ \ \ \ \ \ \text{otherwise}.
\end{aligned}
\right.
\end{align}
\end{lemma}

The proof of this lemma takes much space, so we complete the treatment of Proposition \ref{prok1k2} first and prove the lemma at last.
We apply Poisson summation formula in each of $k_1, k_2 \pmod{q}$ to see that
\begin{align}
\mathop{\sum\nolimits^*}_{\ol{k}_2-\ol{k}_1\equiv z (\bmod q)} W\left(\frac {k_1}{K_1}\right)W\left(\frac {k_2}{K_2}\right)=\frac{K_1K_2}{q^2} \sum_{x}\sum_y \mathcal{T}(x,y,z,q)\wh{W}\left(\frac {x}{q/K_1}\right)\wh{W}\left(\frac {y}{q/K_2}\right),\notag
\end{align}
by which Proposition \ref{prok1k2} follows from Lemma \ref{lemT} directly.
Note that a case analysis we should make in the summation over $x,y$. To be specific, the first term in \eqref{eqk1k2} is due to the case that both $x$ and $y$ equal $0$; the second term is due to the case only one of $x, y$ equals $0$; the third term is due to the case $x\neq0$, $y\neq 0$ and $[(x,y,q),(x+y,z,q/(x,y,q))]>q^{1+\ve}/\max\{K_1,K_2\}$, which means that one of $x,y$ is uniquely determined by the other; the last term is due to otherwise.


Now the rest space is devoted to the proof of the lemma.
\begin{proof}[Proof of Lemma \ref{lemT}]
By the Chinese Remainder Theorem and Lemma \ref{lemmaspe(x)},
there is
\begin{align}\label{eqmTgeneral}
\mathcal{T}(x,y,z,q)=\mathcal{T}(x\ol{q}_2,y\ol{q}_2,z,q_1)\mathcal{T}(x\ol{q}_1,y\ol{q}_1,z,q_2)
\end{align}
for any $q_1,q_2$ satisfying $q_1q_2=q$ and $(q_1,q_2)=1$. Thus, the problem reduces to the evaluation of $\mathcal{T}(x,y,z,q)$
for $q=p^m$, $m\ge1$.

When $q=p$, we treat $\mathcal{T}$ by writing it as a Ramanujan sum if $p\mid z$, and as a Kloosterman sum if $p\nmid z$.
For $p\mid z$, we calculate directly that
\begin{align}\label{eqpmidz}
\mathcal{T}(x,y,0,p)&=\mathop{\sum\nolimits^*}_{a(\bmod p)}e\left(\frac{a(x+y)}{p}\right)=c_{p}(x+y)\ll (x+y,z,p)^{1+\ve}.
\end{align}
For $p\nmid z$,
we solve for $b$ to have
\begin{align}
\mathcal{T}(x,y,z,p)=\mathop{\mathop{\sum\nolimits^*}_{a(\bmod p)}}_{a\neq\ol{z}} e\left(\frac{ax+\ol{(\ol{a}-z)}y}{p}\right).\notag
\end{align}
Using $\ol{(\ol{a}-z)}=a\ol{(1-za)}$ gives
\begin{align}
\mathcal{T}(x,y,z,p)=\mathop{\mathop{\sum\nolimits^*}_{a(\bmod p)}}_{a\neq\ol{z}} e\left(\frac{ax+a\ol{(1-az)}y}{p}\right).\notag
\end{align}
With the variables change $a\rightarrow \ol{z}-\ol{z}a$, it follows
\begin{align}
\mathcal{T}(x,y,z,p)=\mathop{\mathop{\sum\nolimits^*}_{a(\bmod p)}}_{a\neq1} e\left(\frac{-ax\ol{z}+\ol{a}y\ol{z}}{p}\right)e\left(\frac{x\ol{z}-y\ol{z}}{p}\right)= S(-x\ol{z},y\ol{z};p)e\left(\frac{x\ol{z}-y\ol{z}}{p}\right)-1.\notag
\end{align}
Thus, it follows from Weil bound that
\begin{align}\label{eqpnmidz}
\mathcal{T}(x,y,z,p)\ll (x,y,p)^{\frac12}p^{\frac12}.
\end{align}

When $q=p^m$ with $m\ge2$, we apply the orthogonal basis to rewrite $T$ as
\begin{align}\label{eqmT}
\mathcal{T}(x,y,z,q)&=\frac1q\sum_{t(\bmod q)}\mathop{\sum\nolimits^*}_{a(\bmod q)}\mathop{\sum\nolimits^*}_{b(\bmod q)}e\left(\frac{ax+by}{q}\right)e\left(\frac{(\ol{b}-\ol{a}+z)t}{q}\right)\\
&=\frac1qc_q(x)c_q(y)+\frac1qT(x,y,z,q),\notag
\end{align}
where
\[
T(x,y,z,q)=\sum_{t\nequiv0(\bmod q)}e\left(\frac{zt}{q}\right)S(-x,t;q)S(y,t;q).
\]
The evaluation of $T(x,y,z,q)$ will rely on exact expressions for Kloosterman sums to prime power moduli in situations. We will show full details for  $T(x,y,z,q)$ with $q=p^m$, $p>2$, and the special case for $p=2$ could be treated by the same arguments with minor variations. The variations are due to purely technical reasons, such as the solution of the congruence equation $2n\equiv a \pmod p^m$ is slightly different when $p=2$, which is not difficult to treat.

If $x\equiv0 \pmod q$ or $y\equiv0 \pmod q$, at least one of the Kloosterman sums in $T(x,y,z,q)$ reduces to the Ramanujan sum. Taking $x\equiv 0 \pmod q$ for example, with Weil bound for $S(y,t;q)$, we have
\begin{align}\label{eqT1}
T(x,y,z,q)\ll\sum_{t\nequiv0(\bmod q)}c_{q}(t)(y,q)^{\frac12}q^{\frac 12+\ve}\ll (x,y,q)^{\frac12}q^{\frac 32+\ve}.
\end{align}
The same bound holds for $y\equiv0 \pmod q$.


If neither of $x\pmod q$ and $y\pmod q$ is equal to $0$, Lemma \ref{lemKS0} means that the Kloosterman sums vanish unless $(x,q)=(y,q)=(t,q)=p^{m'}$,
in which case we apply Lemma \ref{lemKSrelation} to get
\[
S(-x,t;q)=p^{m'}S(-x',t';q'), \ \ \ \ \ \ \ S(y,t;q)=p^{m'}S(y',t';q')
\]
with $x'=x/p^{m'}$, $y'=y/p^{m'}$, $t'=t/p^{m'}$, and $q'=q/p^{m'}$. For the special case $m'=m-1$, applying Weil bound gives at once
\begin{align}\label{eqT2}
T(x,y,z,q)\ll\sum_{\substack{t\nequiv0(\bmod q)\\ p^{m-1}\mid t}}p^{2m-1}\ll p^{2m}\le(x,y,q)q^{\frac32}.
\end{align}

Now it remains to bound $T(x,y,z,q)$  for $(x,q)=(y,q)=(t,q)=p^{m'}$, $m'\le m-2$. Before doing this, we treat a special case $T'(x,y,z,q)$, defined as
\[
T'(x,y,z,q)=\mathop{\sum\nolimits^*}_{t(\bmod q)}e\left(\frac{zt}{q}\right)S(-x,t;q)S(y,t;q)
\]
with $(x,q)=(y,q)=1$ and $q=p^m\ge p^2$.
Then by Lemma \ref{lemKScoprime}, we have
\[
S(-x,t;q)S(y,t;q)=0
\]
unless that both $-t\ol{x}$ and $t\ol{y}$ are quadratic residues modulo $q$, in which case we write
\[
t\equiv -l^2x=l^2c_{x,y}^2y \pmod q
\]
with $c_{x,y}^2\equiv -x/y$.
Note that
\[
\mathop{\mathop{\sum\nolimits^*}_{t(\bmod q)}}_{t\equiv-l^2x}=\frac12\mathop{\sum\nolimits^*}_{l(\bmod q)}
\]
since every $t$ just appears twice in the sum over $l$. In addition, by \eqref{eqKSexact} we have
\begin{align}
S(-&x,t;q)S(y,t;q)=\left(\frac{-xyc_{x,y}}{q}\right)q\notag\\
& \times\left[\ve_q^2 e\left(\frac{2l(-x+c_{x,y}y)}{q}\right)+e\left(-\frac{2l(x+c_{x,y}y)}{q}\right) +e\left(\frac{2l(x+c_{x,y}y)}{q}\right)+\ve_q^2 e\left(\frac{2l(x-c_{x,y}y)}{q}\right)\right].\notag
\end{align}
So we conclude that
\[
T'(x,y,z,q)=T_1'+T_2'+T_3'+T_4',
\]
where
\[
T_1'=\left(\frac{-xyc_{x,y}}{q}\right)\frac q2 \ve_q^2 \mathop{\sum\nolimits^*}_{l(\bmod q)}e\left(\frac{-l^2xz+2l(-x+c_{x,y}y)}{q}\right),
\]
and there are similar expressions for $T_2'$, $T_3'$, and $T_4'$. We treat $T_1'$ only since the cases of other terms are similar.


To estimate $T_1'$, we have three cases to deal with:
\begin{itemize}
  \item $z\equiv0$ or $-x+c_{x,y}y\equiv0 \pmod q$;
  \item $q\neq(-x+c_{x,y}y,q)\neq(z,q)\neq q$;
  \item $(-x+c_{x,y}y,q)=(z,q)\neq q$.
\end{itemize}
If $z\equiv0 \pmod q$, the exponential sum reduces to the Ramanujan sum, and we have
\[
T_1'\ll qc_q(-x+c_{x,y}y)\ll (x+y,q)q^{1+\ve}=(x+y,z,q)q^{1+\ve},
\]
observing from the definition of $c_{x,y}$ that
\[
(-x+c_{x,y}y,q)\mid (x+y,q), \ \ \ \ \ \ \ \ (x+c_{x,y}y,q)\mid (x+y,q).
\]
If $-x+c_{x,y}y\equiv0 \pmod q$, the exponential sum reduces to a Gauss sum, and it follows by Lemma \ref{lemGauss}
\[
T_1'\ll q^{1+\ve}(z,q)\left(\frac{q}{(z,q)}\right)^{\frac12}\ll (z,q)^{\frac12}q^{\frac32+\ve}=(x+y,z,q)^{\frac12}q^{\frac32+\ve}.
\]

For the second case $q\neq(-x+c_{x,y}y,q)\neq(z,q)\neq q$, we find that $T'$ vanishes. To see this, we take $(-x+c_{x,y}y,q)<(z,q)$ for example, and the treatment of the other case proceeds identically. Let $(-x+c_{x,y}y,q)=p^{m_1}$ and $m=m_1+m_2$ with $m_2\ge 2$. After a reduction of the common factor in the exponential function, it follows
\begin{align}\label{eqT_1'}
T_1'=\left(\frac{-xyc_{x,y}}{q}\right)\frac {p^{m+m_1}}2\ve_q^2\mathop{\sum\nolimits^*}_{l(\bmod p^{m_2})}e\left(\frac{-l^2xz'+2l\omega_{x,y}}{p^{m_2}}\right)
\end{align}
with
$z'=z/p^{m_1}$ and $\omega_{x,y}=(-x+c_{x,y}y,q)/p^{m_1}$.
Put $l=dp^{m_2-1}+c$, where $c$ ranges $\bmod {p^{m_2-1}}$, $c$ is prime to $p$, and $d$ ranges freely $\bmod {p}$.  Then the exponential sum in $T_1'$ evolves into
\begin{align}
\mathop{\sum\nolimits^*}_{c(\bmod {p^{m_2-1}})}\sum_{d(\bmod {p})}e\left(\frac{-(d^2p^{2m_2-2}+2dcp^{m_2-1}+c^2)xz'+2c\omega_{x,y}}{p^{m_2}}\right) e\left(\frac{2d\omega_{x,y}}{p}\right).\notag
\end{align}
Since $(z',p)=p$ and $(\omega_{x,y},p)=1$, the sum over $d$ vanishes, which means
$T_1'=0$.

For $(-x+c_{x,y}y,q)=(z,q)=p^{m_1}\neq q$, we bound the sum over $l$ by writing it as a Gauss sum. With the same notations just as above, we also have the equation \eqref{eqT_1'}, but with $(z',p)=1$ and $(\omega_{x,y},p)=1$. To remove the coprime condition in the exponential sum, we add back all terms with $(l,q)>1$, which are regarded as several sums of type
\begin{align}\label{eqsuml0}
\mathop{\sum\nolimits^*}_{l_0(\bmod p^{m_2-m_0})}e\left(\frac{-l_0^2xz'p^{m_0}+2l_0\omega_{x,y}}{p^{m_2-m_0}}\right),
\end{align}
according to the value of $m_0$. A similar discussion as in the second case shows that the sum in \eqref{eqsuml0} is equal to $0$ unless $m_0=m_2$ or $m_0=m_2-1$, in which cases, the sum is equal to $1$ or $-1$.
Thus, the contribution of these added terms to $T_1'$ is not exceeding
\[
\ll p^{m+m_1}\ll(x+y,z,q)q.
 \]

Now with all $(l,q)>1$ terms added back, the new exponential sum over $l$ in \eqref{eqT_1'} turns out to be a Gauss sum, which is well controlled by Lemma \ref{lemGauss} that
\[
T_1'\ll p^{m+m_1+\frac {m_2}{2}+\ve}\ll(x+y,z,q)^{\frac12}q^{\frac32+\ve}.
\]
Thus, we conclude from the three cases that
\[
T'(x,y,z,q)\ll (x+y,z,q)^{\frac12}q^{\frac32+\ve}.
\]

To deduce $T(x,y,z,q)$ with $(x,q)=(y,q)=(t,q)=p^{m'}$ for $1\le m'<m$, we observe from Lemma \ref{lemKSrelation} that
\begin{align}
T(x,y,z,q)=p^{2m'}T'(x',y',z,q'),\notag
\end{align}
where $x'=x/p^{m'}$, $y'=y/p^{m'}$, and $q'=q/p^{m'}$, and thus it follows
\begin{align}\label{eqT3}
T(x,y,z,q)\ll (x,y,q)^{\frac12}(x+y,z,q/(x,y,q))^{\frac12}q^{\frac32+\ve}.
\end{align}

We apply \eqref{eqT1}, \eqref{eqT2}, and \eqref{eqT3} into \eqref{eqmT} to get that, for $q=p^m$ with $m\ge2$,
\begin{align}\label{eqmTp^m}
\mathcal{T}&(x,y,z,q)\ll\frac{(x,q)(y,q)}{q^{1-\ve}}+\left\{
\begin{aligned}
&(x,y,q)^{\frac12}q^{\frac12+\ve},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x\equiv0~\text{or}~y\equiv0,\\
&(x,y,q)q^{\frac12},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x\equiv y\equiv p^{m-1},\\
&(x,y,q)^{\frac12}(x+y,z,q/(x,y,q))^{\frac12}q^{\frac12+\ve},\ \ \ x\equiv y\equiv p^{m'}\le p^{m-2},\\
&0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{otherwise}.
\end{aligned}
\right.
\end{align}
For convenience, we use the uniform bound
\[
\ll(x,y,q)(x+y,z,q/(x,y,q))^{\frac12}q^{\frac12+\ve}
\]
 to control the last three cases in the second term of \eqref{eqmTp^m}. Consequently, the first term in \eqref{eqmTp^m} is dominated by the second term since $(x,q)(y,q)\le (x,y,q)q$. That is,
\begin{align}\label{eqmTp^m2}
\mathcal{T}(x,y,z,q)\ll\left\{
\begin{aligned}
&(x,y,q)^{\frac12}q^{\frac12+\ve},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x\equiv0~\text{or}~y\equiv0 \pmod q,\\
&(x,y,q)(x+y,z,q/(x,y,q))^{\frac12}q^{\frac12+\ve},\ \ \ \ \ \text{otherwise},
\end{aligned}
\right.
\end{align}
for $q=p^m$ with $m\ge2$.


Now we come to deduce an upper bound for $\mathcal{T}(x,y,z,q)$ for general $q$ by multiplying all distinct prime powers. By the Prime Number Theory, the number of different prime powers in $q$ is not more than $\log q/\log\log q$, so the product of all absolute constants in $\ll$ is bounded by $q^{\ve}$. Note that \eqref{eqpmidz} and \eqref{eqpnmidz} are dominated by \eqref{eqmTp^m2}, and the first bound in \eqref{eqmTp^m2} for $x\equiv0~\text{or}~y\equiv0 \pmod q$ is also dominated by the second one. So we multiply all prime powers with \eqref{eqmTp^m2} to get the same bound for general $q$, completing the proof.
\end{proof}









