One way to think about this sort of problem is to embed in continuous time. Take $N$ independent Poisson processes of rate 1. (Think of $N$ independent Geiger counters, each going off at rate 1, if you like). A point in the $i$th process corresponds to picking the $i$th ball. Since the processes are independent and all have the same rate, the sequence of ball selections is just a sequence of independent uniform choices, as we desire.

Let $M_i(x)$ be the number of points in the $i$th Poisson process up to time $x$.
Then the distribution of $M_i(x)$ is Poisson($x$). In particular,
$P(M_i(x)\geq 2)= 1-(1+x)e^{-x}$. The time of the first point in such a process has exponential(1) distribution, so its probability density function is $e^{-x}$.

So fix one ball, say ball 1. Consider the event that when ball 1 is first chosen, all the other $N-1$ balls have each been chosen at least twice. To get the probability of this event, integrate over the
time that ball 1 is first chosen (i.e. the time of the first event in process 1):

$\int_0^\infty e^{-x} P\big(M_i(x)\geq 2 \text{ for } i=2,3,\dots,N\big) dx$

$=\int_0^\infty e^{-x} \big(1-(1+x)e^{-x}\big)^{N-1} dx$.

The same applies for any ball, and the events are disjoint, so an exact answer to your question is

$N\int_0^\infty e^{-x} \big(1-(1+x)e^{-x}\big)^{N-1} dx$.

I don't know if it's possible to get an exact expression for this integral, but it's easy enough to bound it. For any $K$

$N\int_0^\infty e^{-x} \big(1-(1+x)e^{-x}\big)^{N-1} dx$

$\leq N \int_0^K e^{-x} \big(1-e^{-x}\big)^{N-1} dx
+ N\int_K^\infty e^{-x} \big(1-Ke^{-x}\big)^{N-1} dx$.

Take $K=\frac12 \log N$, for example; it's easy to evaluate both integrals
exactly (substitute $u=e^{-x}$) and to show that they both tend to 0 as $N\to\infty$.