We present a construction of the parametrix for the fractional Helmholtz equation, $[(-\Delta)^s - \tau^{2s} r(x)^{2s} + q(x)] u = 0$, making use of geometrical optics solutions. We show that the associated eikonal equation is the same as in the classical case, while in the first transport equation the effect of nonlocality is visible only in the zero-th order term, which depends on s. Moreover, we show that the approximate geometrical optics solutions present different behaviors in the regimes 0= 1/2. While the latter case is rather similar to the classical one (s = 1), in the former case we find that the potential is a strong perturbation, which changes the propagation of singularities. As an application, we study the inverse problem consisting in recovering the potential q from Cauchy data when the refraction index r is fixed and simple. Using our parametrix, we prove that Holder stability holds for this problem. This is an improvement over the state of the art for fractional wave equations, for which the usual Runge approximation argument can provide only logarithmic stability. The introduction of the fractional Helmholtz equation is motivated by the following: The solutions of the (static) Dirichlet problem for the Laplacian in a domain with a polyhedral boundary present so-called corner singularities; when the domain is an infinite wedge, a polyhedral cone, or a polygon in the two-dimensional case, the solutions exhibit a similar behavior. Such singularities are mitigated in the fractional Laplacian case with s > 1. This study is further motivated by applications in nonlocal or fractional gradient elasticity models emerging from geophysics.