  We study stable solutions to the equation $(-\Delta)^{1/2} u = f(u)$, posed
in a bounded domain of $\mathbb{R}^n$. For nonnegative convex nonlinearities,
we prove that stable solutions are smooth in dimensions $n\leq 4$. This result,
which was known only for $n=1$, follows from a new interior H\"older estimate
that is completely independent of the nonlinearity $f$.
  A main ingredient in our proof is a new geometric form of the stability
condition. It is still unknown for other fractions of the Laplacian and,
surprisingly, it requires convexity of the nonlinearity. From it, we deduce
higher order Sobolev estimates that allow us to extend the techniques developed
by Cabr\'e, Figalli, Ros-Oton, and Serra for the Laplacian. In this way we
obtain, besides the H\"older bound for $n\leq 4$, a universal $H^{1/2}$
estimate in all dimensions.
  Our $L^\infty$ bound is expected to hold for $n\leq 8$, but this has been
settled only in the radial case or when $f(u) = \lambda e^u$. For other
fractions of the Laplacian, the expected optimal dimension for boundedness of
stable solutions has been reached only when $f(u) = \lambda e^u$, even in the
radial case.
