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