$$\frac{\pi^2}{6} = \zeta(2) = \sum_{k=1}^\infty \frac{1}{k^2}$$ I hope we agree.
Now how do I get a grip on the tail end $\sum_{k \geq N} \frac{1}{k^2}$ which is the tail end which goes to zero?
I want to show that $\sqrt{x}\cdot \mathrm{tailend}$ is bounded as $x \to \infty$.
All this to show that $x\cdot \mathrm{tailend} = \mathcal O\sqrt{x}$
The purpose is to get the asymptotic formula for the distribution of square free integers.
p.269 Exercise 8 Stewart and Tall.
You can show that
$$\int_0^{\infty} dt \frac{t \, e^{-N t}}{e^t-1} = \sum_{k=N+1}^{\infty} \frac1{k^2}$$
For large $N$, apply integration by parts to get an expansion in $1/N$:
$$\begin{align} \int_0^{\infty} dt \frac{t \, e^{-N t}}{e^t-1} &= \frac1{N} + \frac1{N} \int_0^{\infty} dt \, e^{-N t} \frac{e^t (1-t)-1}{\left(e^t-1\right)^2} \\ &= \frac1{N} - \frac1{2 N^2} + \frac1{N^2} \int_0^{\infty} dt \, e^{-N t} \frac{e^t \left(e^t (t-2)+t+2\right)}{\left(e^t-1\right)^3}\\\sum_{k=N+1}^{\infty} \frac1{k^2} &= \frac1{N} - \frac1{2 N^2} + \frac1{6 N^3} + O\left (\frac1{N^4}\right )\end{align}$$
Generally,
$$\sum_{k=1}^{N} \frac1{k^2} = \frac{\pi^2}{6} - \sum_{k=1}^{K} \frac{a_k}{N^k} + O\left (\frac1{N^{K+1}}\right )$$
where
$$a_k = \left [\frac{d^{k-1}}{dt^{k-1}} \frac{t}{e^t-1} \right ]_{t=0} = B_{k-1}$$
i.e., the $k-1$st Bernoulli number.
As an example, if we compare the remainder sum as a function of $N$ with the three-term sum above in a log-log plot as follows:
we see that the slope is $-5$, meaning that the next nonzero term in the series is $O(1/N^5)$, not $1/N^4$ as one might expect. The reason for this is that $B_3=0$. Note also the noise for very high values of $N$, as we have run against some machine precision for this calculation.
We can actually compute the complete asymptotic expansion of the remainder term. Introduce $$S(x) = \sum_{k\ge 1} \left(\frac{1}{k^2}-\frac{1}{(x+k)^2}\right)$$ so that our answer is given by $$\frac{\pi^2}{6}-S(N).$$
Re-write $S(x)$ as follows: $$S(x) = \sum_{k\ge 1} \frac{1}{k^2} \left(1-\frac{1}{(x/k+1)^2}\right).$$ The sum term is harmonic and may be evaluated by inverting its Mellin transform.
Recall the harmonic sum identity $$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) = \left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$ where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have $$\lambda_k = \frac{1}{k^2}, \quad \mu_k = \frac{1}{k} \quad \text{and} \quad g(x) = 1-\frac{1}{(x+1)^2}.$$ We need the Mellin transform $h^*(s)$ of $h(x)=g(x)-1$ which is $$\int_0^\infty -\frac{1}{(x+1)^2} x^{s-1} dx = \left[\frac{1}{1+x} x^{s-1} \right]_0^\infty - (s-1) \int_0^\infty \frac{1}{1+x} x^{s-2} dx.$$
The fundamental strips of these three components are $\langle 0, 2\rangle$, followed by $\langle 1,2 \rangle$ and finally, $\langle 1,2\rangle$.
The contribution in the square brackets disappears in the fundamental strip $\langle 1, 2\rangle.$
Now this last Mellin transform can be evaluated using a keyhole contour with the slot on the positive real axis, which gives $$(1-e^{2\pi i (s-1)})\times \int_0^\infty \frac{1}{1+x} x^{s-1} dx = 2\pi i \times \mathrm{Res}\left(\frac{1}{1+x} x^{s-1}; x= -1\right)$$ which yields in turn $$\int_0^\infty \frac{1}{1+x} x^{s-1} dx = 2\pi i \frac{ e^{\pi i (s-1)}}{1-e^{2\pi i s}} = - 2\pi i \frac{1}{e^{-\pi i s}-e^{\pi i s}} = \frac{\pi}{\sin(\pi s)}.$$ It follows that the Mellin transform of $g(x)$ is $$g^*(s) = - (s-1) \frac{\pi}{\sin(\pi (s-1))} = (s-1)\frac{\pi}{\sin(\pi s)}$$ with fundamental strip $\langle -1, 0\rangle$ as can be seen by expanding $g(x)$ about zero and about infinity and testing convergence of the integral.
Therefore the transform $Q(s)$ of $S(x)$ is $$ Q(s) = (s-1)\frac{\pi}{\sin(\pi s)} \zeta(2-s) \quad\text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \sum_{k\ge 1}\frac{1}{k^2} k^s = \zeta(2-s).$$ The half plane of convergence of the zeta term is $\Re(s)<1.$
We thus obtain the Mellin inversion integral $$ S(x) = \frac{1}{2\pi i} \int_{-1/2-i\infty}^{-1/2+i\infty} Q(s)/x^s ds$$ which we evaluate by shifting it to the right for an expansion about infinity.
We get $$\mathrm{Res}(Q(s)/x^s; s = 0) = -\frac{\pi^2}{6}$$ and $$\mathrm{Res}(Q(s)/x^s; s = 1) = \frac{1}{x}$$ and $$\mathrm{Res}(Q(s)/x^s; s = 2) = -\frac{1}{2x^2}.$$ For the remaining poles at $s=q$ where $q>2$ we obtain $$\mathrm{Res}(Q(s)/x^s; s = q) = (q-1) (-1)^q \zeta(2-q)\frac{1}{x^q} \\= (-1)^q (q-1) \zeta(-(q-2))\frac{1}{x^q} = (-1)^{q-1} (q-1)\frac{B_{q-1}}{q-1} \frac{1}{x^q} = (-1)^{q-1} B_{q-1}\frac{1}{x^q} .$$ Now this only contributes when $q$ is odd so that we may simplify it to $B_{q-1}\frac{1}{x^q} .$ (Note that those zero values from the Bernoulli numbers at odd indices cancel the poles from the sine term, so it is in fact correct to let the summation for the expansion range over all residues at poles at $q>2.$) This gives for the asymptotic expansion $$S(x) \sim \frac{\pi^2}{6} - \frac{1}{x} + \frac{1}{2x^2} - \sum_{q>2} \frac{B_{q-1}}{x^q}.$$
We finally have $$\frac{\pi^2}{6} - S(N) \sim \frac{1}{N} - \frac{1}{2N^2} + \sum_{q>2} \frac{B_{q-1}}{N^q}$$ as discovered by Ron Gordon in his excellent answer.
The tail is $$\frac{1}{N^2}+\frac{1}{(N+1)^2}+\frac{1}{(N+2)^2}+\cdots.$$ For $N\gt 1$ this is less than $$\frac{1}{(N-1)(N)}+ \frac{1}{(N)(N+1)}+ \frac{1}{(N+1)(N+2)}+\cdots.$$ Note that $\frac{1}{t(t+1)}=\frac{1}{t}-\frac{1}{t+1}$. Using this we find that the tail is less than $$\left(\frac{1}{N-1}-\frac{1}{N}\right)+ \left(\frac{1}{N}-\frac{1}{N+1}\right)+ \left(\frac{1}{N+1}-\frac{1}{N+2}\right)+\cdots.$$ Open up the brackets and note the massive cancellation (telescoping). We find that the tail is less than $\frac{1}{N-1}$.
Another way: Draw a picture of the curve $y=\frac{1}{x^2}$. Note that the area of the rectangle of height $\frac{1}{n^2}$ that goes from $x=n-1$ to $x=n$ is $\frac{1}{n^2}$.
This rectangle is entirely below the curve on the interval from $n-1$ to $n$. Thus the tail is less than $$\int_{N-1}^\infty \frac{dx}{x^2}.$$ Evaluate the integral. We get $\frac{1}{N-1}$.
Use either the "Euler-Maclaurin summation formula" or "Abel summation formula" applied to the function $f(x) = 1/x^2$.
Maybe this should be a comment, but it is a bit better to use $$\frac{\pi^2}{12}=\eta(2)=\sum_{n\geqslant 1}(-1)^{n-1}\frac 1 {n^2}$$ and the usual bound for alternating series in this case.
Since the question includes a hint that the goal is to study the distribution of square-free integers, here is an observation.
The Euler product for the Dirichlet series $$L(s) = \sum_{n\ge 1} \frac{q_n}{n^s}$$ with $q_n$ the indicator function of positive integers that are square-free is given by $$L(s) = \prod_p \left(1+\frac{1}{p^s}\right) = \prod_p \frac{1-1/p^{2s}}{1-1/p^s} = \frac{\zeta(s)}{\zeta(2s)}.$$ It follows by the Wiener-Ikehara theorem that $$\sum_{k=1}^{N} q_k \sim \mathrm{Res}(L(s); s=1) \frac{N^1}{1} = \frac{1}{\zeta(2)} N = \frac{6}{\pi^2} N.$$
The actual question is to show that
$$ xD(x) = x \prod_{p_j}\left(1-\frac{1}{1-1/p_j}\right) = x\sum_{k \geq 1} \frac{1}{k^2} = \frac{6x}{\pi^2} + O(\sqrt x)$$
I get $o(\sqrt x)$.
$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\on}[1]{\operatorname{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} &\bbox[5px,#ffd]{\sum_{k\ \geq N}\,{1 \over k^{2}}} = \overbrace{\sum_{k = 1}^{\infty}{1 \over k^{2}}} ^{\ds{\zeta\pars{2}}}\ -\ \sum_{k = 1}^{N - 1}{1 \over k^{\color{red}{2}}} \\[5mm] = &\ \zeta\pars{2} - \bracks{\zeta\pars{\color{red}{2}} + {\pars{N - 1}^{1 - \color{red}{2}} \over 1 - \color{red}{2}} + \color{red}{2}\int_{N - 1}^{\infty}\,\,{\braces{x} \over x^{\color{red}{2} + 1}}\,\dd x} \end{align}
The last term is an identity related to the $\ds{\zeta}$-function.
Then, \begin{align} &\bbox[5px,#ffd]{\sum_{k\ \geq N}\,{1 \over k^{2}}} = {1 \over N - 1} - 2\int_{N - 1}^{\infty}\,\,{\braces{x} \over x^{3}}\,\dd x \end{align} Note that $\ds{0 < 2\int_{N - 1}^{\infty}\,\,{\braces{x} \over x^{3}}\,\dd x < 2\int_{N - 1}^{\infty}\,\,{\dd x \over x^{3}} = {1 \over \pars{N - 1}^{2}}}$ which yields $$ \bbx{{1 \over N - 1} - {1 \over \pars{N - 1}^{2}} < \bbox[5px,#ffd]{\sum_{k\ \geq N}\,{1 \over k^{2}}} < {1 \over N - 1}} \\ $$
