Denote $g(n)=\sum_{a,b=1}^n\gcd(a,b)$, can we prove that
$$g(n)=\frac6{\pi^2}n^2\ln n+Cn^2+O(n\ln n)$$, where $C=-\frac12+\frac{6}{\pi^2}(-\frac12+\gamma-\ln(2\pi)+12\ln A),$ where $\gamma$ denotes Euler's constant and $A$ denotes Glaisher's constant?

**Failed Attempt**

$$g(n)=\sum_{k=1}^n\varphi(k)\left[\frac nk\right]^2\\
=\sum_{k=1}^n\varphi(k)\left[\frac nk\right]^2=\sum_{k=1}^n\varphi(k)\left(\frac nk+O(1)\right)^2\\=n^2\sum_{k=1}^n\frac{\varphi(k)}{k^2}+O\left(n\sum_{k=1}^n\frac{\varphi(k)}k\right)\\
=\frac{6}{\pi^2}n^2\ln n+O(n^2)$$
**Closer Attempt to $O(n^{3/2}\ln n)$**

$$g(x)=2\sum_{k\le n\le x}\gcd(n,k)-\frac12x^2+O(x)\\
=2\sum_{dl\le x}d\varphi(l)-\frac12x^2+O(x)\\
=2\sum_{d\le \sqrt x}d\sum_{l\le x/d}\varphi(l)+
2\sum_{l\le \sqrt x}\varphi(l)\sum_{d\le x/l}d-2\left(\sum_{d\le \sqrt x}\varphi(d)\right)\left(\sum_{l\le \sqrt x}l\right)-\frac12x^2+O(x)\\
=2\sum_{d\le \sqrt x}d\left(\frac3{\pi^2}\frac{x^2}{d^2}+O\left(\frac xd\ln x\right)\right)+2\sum_{l\le \sqrt x}\varphi(l)\frac{x^2}{2l^2}+O\left(\frac xl\right)-2\left(\frac 3{\pi^2}x+O(\sqrt x\ln x)\right)\left(\frac x2+O(\sqrt x)\right)\\
=\frac{6}{\pi^2}x^2(\ln\sqrt x+\gamma)+\frac{6}{\pi^2}x^2(\ln\sqrt x+\gamma)-\frac{36}{\pi^4}\zeta'(2)x^2-\left(\frac 3{\pi^2}+\frac12\right)x^2+O(x^{3/2}\ln x)\\
=\frac6{\pi^2}x^2\ln x+Cx^2+O(x^{3/2}\ln x)$$
How can I improve the big-$O$ term?

The asymptotic you want does not hold just because the "last-term fluctuation" $$ g(n)-g(n-1) = 2\sum_{a=1}^n \gcd(a,n)-n $$ is too large. Indeed, denoting the sum in the right-hand side by $\sigma(n)$, we have $$ \sigma(n) =\sum_{d\mid n} d\varphi(n/d) = n \sum_{d\mid n} \prod_{p\mid(n/d)} \Big(1-\frac1p\Big) = n \sum_{d\mid n} \prod_{p\mid d} \Big(1-\frac1p\Big). $$

If now $n=p_1\dotsb p_k$ is the product of the first $k$ primes, then $$ \sigma(n) = \Big(2-\frac1{p_1}\Big)\dotsb\Big(2-\frac1{p_k}\Big)\,n \ge (3/2)^k n $$ while, using the prime number theorem, it is not difficult to see that $k\ge c\ln n/\ln\ln n$ with an absolute constant $c>0$; hence, $$ \sigma(n) \ge \exp(c'\ln n/\ln\ln n)\,n,\ n=p_1\dotsb p_k. $$ On the other hand, for $n$ prime we have $$ \sigma(n) = 2n-1. $$ It follows that one cannot get an asymptotic for $g(n)$ with the remainder term better than $O(\exp(c'\ln n/\ln\ln n)\,n)$.