In this paper, we investigate the large time behavior of the solution to the fractional heat equation $$\begin{aligned} \frac{\partial u}{\partial t}(t,x)=-(-\Delta )^{\beta /2}u(t,x)+u(t,x)\frac{\partial ^{d+1} W}{\partial t\partial x_1\cdots \partial x_d},\quad t>0,\quad x\in \mathbb {R}^d, \end{aligned}$$ ∂ u ∂ t ( t , x ) = - ( - Δ ) β / 2 u ( t , x ) + u ( t , x ) ∂ d + 1 W ∂ t ∂ x 1 ⋯ ∂ x d , t > 0 , x ∈ R d , where $$\beta \in (0,2)$$ β ∈ ( 0 , 2 ) and the noise W(t, x) is a fractional Brownian sheet with indexes $$H_0, H_1,\ldots ,H_d\in (\frac{1}{2},1)$$ H 0 , H 1 , … , H d ∈ ( 1 2 , 1 ) . By using large deviation techniques and variational method, we find a constant $$M_1$$ M 1 such that for any integer $$p\ge 1$$ p ≥ 1 and $$\alpha _0\beta +\alpha <\beta ,$$ α 0 β + α < β , $$\begin{aligned} \lim _{t\rightarrow \infty }t^{-\frac{2\beta -\beta \alpha _0-\alpha }{\beta -\alpha }} \log Eu(t,x)^p=p^{\frac{2\beta -\alpha }{\beta -\alpha }}\left( \frac{\alpha _H}{2} \right) ^{\frac{\beta }{\beta -\alpha }}M_1, \end{aligned}$$ lim t → ∞ t - 2 β - β α 0 - α β - α log E u ( t , x ) p = p 2 β - α β - α α H 2 β β - α M 1 , where $$\alpha _0=2-2H_0$$ α 0 = 2 - 2 H 0 , $$\alpha =\sum \nolimits _{j=1}^d(2-2H_j)$$ α = ∑ j = 1 d ( 2 - 2 H j ) and $$\alpha _H=\prod \nolimits _{i=0}^dH_i(2H_i-1)$$ α H = ∏ i = 0 d H i ( 2 H i - 1 ) .