We then extend our microscopic approach to include the effect of finite temperatures.
Experimentally it has been shown that a flat temperature dependence of the FES mode is expected in the quantum paraelectric phase (\(T<4\) K), [1]}, [2]} whereas for increasing temperatures a stiffening of the FES mode and a drop of the dielectric function are observed [3]}, [2]}, [5]}, [6]}.
We first evaluate the temperature dependent FES frequency via ab initio molecular dynamics simulations with a thermostat [7]}.
As shown in Fig. 3(a), the results with PBE provide values which are comparable with the experimental observations at high temperature, where the quantum effects are overshadowed by thermal fluctuations (see Supplementary Materials for details [8]}).
This indicates that PBE provides a realistic potential energy surface and that the effect of a slightly overestimated lattice is negligible [9]}.
To include temperature in our quantum lattice model we apply Kubo's formula for the linear response of a thermal state to a perturbation \(\hat{H}^{\prime }(t) = -Z^{*}\hat{Q_{\rm {f}}} E(t)\) , where \(Z^{*}\)  is the FES mode effective charge that we assume to be temperature independent.
The resulting polarizability takes the form: \(\alpha (\omega ,T)=-\sum _{i,j} \rho _i(T)Z^{*}|D_{ij}|^2 \times (\frac{1}{(\epsilon _j - \epsilon _i)-\omega -i\delta }+\frac{1}{(\epsilon _j - \epsilon _i)+\omega +i\delta }) \) , where the dipole matrix and the thermal density matrix are defined as \(D_{ij}=\langle \psi _i|\hat{Q}_f|\psi _j\rangle \)  and \(\rho _i(T)=e^{-(\epsilon _i-\epsilon _0)/k_BT}/\sum _j e^{-(\epsilon _j-\epsilon _0)/k_BT}\) , respectively.
The temperature dependent frequency of the FES mode, is then evaluated by averaging over the polarizability as \(\omega (T)=\frac{\int \omega Im[\alpha (\omega ,T)] d\omega }{\int Im[\alpha (\omega ,T)] d\omega }\) ; the results of this procedure with PBE are depicted in Fig. 3(a).
Even though the 2D potential only includes two degrees of freedom (\(Q_{\rm {f}}\)  and \(Q_{\rm {c}}\) ), the temperature dependence behavior is well reproduced and the typical flattening at low temperatures (\(<10\)  K) is evident.
We assigned the observed deviations of our model from the experiment to the effect of the phonon degrees of freedom that are not included in the model [3]}, [5]}, [6]}.
An attempt on improving the description of our method by including the effect of temperature on the lattice parameter bearing no significant change to the picture is discussed in Supplementary Material [8]}.
<FIGURE>