We give sharp two-sided estimates of the semigroup generated by the fractional Laplacian plus the Hardy potential on R d , including the case of the critical constant. We use Davies' method back-to-back with a new method of integral analysis of Duhamel's formula.The subject of the paper can be tracked down to Baras and Goldstein [3], who proved the existence of nontrivial nonnegative solutions of the classical heat equation ∂ t = ∆ + κ|x| −2 in R d for 0 ≤ κ ≤ (d − 2) 2 /4, and nonexistence of such solutions, that is explosion, for bigger constants κ. Vazquez and Zuazua [36] studied the Cauchy problem and spectral properties of the operator in bounded subsets of R d using the improved Hardy-Poincaré inequality, and they used weighted Hardy-Poincaré inequality in the more delicate case of the whole of R d . Sharp upper and lower bounds for the heat kernel of the Schrödinger operator ∆ + κ|x| −2 were obtained by Liskevich and Sobol [26, p. 365, Example 3.8, 4.4 and 4.10] for 0 < κ < (d − 2) 2 /4. Milman and Semenov proved the upper and lower bounds for κ ≤ (d − 2) 2 /4, see [28, Theorem 1], [29] and [22, Section 10.4]. Note that Moschini and Tesei [30, Theorem 3.10] gave estimates for the subcritical case in bounded domains, and Filippas, Moschini and Tertikas [18] obtained the critical case in bounded domains. Recently Metafune, Sobajima and Spina [27] extended the results of Milman and Semenov to sigular gradient-and-Schrödinger perturbations of ∆. Because of the borderline singularity of the function R dx → κ|x| −2 at the origin, the choice of κ influences the growth rate of the heat kernel at the origin.The operators ∆+κ|x| −2 play distinctive roles in limiting and self-similar phenomena in probability [32] and partial differential equations [33]. This results in part from the scaling of the corresponding heat kernel, which is similar to that of the Gauss-Weierstrass kernel. Analogous applications are expected for the Hardy perturbations of the fractional Laplacian, see [34] for first attempts. We also note that the effect of Schrödinger perturbations of ∆ is much milder if κ|x| −2 is replaced by functions in appropriate Kato classes. We refer to Bogdan and Szczypkowski [15, Section 1, 4] for references and results on Gaussian bounds for Schrödinger heat kernels along with an approach based on the so-called 4G inequality. The case when even the Gaussian constants do not deteriorate is discussed by Bogdan, Dziubański and Szczypkowski [11]. Estimates of Schrödinger perturbations of general operators and their transition semigroups are given by Bogdan, Hansen and Jakubowski in [12], with focus on situations with 3G inequality. Bogdan, Jakubowski and Sydor [14] and Bogdan, Butko and Szczypkowski [8] estimate Schrödinger perturbations of integral kernels which are not necessarily semigroups. The results of these paper give comparability or near comparability of the perturbed kernel with the original one. In fact, the explicit estimate in [14, Theorem 3] suffices for many applications.The operator (1.1) cannot be ...