The heuristic is probably that
$ \lambda^* = \frac{\int_{-1}^1 1 \mathrm{d}x}{\int_{-1}^1 e^{-x^2}\mathrm{d}x} = \frac{2}{\sqrt{\pi} \mathrm{erf}(1)} \approx 1.339 $
Note that if $\lambda^*$ is defined as the above you have that
$ \int_{-1}^1 \lambda^* e^{-x^2} - 1 \mathrm{d}x = 0 $
For larger $\lambda$ the integral is positive, while for smaller ones the integral is negative.
So integrating your equation from $-1$ to $1$ you have that
$ \partial_t \int_{-1}^1 u\mathrm{d}x = \int_{-1}^1 \partial_x^2u \mathrm{d}x + \int_{-1}^1 f(x) u(x) \mathrm{d}x $
The first term on the right hand side is 0, using Neumann's BC and the fundamental theorem of calculus.
Writing $\bar{u}$ for the mean value of $u$ on the interval $[-1,1]$, the equation becomes
$ \frac{d}{dt} \bar{u} = 2\bar{f}\bar{u} + \int_{-1}^1 (\bar{f}-f)(\bar{u}-u) \mathrm{d}x $
The second term is roughly speaking small: initially $\bar{u} - u$ is small, and the heat equation forces an exponential decay in all higher modes: the only mode that is not expected to decay is the 0 mode, which is captured by the mean $\bar{u}$. If we just heavyhandedly set the second term to be negligible, we have that the equation is roughly governed by the ODE
$ \frac{d}{dt}\bar{u} = \bar{f} \bar{u} $
from which we get that if $\bar{f}$ is positive ($\lambda > \lambda^*$) you get exponential increase, and if $\bar{f}$ is negative ($\lambda < \lambda^*$) you get exponential decrease.