*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})$.

For the torus, all you need to do is to periodize the heat kernel: if we write $$\tilde p_t(x) = \sum_{n \in \mathbb{Z}^d} p_t(x + n),$$ then the solution on the torus is given by $$u(t, x) = u_0 * \tilde p_t(x),$$ where the convolution on the torus is defined as $$ \int_{[0,1)^d} u_0(y) \tilde p_t(x - y) dy .$$

You can find a number of references for the $\mathbb{R}^d$ case in my two survey papers: a more probabilistic view in Section 4 of

M. Kwaśnicki, *Fractional Laplace Operator and its Properties*, in: A. Kochubei, Y. Luchko, *Handbook of Fractional Calculus with Applications. Volume 1: Basic Theory*, De Gruyter Reference, De Gruyter, Berlin, 2019

or an analytical perspective in Section 2.6 of

M. Kwaśnicki, *Ten equivalent definitions of the fractional Laplace operator*, Frac. Calc. Appl. Anal. 20(1) (2017): 7–51

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)$.