  We give the upper and the lower estimates of heat kernels for Schr\"odinger
operators $H=-\Delta+V$, with nonnegative and locally bounded potentials $V$ in
$\mathbb{R}^d$, $d \geq 1$. We observe a factorization: the contribution of the
potential is described separately for each spatial variable, and the interplay
between the spatial variables is seen only through the Gaussian kernel -
optimal in the lower bound and nearly optimal in the upper bound. In some
regimes we observe the exponential decay in time with the rate corresponding to
the bottom of the spectrum of $H$. Our estimates identify in a fairly
informative and uniform way the dependence of the potential $V$ and the
dimension $d$. The upper estimate is more local; it applies to general
potentials, including confining and decaying ones, even if they are non-radial,
and mixtures. The lower bound is specialized to confining case, and the
contribution of the potential is described in terms of its radial upper
profile. Our results take the sharpest form for confining potentials that are
comparable to radial monotone profiles with sufficiently regular growth - in
this case they lead to two-sided qualitatively sharp estimates. In particular,
we describe the large-time behaviour of nonintrinsically ultracontractive
Schr\"odinger semigroups - this problem was open for a long time. The methods
we use combine probabilistic techniques with analytic ideas. We propose a new
effective approach which leads us to short and direct proofs.
