We study gradient estimates of q-harmonic functions u of the fractional
Schr{\"o}dinger operator Δα/2+q, α∈(0,1] in bounded
domains D⊂Rd. For nonnegative u we show that if q is H{\"o}lder
continuous of order η>1−α then ∇u(x) exists for any x∈D and |\nabla u(x)| \le c u(x)/ (\dist(x,\partial D) \wedge 1). The
exponent 1−α is critical i.e. when q is only 1−α H{\"o}lder
continuous ∇u(x) may not exist. The above gradient estimates are well
known for α∈(1,2] under the assumption that q belongs to the Kato
class \calJ^{\alpha - 1}. The case α∈(0,1] is different. To obtain
results for α∈(0,1] we use probabilistic methods. As a corollary, we
obtain for α∈(0,1) that a weak solution of Δα/2u+qu=0 is in fact a strong solution