Given $N\gg0$, small $\epsilon>0$ and $\alpha\in(\frac12,\frac34)$ how many pairs $$(a,b)\in[-N^{\alpha+\epsilon},N^{\alpha+\epsilon}]\times[-N^{\alpha+\epsilon},N^{\alpha+\epsilon}]$$ of integers of size $N^{\alpha+\epsilon}$ have $$gcd(|a|,|b|)>N^{1/3}?$$
Probability that $d|a$ and $d|b$ is $\frac1{d^2\zeta(2)}$ and so probability that every $d$ that divides $a$ and that divides $b$ is $\leq N^{1/3}$ is $$\int_{1}^{N^{1/3}}\frac1{x^2\zeta(2)}dx=\frac1{\zeta(2)}\Bigg(1-\frac1{N^{1/3}}\Bigg)$$ and so probability that there is at least one divisor $>N^{1/3}$ is $$1-\frac1{\zeta(2)}\Bigg(1-\frac1{N^{1/3}}\Bigg)=1-\frac6{\pi^2}+\frac1{\zeta(2)N^{1/3}}\underbrace{\approx}_{N\rightarrow\infty}0.3920728981459733713367232207\dots.$$
Is this correct estimation?
No, your estimation is not correct. The integral you are computing is the probability that some $d \le N^{1/3}$ divides $a$ and $b$, since you are adding the probabilities for each $d \le N^{1/3}$.
What we want is the probability that all $d$ which divide $a$ and $b$ are $\le N^{1/3}$. Conversely, we want no $d > N^{1/3}$ dividing $a$ and $b$, which means we will multiply the probabilities of each individual $d > N^{1/3}$ dividing $a$ and $b$. This means the probability we want is: $$ \prod_{d = N^{1/3}}^{N^{\alpha + \epsilon}} \left( 1 - \frac{1}{\zeta(2)d^2} \right) $$ Let this product be called $S$. Then: $$ \ln S = \sum_{d = N^{1/3}}^{N^{\alpha + \epsilon}} \ln \left( 1 - \frac{1}{\zeta(2)d^2} \right) $$ Since the function inside the sum is monotone increasing, it is asymptotically equal to: $$ \int_{N^{1/3}}^{N^{\alpha + \epsilon} + 1} \ln \left( 1 - \frac{1}{\zeta(2)x^2} \right) \mathrm{d}x $$ To make the formulas simpler, let $\beta = N^{\alpha + \epsilon} + 1$ and $\eta = N^{1/3}$. Then the integral above evaluates too: $$ \beta \ln \left( 1 - \frac{1}{\zeta(2)\beta^2} \right) - \eta \ln \left( 1 - \frac{1}{\zeta(2)\eta^2} \right) + \frac{2}{\sqrt{\zeta(2)}} \tanh^{-1}\left(\sqrt{\zeta(2)}\beta\right) - \frac{2}{\sqrt{\zeta(2)}} \tanh^{-1} \left(\sqrt{\zeta(2)} \eta \right) $$ This is asymptotically equal to $\ln S$, so $S$ is asymptotically: $$ \frac{\left(1 - \frac{1}{\zeta(2)\beta^2}\right)^\beta}{\left(1 - \frac{1}{\zeta(2)\eta^2}\right)^\eta} \cdot \left( \frac{1 - \zeta(2)\beta^2}{1 - \zeta(2)\eta^2} \right)^{1/\sqrt{\zeta(2)}} $$ This is (asymptotically) the probability that $gcd(a,b) \le N^{1/3}$. Therefore, the probability we actually want is one minus this. However, as $N \to \infty$, the above goes to $\infty$ (as far as I am able to determine), and so the probability that $gcd(a,b) > N^{1/3}$ appears to approach $0$.