Beryllium (Be) is an important material with wide applications ranging from aerospace components to X-ray equipments. Yet a precise understanding of its phase diagram remains elusive. We have investigated the phase stability of Be using a recently developed hybrid free energy computation method that accounts for anharmonic effects by invoking phonon quasiparticles. We find that the hcp → bcc transition occurs near the melting curve at 0 < P < 11 GPa with a positive Clapeyron slope of 41 ± 4 K/GPa. The bcc phase exists in a narrow temperature range that shrinks with increasing pressure, explaining the difficulty in observing this phase experimentally. This work also demonstrates the validity of this theoretical framework based on phonon quasiparticle to study structural stability and phase transitions in strongly anharmonic materials.PACS numbers: 63.20. Ry, 81.30.Bx, 61.50.Ks,Elemental solids usually undergo a series of phase transitions from ambient conditions to extreme conditions [1][2][3]. Knowledge of their phase diagrams is a prerequisite for establishing their equations of state (EOS), a fundamental relation for determining thermodynamics properties and processes at high pressures and temperatures (PT). However, resolving phase boundaries is challenging for experiments given the uncertainties from several sources, especially at very high PT. Beryllium (Be) is a typical system whose phase diagram remains an open problem despite intense investigations. It assumes a hexagonal close-packed (hcp) structure at relatively low T [1]. Near the melting temperature T M (∼ 1, 550 K at 0 GPa), a competing phase with the body-centered cubic symmetry (bcc) seems to emerge [4][5][6]. However, not all experiments [7][8][9][10][11][12][13] have observed this bcc phase, causing confusion and controversies. Be is important for both fundamental research [14][15][16][17] and practical applications. Being a strong and light-weight metal, it has been widely used in a broad range of technological applications in harsh environments and extreme PT conditions, e.g., up to T > 4,000 K and P > 200 GPa [18][19][20][21][22].The bcc phase of Be was directly observed [4,5] only at T > 1, 500 K around ambient pressure before melting. Measurements of the temperature dependent resistance suggested that bcc Be is a high pressure phase and the hcp/bcc phase boundary between 0 < P < 6 GPa has a negative Clapeyron slope (−52 ± 8 K/GPa) [6]. However, recent experiments have challenged this conclusion [8,9,[11][12][13]. For example, it was reported that bcc Be was not observed for 8 < P < 205 GPa and 300 < T < 4, 000 K [13]. On the theory side, the study of Be's phase diagram using conventional methods encounters significant difficulties. The lattice dynamics of bcc Be is highly anharmonic, and the widely used quasi-harmonic approximation (QHA) and Debye model are not able to capture such effect [23][24][25][26][27][28][29]. For this reason, bcc Be and the associated hcp/bcc phase transition remain poorly understood for P < 11 GPa (density < 2.1g/...