We provide sufficient conditions for the existence of viscosity solutions of fractional semilinear elliptic PDEs of index $α\in (1,2)$ with polynomial gradient nonlinearities on $d$-dimensional balls, $d\geq 2$. Our approach uses a tree-based probabilistic representation of solutions and their partial derivatives using $α$-stable branching processes, and allows us to take into account gradient nonlinearities not covered by deterministic finite difference methods so far. In comparison with the existing literature on the regularity of solutions, no polynomial order condition is imposed on gradient nonlinearities. Numerical illustrations demonstrate the accuracy of the method in dimension $d=10$, solving a challenge encountered with the use of deterministic finite difference methods in high-dimensional settings.