  We provide sufficient conditions for the existence of viscosity solutions of
fractional semilinear elliptic PDEs of index $\alpha \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 $\alpha$-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.
