EDIT: The following is (mostly) for $\alpha < 2$; scroll to the bottom for more on general $\alpha$. Kudos to Abdelmalek Abdesselam (again).
As for $\mathbb{R}^d$, this is classical: the solution is given by the "fractional heat kernel": $$u(t,x) = u_0 * p_t(x),$$ where $p_t(x)$ is the inverse Fourier transform of $\exp(-t |\xi|^\alpha)$. Since $p_t$ is a probability density function, convolution with $p_t$ is a contraction on every $L^p(\mathbb{R}^d)$. The kernel $p_t$ has several alternative representations; in particular, we have Bochner's subordination formula, which asserts that $p_t(x)$ is a mixture of Gaussians: $$p_t(x) = \int_0^\infty q_s(x) \eta_t(s) ds,$$ where $$q_s(x) = (4 \pi s)^{-d/2} \exp(-|x|^2 / (4t))$$ is the Gaussian and $\eta_t(s)$ is a probability density distribution with Laplace transform $\exp(-t \xi^{\alpha/2})$.
I do not know any reference specifically for the torus. You may search for papers on the fractional heat equation on manifolds (or even "fractals"/"$d$-sets"/"metric measure spaces"), but this will likely be too general and abstract for your needs. (I only know of a paper Fractional Laplacian on the torus by Luz Roncal and Pablo Raúl Stinga, but this one is about the extension technique, not very useful for the heat equation.)
EDIT: what about general $\alpha > 0$?
For general $\alpha > 0$, in $\mathbb{R}^d$, a solution is again given by the convolution with $p_t$ given as an inverse Fourier transform of $\exp(-t |\xi|^\alpha)$. This is no longer a positive function if $\alpha > 2$, but it is anyway an integrable function. Here is a short (but perhaps not the most elementary) proof of this fact.
If $\alpha$ is an even integer, then the Fourier transform of $p_t$ is a Schwartz class function, and hence $p_t$ is Schwartz class. If $\alpha$ is not an even integer, then $p_t$ is still smooth, but it no longer decays rapidly. Now the result of:
K. Soni, R.P. Soni, Slowly Varying Functions and Asymptotic Behavior of a Class of Integral Transforms I, II, III. J. Anal. Appl. 49 (1975): 166--179; 477--495; 612--628
applied to the $d$-dimensional Hankel transform provides an asymptotic expansion of $p_t(x)$ at infinity, which in particular implies that $p_t$ is of constant sign in a neighbourhood of infinity. This easily leads to the conclusion that $p_t$ is integrable: otherwise, its Fourier transform would necessarily diverge at zero.
The convolution with $p_t$ is therefore again a bounded operator on $L^p(\mathbb{R}^d)$ for every $p \in [1, \infty]$. With no doubt it is written somewhere, but I do not have a reference at hand.
Of course, periodization leads to similar results on the torus. However, for general $\alpha > 0$ it is way easier to simply note that $\tilde{p}_t$ is given by a Fourier series with rapidly decreasing coefficients, and hence it is infinitely smooth. For this reason, the convolution with $\tilde p_t$ is a bounded operator on $L^p(\mathbb{T}^d)$.