Using the method of contour integration, one can get

$$P(N)=\frac{\log N}{\zeta(2)N}+\frac{A}{N}+O(N^{-23/22+o(1)}),$$

where $P(N)$ is our probability and $A=\frac{3\gamma}{\zeta(2)}-1-\frac{2\zeta'(2)}{\zeta(2)^2}-\frac{1}{\zeta(2)}$. To prove this formula, lets note that if

$$N^2P(N)-(N-1)^2P(N-1)=a(N)$$

then $a(N)-1$ is equal to the number of pairs of the form $(x,N)$ or $(N,x)$ with $xN=y^2$ for some $y$, $x<N$. One can easily see that if $D(N)^2$ is a maximal square divisor of $N$ then $a(N)-1=2D(N)-2$. This is because $N=MD(N)^2$ with squarefree $M$ and if $xN=y^2$, then we should have $x=Mz^2$ for some integer $z$. There are exactly $D(N)-1$ admissible values for $z$, hence there are $2(D(N)-1)$ pairs of the form $(x,N)$ or $(N,x)$ with $xN=y^2$ and $x<N$. So we get

$$a(N)=2D(N)-1.$$

Now let us consider the generating Dirichlet series of the function $a(N)$:

$$f(s)=\sum_{n\in \mathbb N} \frac{a(n)}{n^s}=2\sum_{n \in \mathbb N} D(n)/n^s-\zeta(s).$$

The function $D(n)$ is multiplicative and for any prime $p$ we have $D(p^k)=p^{[k/2]}$. Therefore,

$$\sum_{n \in \mathbb N} D(n)n^{-s}=\prod_p (1+p^{-s}+p^{1-2s}+p^{1-3s}+\ldots)=\prod_p(1+p^{-s})(1+p^{1-s}+p^{2-2s}+\ldots)=\frac{\zeta(s)\zeta(2s-1)}{\zeta(2s)}.$$

Next, using the truncated version of Perron's formula and the fact that $D(n) \leq \sqrt n$, we obtain for any $b>1$ (and fixed) and any $T$ the following:

$$\sum_{n \leq N} D(n)=\frac{1}{2\pi i}\int_{b-iT}^{b+iT} N^s\frac{\zeta(s)\zeta(2s-1)}{s\zeta(2s)}ds+O\left(\frac{N^{3/2}\log N}{T}\right).$$

Moving the contour of integration to the line $\mathrm{Re}\,s=1/2+\varepsilon$ for any $\varepsilon>0$ and using the estimates $\zeta(s) \ll |s|^{1/6}$, $\zeta(2s-1) \ll |s|^{1/2}$ and $\zeta(s)^{-1} \ll 1$ one can prove that

$$\sum_{n \leq N} D(n)=\mathrm{Res}_{s=1}N^s\frac{\zeta(s)\zeta(2s-1)}{s\zeta(2s)}+O(N^\varepsilon(\sqrt N T^{5/6}+NT^{-1/3}+N^{3/2}T^{-1})).$$

Choosing $T=N^{6/11}$, we finally obtain

$$\sum_{n \leq N} D(n)=\mathrm{Res}_{s=1}N^s\frac{\zeta(s)\zeta(2s-1)}{s\zeta(2s)}+O(N^{21/22+\varepsilon}).$$

It remains to compute the residue.
For $s \to 1$ we have

$$\frac{N^s}{s}=N+(N\log N-N)(s-1)+O((s-1)^2),$$

$$\zeta(s)=\frac{1}{s-1}+\gamma+O(s-1),$$

$$\zeta(2s-1)=\frac{1}{2(s-1)}+\gamma+O(s-1)$$

and

$$\frac{1}{\zeta(2s)}=\frac{1}{\zeta(2)}-\frac{2\zeta'(2)}{\zeta(2)^2}(s-1)+O((s-1)^2).$$

Multipying this gives the result

$$\frac{N}{2\zeta(2)^2(s-1)^2}+\frac{N\log N/2\zeta(2)+N(A+1)/2}{s-1}+O(1),$$

so

$$\sum_{n \leq N} D(n)=\frac{N\log N}{2\zeta(2)}+N(A+1)/2+O(N^{21/22+o(1)})$$

and

$$N^2P(N)=\sum_{n \leq N} a(N)=2\sum_{n \leq N} D(n)-N=\frac{N\log N}{\zeta (2)}+AN+O(N^{21/22+o(1)})$$

which is equivalent to the stated result.

Analogous estimates for $k$-th powers with $k>1$ will give us something like

$$\frac{2\zeta(k-1)}{\zeta(k)N}-\frac{1}{N}+O(N^{-1-e_k})$$

with $e_k>0$, as the corresponding generating function is $2\frac{\zeta(s)\zeta(ks-1)}{\zeta(ks)}-\zeta(s)$. Also, I think that the error term could be significantly improved under the assumption of the Riemann Hypothesis.