Summary: The Liouville Brownian motion (LBM), recently introduced by C. Garban et al. [Ann. Probab. 44, No. 4, 3076–3110 (2016; Zbl 1393.60015)] and in a weaker form also by N. Berestycki [Ann. Inst. Henri Poincaré, Probab. Stat. 51, No. 3, 947–964 (2015; Zbl 1325.60125)], is a diffusion process evolving in a planar random geometry induced by the Liouville measure \(M_\gamma \), formally written as \(M_\gamma (dz)=e^{\gamma X(z)-{\gamma^2} \mathbb E[X(z)^2]/2}\, dz\), \(\gamma \in (0,2)\), for a (massive) Gaussian free field \(X\). It is an \(M_\gamma \)-symmetric diffusion defined as the time change of the two-dimensional Brownian motion by the positive continuous additive functional with Revuz measure \(M_\gamma \). In this paper we provide a detailed analysis of the heat kernel \(p_t(x,y)\) of the LBM. Specifically, we prove its joint continuity, a locally uniform sub-Gaussian upper bound of the form \(p_t(x,y)\leq C_{1} t^{-1} \log (t^{-1}) \exp \bigl (-C_{2}((|x-y|^{\beta}\wedge 1)/t)^{\frac{1}{\beta -1}}\bigr )\) for \(t\in \bigl (0,\frac{1}{2}\bigr ]\) for each \(\beta >\frac{1}{2}(\gamma +2)^2\), and an on-diagonal lower bound of the form \(p_{t}(x,x)\geq C_{3}t^{-1}\bigl (\log (t^{-1})\bigr )^{-\eta}\) for \(t\in (0,t_{\eta}(x)]\), with \(t_{\eta}(x)\in \bigl (0,\frac{1}{2}\bigr ]\) heavily dependent on \(x\), for each \(\eta >18\) for \(M_{\gamma}\)-almost every \(x\). As applications, we deduce that the pointwise spectral dimension equals 2 \(M_\gamma \)-a.e. and that the global spectral dimension is also 2.