We present a first-principles study, within the quasiharmonic approximation, of the thermal expansion of Be (0001) surface. The free energy is obtained from full vibrational dispersions computed by density-functional perturbation theory. Our calculation describes very well the thermal expansion in the bulk and has been checked against first-principles molecular dynamics simulations in the surface case. We do not find the large thermal expansion recently observed experimentally and we argue that the morphology of the actual surface could be less ideal than assumed. [S0031-9007(98) PACS numbers: 63.20. Ry, 68.35.Bs, 71.15.Mb, 82.65.Dp Some metal surfaces displaying anomalously large surface thermal expansion have been recently studied by combining first-principles calculations with a simplified quasiharmonic approximation [1,2]. The validity of such an approach, however, has been criticized by some authors [3] due to the very poor sampling of vibrational modes adopted in the free-energy calculation and needs to be tested against more accurate treatments. Because of its large outward relaxation [4], beryllium (0001) surface has attracted much experimental [4][5][6][7] and theoretical [8][9][10][11][12] interest. Until very recently, there was a substantial disagreement between experimental and first-principles theoretical [8-12] results on the amount of topmost interlayer expansion, theoretical calculations giving roughly half of the observed value. In a recent Letter, Pohl et al.[2] reconcile experiment and theory on this point showing that low temperature (110 K) low-energy electron diffraction (LEED) determinations of the first interlayer separation are in agreement with first-principles (zero temperature, static) calculations and that reported discrepancies at room temperature originate from a large thermal expansion of the top layer, reaching 6.7% at 700 K. These findings are very interesting and puzzling, since surface phonons show no sign of enhanced anharmonicity [2,6]. The calculation of the surface thermal expansion [2], within the simplified quasiharmonic approach of Ref.[1], agrees with experimental findings. However, the accuracy of this approach still needs to be established.In this Letter the thermal expansion of Be (0001) surface is studied from first principles, minimizing the free energy of the surface as a function of the interlayer separation in the quasiharmonic approximation (QHA). The QHA is applied by avoiding additional simplifications [1,2] thanks to the efficient calculation by density-functional perturbation theory [13] of the needed vibrational frequencies. The ability of QHA to deal with thermal effects in bulk solids has been demonstrated in recent studies [14,15]. Here its validity in the more delicate surface case, where larger anharmonicity can take place [16], is explicitly confirmed by comparing it with a first-principles molecular dynamics (MD) simulation at the highest temperature considered (700 K). The approximations involved in the determination of the surface thermal exp...