Barium zirconate (BaZrO$_3$) is one of the very few perovskites that is claimed to retain an average cubic structure down to \SI{0}{\K}, while being energetically very close to an antiferrodistortive phase obtained by condensation of a soft phonon mode at the R point of the Brillouin zone boundary. In this work, we report a combined experimental and theoretical study of the temperature dependence of this soft phonon mode. Inelastic neutron and x-ray scattering measurements on single crystals show that it softens substantially from \SI{9.4}{\meV} at room temperature to \SI{5.6}{\meV} at \SI{2}{\K}. In contrast, the acoustic mode at the same R point is nearly temperature independent. The effect of the anharmonicity on the lattice dynamics is investigated non-perturbatively using direct dynamic simulations as well as a first-principles based self-consistent phonon theory, including quantum fluctuations of the atomic motion. By adding cubic and quartic anharmonic force constants, quantitative agreement with the neutron data for the temperature dependence of the antiferrodistortive mode is obtained. The quantum fluctuations of the atomic motion are found to be important to obtain the proper temperature dependence at low temperatures. The mean squared displacements of the different atoms are determined as function of temperature and are shown to be consistent with available experimental data. Adding anharmonicity to the computed fluctuations of the Ba-O distances also improves the comparison with available EXAFS data at \SI{300}{\K}.
