You are on the right track. Suppose there is at least one root, $n$. As you have noted, there is a factorization $x^3 - a = (x-n)(x^2 + xn + n^2)$. Now, to see that the quadratic has a root, you need to show that the discriminant $(-3n^2)$ is a square modulo $p$, or equivalently that $-3$ is a square modulo $p$. And remember, you need to use the mod $3$ condition. Does quadratic reciprocity ring a bell?
Assume $N=pq$, with $p,q$ primes. Then $(\Bbb{Z}/N\Bbb{Z})^\times = (\Bbb{Z}/p\Bbb{Z})^\times \times (\Bbb{Z}/q\Bbb{Z})^\times$. In this setting, we have the four square roots of unity $\pm 1=\pm(1,1),\:\pm x=\pm (1,-1)$.
I think your algorithm is a bit unclear. What I think you mean is the following. Let $g\in (\Bbb{Z}/N\Bbb{Z})^\times$.
Compute $g^{a^i r}$ for $i=1,\dots,a$ until you find $i_0$ such that $g^{a^{i_0} r}=1$.
If $g^{a^{i_0-1} r}$ is not $\pm 1$, you have found $\pm x$, and this gives you a factorization of $N$.
The claim is now that the procedure above, for $\text{gcd}(g,N)=1$, finds $\pm x$ with probability at least $\frac12$.
The paper you link claim it is a "straightforward argument", but my proof is slightly complex, unfortunately. Let me provide it anyway.
Let $g=(a,b)\in (\Bbb{Z}/p\Bbb{Z})^\times \times (\Bbb{Z}/q\Bbb{Z})^\times$. It should be clear that $g^{2^i r} = -1$ if and only if $c^{2^i r} = d^{2^i r} =-1$, with the same $i$. In particular, they are both roots in the polynomial $x^{2^i r} -1$ in their respective rings.
Pick $s,t$ such that $p-1=2^e s$ and $q=2^f t$ for some $e,f\geq 1$. Note that $s=\text{gcd}(r,p-1)$, $t=\text{gcd}(r,q-1)$. Then $x^{2^i r}-1$ has at most $2^is$ roots in $\Bbb{Z}/p\Bbb{Z}$ and $2^it$ roots in $\Bbb{Z}/q\Bbb{Z}$. This is because the image of $x^{s}$ in $\Bbb{Z}/p\Bbb{Z}$ are just the elements of even multiplicative order, and further taking this image to some odd power simply permutes it. The same story holds for $q$.
Let us wlog assume $e\leq f$. We can now count all the elements $g$ such that $g^{2^i r}= -1$ for some $i$. By the above, we only need to consider $i\leq e-1$, since any higher $i$ will not allow $c^{2^i r} = -1$. Since $c$ and $d$ must both be roots of the same polynomial, in their respective fields, we find an upper bound on the number of elements satisfying $g^{2^i r} = -1$ as $(2^i s)(2^i t)$. Summing over $i$, we get
$$\sum_{i=0}^{e-1} (2^i s)(2^i t) = st\sum_{i=0}^{e-1} 4^i = st \frac{4^e-1}{3} = \frac13(2^e s)(2^e t)-\frac{st}{3} = \frac{1}{3\cdot 2^{f-e}}pq-\frac{st}{3}.$$
One should then compare it to $\frac12\phi(N) = \frac12 pq-\frac12(p+q)$, and one will see that the latter dominates the former (with a little work), thus it is true that less than half the elements in $(\Bbb{Z}/N\Bbb{Z})^\times$ hit $-1$ under the procedure above, and thus, more than half find $\pm x$.
Best Answer
The answer is yes, because Chebotarev theorem tells us the asymptotic $$\#\{ p \le N, x^3-a\bmod p \text{ has a root }\}\sim \frac{|H|}{|G|} \frac{N}{\log N}$$ where $G=Gal(\Bbb{Q}(a^{1/3},\zeta_3)/\Bbb{Q}), |G|=6$ and $H$ is the set of $\sigma \in G$ such that $\sigma(a^{1/3})=a^{1/3}$, ie. $|H| = 2$ independently of $a$.