First, we calculate the surface energy and the electronic structure of the [0001], [1 100] and [11 20] facets with different stoichiometries, which were previously reported as the most stable ones. Table presents the calculated surface energies corresponding to the cases shown in Figure , obtained for T=300 K and P=1 atm, using both GGA and HSE. The γ energies are in good agreement with previous reports. From Table , it can be seen how GGA and HSE provide a qualitatively different thermodynamic description of the trend in the stability of the surfaces. Within GGA the most stable orientation is the [0001] O-rich-1 while within HSE it is the [1 100] ST-3 one. The discrepancy between GGA and HSE arises from the fact that GGA is not able to correctly describe the localization of the non-compensated charge (the hole polarons) due to the intrinsic self-interaction error. The functional HSE partially corrects this error by adding a fraction of exact exchange and it is able to achieve the formation of hole polarons. These polarons have an elastic cost which rises the surface energy γ, therefore, it is expected that the O-rich facets have higher γ than the ST ones within HSE. The corresponding projected DOS shows that the polaron states lye in the band-gap of the 1st layer (solid red in Figure ).