Temperature dependence of the topological phase transition of BiTeI from first principles

A topological phase transition from a trivial insulator to a $\mathbb{Z}_2$ topological insulator requires the bulk band gap to vanish. In the case of noncentrosymmetric materials, these phases are separated by a gapless Weyl semimetal phase. However, at finite temperature, the gap is affected by atomic motion, through electron-phonon interaction, and by thermal expansion of the lattice. As a consequence, the phase space of topologically nontrivial phases is affected by temperature. In this paper, the pressure and temperature dependence of the indirect band gap of BiTeI is investigated from first principles. We evaluate the contribution from both electron-phonon interaction and thermal expansion, and show that their combined effect drives the topological phase transition towards higher pressures with increasing temperature. Notably, we find that the sensitivity of both band extrema to pressure and topology for electron-phonon interaction differs significantly according to their leading orbital character. Our results indicate that the Weyl semimetal phase width is increased by temperature, having almost doubled by 100 K when compared to the static lattice results. Our findings thus provide a guideline for experimental detection of the nontrivial phases of BiTeI and illustrate how the phase space of the Weyl semimetal phase in noncentrosymmetric materials can be significantly affected by temperature.


I. INTRODUCTION
Topological phases of matter have become a thriving field of condensed matter physics, for both fundamental and applied research [1]. In three-dimensional (3D) materials, these phases are characterized by the existence of metallic surface states with peculiar properties, such as protection against backscattering, spin-momentum locking and dissipationless spin-polarized currents, as well as inverted bulk band gaps [1,2]. The discovery of experimentally tunable topological phases, whether through stoichiometric doping [3][4][5][6], hydrostatic pressure [7,8], strain [9,10], external electric fields [11,12] or interaction with light [13,14], has led to a continually growing number of proposals for promising and innovative applications relying on the refined engineering of these robust states and their associated phase transitions [15][16][17][18].
Another widely studied class of materials is the bulk Rashba semiconductors, in which a strong spin-orbit interaction combined with the absence of inversion symmetry leads to a splitting of electronic bands of opposite spin polarization [19,20]. The band extrema are shifted away from the time-reversal invariant points in the Brillouin zone, both in energy and momentum, in the plane perpendicular to the potential gradient [21]. This Rashba effect gives rise to many quantum phenomena, such as the Edelstein, spin Hall and spin galvanic effects, as well as noncentrosymmetric superconductivity [see, for example, Bahramy and Ogawa [21] and references therein]. Just * veronique.brousseau.couture@umontreal.ca like the topological surface states, the bulk Rashba split bands have their spin orientation locked perpendicular to their momentum, leading to a helical spin texture [22]. This feature makes these materials promising candidates for realizing devices involving the active manipulation of the spin degree of freedom, for both spintronics [23] and quantum computing [24].
While at first it was thought that the TPT of BiTeI occurred at a single critical pressure P C , it was later demonstrated that the lack of inversion symmetry imposed the existence of an intermediate Weyl semimetal (WSM) phase [41] between the trivial band insulator and topological insulator phases, yielding two critical pressures. Tight-binding [41] and first-principles calculations [42,43] predicted that this WSM phase could exist only within a narrow pressure range of 0.1-0.2 GPa, making its experimental detection technically challenging.
The precise value of the critical pressure of BiTeI is still elusive to this day. Experiments have located it between 2.0 and 3.5 GPa [7,[34][35][36]38], while first-principles calculations have predicted it in a slightly wider range of pressures, from 1.6 to 4.5 GPa [31,[42][43][44]. From a theoretical point of view, predicting the critical pressure comes down to finding the gap closing pressure, which is inherently dependent on the accuracy of the calculated band gap at ambient pressure, which can be biased by the well-established underestimation of the band gap by density functional theory (DFT) [45].
On the other hand, one should not forget that most first-principles calculations are done under the assumption of a static lattice, while experiments are inherently done at finite temperature. This completely overlooks the fact that electrons can interact with thermally activated phonons, resulting in a temperaturedependent shift in the electronic eigenenergies. Moreover, at T = 0 K, it neglects the contribution of the zeropoint renormalization (ZPR) to the eigenenergies, due to the zero-point motion of the ions. For narrow gap materials, this renormalization could, in principle, push the system towards a band inversion. In the case of BiTeI, this shift will directly affect both critical pressures at which the band gap will effectively close or reopen at a given temperature.
The ability of electron-phonon interaction (EPI) to induce a change in the bulk topology was first demonstrated with model Hamiltonians [46][47][48], and later analyzed by first-principles calculations for BiTl(S 1−x Se x ) 2 [49] and for the Bi 2 Se 3 family [50]. These studies were able to capture the whole complexity of the EPI, with fewer approximations. They demonstrated that, depending on the strength and sign of the different types of couplings, nontrivial topology could either be promoted or suppressed by temperature in real materials. The suppression of the topological surface state signatures at sufficiently high temperatures was observed experimentally in Pb 1−x Sn x Se [5,51,52]. The temperature-dependent band structure of BiTeI was previously investigated by focusing on the variation of the Rashba splitting rather than on tracking the TPT [53].
In this paper, we study the temperature dependence of the topological phase transition in BiTeI using firstprinciples methods, assessing for both the EPI and thermal expansion (TE) contributions. We explicitly assess the temperature dependence of the band-gap renormalization as a function of pressure, both in the trivial and the topological phases. We observe that the EPI contribution changes sign as the system undergoes the TPT, and link this behavior to the band inversion phenomenon, showing that the band extrema exhibit a distinct pressure dependence related to their leading orbital character. In turn, the TE contribution is unaffected by the leading orbital character of the band extrema, but manifests sensitivity to the topological invariant by changing sign as the system undergoes the TPT. We finally evaluate how the pressure width of the WSM phase is affected by temperature by extrapolating the total renormalization trends of each band towards the TPT, and find that the stronger renormalization of the band gap in the TI phase widens the WSM with increasing temperature.
The remainder of this paper is organized as follows. After presenting the theoretical formalism used to obtain the EPI and TE contributions to the electronic structure in Sect. II, we summarize the computational details in Sect. III. Sect. IV A presents our analysis of the TPT and WSM phase for a static lattice at T = 0 K, while Sect. IV B and IV C focus on the effects of temperature on both the band gap and the TPT. We finally discuss the implications and limitations of our results in Sect. IV D and summarize our findings.

II. METHODOLOGY
The temperature dependence of the band-gap energy, E g , for a given constant pressure, can be approximated at first order by the sum of two distinct contributions [54], namely the thermal expansion, ∆E TE g , and the renormalization of the electronic eigenenergies by electron-phonon interactions at constant volume, ∆E EPI g :

A. Electron-phonon interaction
In many-body perturbation theory [55], the temperature dependence of the electronic eigenenergies, ε kn , induced by electron-phonon interaction can be captured, under certain assumptions [56], by the real part of the electron-phonon self-energy, Σ EPI , evaluated at the quasiparticle energy: where ε 0 kn is the unperturbed eigenvalue of a Kohn-Sham (KS) eigenstate with wavevector k and band index n, computed at the fixed, relaxed geometry. The Hartree atomic unit system is used throughout this paper ( = m e = e = 1).
In the nonadiabatic Allen-Heine-Cardona (AHC) framework [54,57,58], the self-energy at the second-order of perturbation is the sum of the Fan and Debye-Waller (DW) contributions [ Fig. 1]: where the former captures a frequency-dependent interaction with two first-order EPI vertices, and the latter captures a static one with a single second-order vertex. Evaluating the previous equation at absolute zero temperature yields the ZPR, which arises from the zero-point motion of the ions. For a detailed review of the AHC methodology, we refer our readers to the works of Poncé et al. [59,60] and to the extensive review by Giustino [61].
Assuming that the quasiparticle energy is close to the unperturbed electronic energy, one can apply the onthe-mass-shell approximation [62], in which the Fan selfenergy is evaluated at the poles of the Green's function, namely at the frequency ω corresponding to the bare eigenvalue ε 0 kn , yielding Furthermore approximating the fully interacting electronic wavefunction by the noninteracting KS-DFT wavefunction, one obtains the standard result for the retarded Fan self-energy [59,61]: in which the contributions from all phonon modes with frequency ω qν are summed for all wavevectors q and branches ν in the Brillouin zone (BZ). The whole temperature dependence of this expression is captured by the bosonic and fermionic occupation factors, respectively n qν and f k+qn . The parameter η k = η sgn(ε 0 kn − µ), where η is real and positive, preserves causality by correctly shifting the poles of the Green's function in the complex plane. Equation (5) implies the limit of an infinite number of phonon wavevectors (N q → ∞), leading to a vanishingly small value of η. V (1) qν is the first-order self-consistent change of the local KS potential induced by the collective atomic displacement R ν κα (q) along a given phonon mode [56,59,61]: with In these expressions, R lκα denotes the displacement of atom κ, located in unit cell l with lattice vector R l , in cartesian direction α; R ν κα (q) therefore describes the collective atomic displacement of along the qν-phonon mode. V (1) κα (qν) refers to the contribution of a displacement of atom κ in direction α to the full first-order potential V (1) qν . U ν,κα (q) is the qν-phonon displacement vector, which is related to the phonon eigenvector, ξ ν κα (q), through: with M κ being the mass of atom κ. The static Debye-Waller self-energy is defined as the second derivative of the potential with respect to two atomic displacements, evaluated at the first order in perturbation theory [56] : where the second-order perturbation potential is with the derivatives defined as in Eq. (7). Only the phonon occupation factor contributes to the temperature dependence of this term since it does not involve any intermediate electronic state [as can be seen from Fig. 1(b)]. The constant term inside the brackets accounts for the ZPR contribution. The numerical evaluation of this second-order derivative is a computational bottleneck of density-functional perturbation theory (DFPT) [60]. It can, however, be circumvented by applying the rigid-ion approximation (RIA), which supposes that the potentials created by each nucleus are independent of each other. Within this approximation, thanks to the translational invariance of the lattice, V (2) qν can be expressed in terms of the firstorder derivatives entering the Fan self-energy [59] [see Eq. (6)]. The consequences of this approximation on the ZPR of diamond have been discussed by Poncé et al. [59]. In this framework, the Debye-Waller self-energy becomes [49,59]: with and We finally evaluate the sum on band index n in Eq. (5) and (11) using the semi-static approach described in Poncé et al. [59] and Antonius et al. [56]. We compute the full, nonadiabatic contribution for all bands up to a band index M , which we choose to be well-separated in energy (i.e., more than 10 eV) from the first conduction band. For the high-energy bands with band index n > M , the phonon frequencies in the denominators of Eq. (5) and (11) are negligible with respect to the electronic eigenenergy difference. These bands can, therefore, be treated within the adiabatic approximation, neglecting the phonon frequencies. Furthermore, this sum over an arbitrarily large number of empty states can be replaced by the solution of a linear Sternheimer equation for the subspace orthonormal to the low-energy states (n ≤ M ), thus significantly reducing the numerical cost of the calculation [63]. We estimate the relative error on the renormalization induced by this semi-static treatment of the high energy bands to be 1-2% at most.

B. Thermal expansion
The thermal expansion contribution was evaluated through the quasi-harmonic approximation (QHA) [64,65]. Within this framework, the only departure from harmonicity occurs through the explicit volume dependence acquired by the phonon frequencies. It should be understood that the purpose of the QHA is to deliver the leading order expression of the thermal expansion coefficient α, rather than capturing the full anharmonic effects present in the crystal [66], as would do nonperturbative methods [56].
The thermal expansion coefficient, α, is a rank 2 tensor that relates a small temperature increment ∆T to the strain, , it induces on the lattice: where i, j = 1, 2, 3 are cartesian directions. Hence, the most general definition of α ij is where the derivative is evaluated at constant stress σ, which is usually taken to be a constant pressure P . The diagonal components of the strain tensor describe the relative change in the length of the lattice parameters: ii = ∆a i /a i . The volumic expansion coefficient β can therefore be obtained by taking the trace of α: In the most simple case of cubic symmetry, all α ii are equal and can be expressed by [64] where κ 0 is the bulk compressibility at equilibrium volume, V 0 is the primitive cell equilibrium volume, N q is the number of phonons in the homogenous grid used to sample the BZ and c qν (T ) is the specific heat of the qνphonon mode at temperature T . We have also introduced the volumic mode Grüneisen parameters, defined as [64] γ The bulk Grüneisen parameter, often referred to in the literature, is defined as Thus, for cubic systems, the linear thermal expansion coefficient, α(T ), Eq. (17), is simply proportional to two scalar parameters, γ V and κ 0 , the former governing the sign of the thermal expansion.
For materials belonging to noncubic symmetry groups, one must consider the most general case, in which the mode Grüneisen parameters take a tensor form and capture the variation of the phonon frequencies with respect to a given strain: One can also define bulk Grüneisen parameters γ ij following the same procedure as in Eq. (19). For axial crystals, which include the case of BiTeI's trigonal symmetry, the thermal expansion coefficient tensor has two distinct nonvanishing components: α 11 = α 22 and α 33 . Because of this anisotropy, the resulting linear thermal expansion coefficients along crystallographic directionsâ andĉ capture a subtle interplay between the vibrational and elastic properties of the material [64,67,68]: where C σ is the heat capacity at constant stress and s ij is the ij-coefficient of the elastic compliance tensor. In these expressions, the subscripts a and c refer to the length of the unit cell along (â 1 ,â 2 ) andĉ. γ a1 formally refers to a strain applied in only one direction perpendicular in the crystallographic axis.
From Eq. (20), the directional mode Grüneisen parameters γ a1 qν and γ c qν take the form In the previous expressions, the phonon frequency derivatives are evaluated with respect to the variation of only one lattice parameter, the other one remaining fixed at its static equilibrium value. One should finally note that, in contrast to Eq. (20), the derivative entering γ a1 qν is taken by changing the a cell dimension, which affects both lattice vectors in the basal plane, hence the factor 1/2. We finally note that the derivation of the Grüneisen formalism [64] neglects the zero-point energy of the phonons, such that there is no thermal expansion at T = 0 K.

A. First-principles calculations
All first-principles calculations were performed using the ABINIT software package [69]. The bulk electronic structure was obtained within DFT, while response function and electron-phonon coupling calculations were performed within DFPT [70,71]. We used a maximum plane-wave energy of 35 hartree and sampled the BZ using a 8 × 8×8 Monkhorst-Pack k-point grid. Spin-orbit coupling was taken into account as it is necessary to obtain the Rashba effect and since it has been shown to strongly affect both electronic [29] and vibrational [72] properties of BiTeI. We used Hartwigsen-Goedecker-Hutter (HGH) fully relativistic norm-conserving pseudopotentials [73], including explicitly the semi-core 5d electrons for Bi. The electron-phonon self-energy was computed with the ElectronPhononCoupling module.
Throughout this paper, we rely on the generalized gradient approximation of the Perdew-Burke-Ernzerhof (PBE-GGA) functional [74], although it has been shown to overestimate the lattice parameters at ambient pressure for this particular material [42,53,75]. Since the purpose of this paper is to investigate the temperature dependence of the TPT, we chose the functional that gives us the best overall agreement with experiment for the bare band gap and the projected orbital character of the band extrema at ambient pressure, as well as for the predicted critical pressure P C1 . For detailed results and a complete discussion about our choice of PBE-GGA functional, see Appendix A.

B. Structural properties
The unit-cell geometry has been optimized until the resulting forces on all atoms were lower than 10 −5 hartree/bohr 3 . Due to the layered nature of BiTeI, the in-plane and normal lattice parameters of BiTeI do not vary isotropically under hydrostatic pressure [7]: c decreases more rapidly than a within the first GPa applied. In order to allow for this nonhomogenous variation of the lattice parameters, the external pressure was modelized by applying an isotropic stress tensor. The lattice was fully optimized for nine different pressures between 0 and 5 GPa. The resulting cell volumes were then fitted with a Murnaghan equation of state (EOS) [76] in order to validate the theoretically applied pressure. For all pressures below 3.5 GPa, the discrepancy between the EOS pressure and the applied pressure was less than 0.1 GPa, and less than 0.06 GPa for pressures surrounding the TPT. The lattice structure at intermediate pressures was obtained by interpolating from the results of those nine calculations. The optimized lattice structure at all pressures considered in this paper can be found in Table II of Appendix A.

C. Lattice dynamics and electron-phonon coupling
Phonon frequencies and electron-phonon matrix elements were computed with a 12 × 12 × 12 homogeneous, Γ-centered Monkhorst-Pack q-point grid. Explicit calculations were done at seven different pressures, spanning both the trivial (0.0, 0.5, 1.0, and 1.5 GPa) and topological (3.0, 3.5 and 5.0 GPa) phases. The pressure dependence of the phonon frequencies for the Raman-active modes is in good agreement with experimental data [see Tables IV and V of Appendix A].
In a typical nonadiabatic AHC calculation, one aims for a value of η smaller than the highest phonon energy (typically around η = 0.01 eV, see, for example, Nery et al. [77]). This however requires a very dense q-point sampling. As discussed in the Appendix of Antonius et al. [56], when using a coarser q-point grid, one risks artificially emphasizing the contribution of a discrete number of terms in Eq. (5) where, for a given qν-phonon, the value of ε 0 kn − ε 0 k+qn ± ω qν gets vanishingly small, hence inducing rapid fluctuations of the self-energy [56]. These peaks are a numerical artefact of the q-point sampling, and the presence of such fluctuations in the vicinity of the bare electronic eigenvalue ε 0 kn could lead to an overestimation of the renormalization. In contrast, increasing η without caution can lead to an unphysical flattening of the self-energy, suppressing physically relevant features. We therefore chose η=0.1 eV, which allows the self-energy to be a smooth function of frequency, without suppressing physically relevant features. We verified that this choice was suitable for the claims reported in this paper.
Lastly, our first-principles methodology relies on the assumption that the electron-phonon interaction can modify the electronic eigenenergies, but cannot by itself promote an electron to an excited state. Such a possibility would modify the resulting electronic density and would require one to evaluate the nondiagonal elements of the self-energy, Σ kn,n [see Eq. (4) of Antonius et al. [56]], allowing n = n. This procedure was, however, out of reach for the present work. Thus, we restricted our investigation of electron-phonon coupling to pressures where the predicted bare band gap was greater than the highest phonon energy (E g > 20 meV).

D. Thermal expansion
The directional Grüneisen parameters γ a1 and γ c were evaluated using an 8 × 8 × 8 Γ-centered Monkhorst-Pack q-point grid. Note that the Grüneisen parameters require a smaller q-point grid than the EPI interaction since they only depend on the phonon frequencies, which converge faster with the q-point sampling.
The derivatives in Eq. (22) were computed by central finite difference using three different volumes per direction, per pressure. For each calculation, the internal atomic coordinates in theĉ direction were optimized at fixed volume, using the same convergence criterion as for the equilibrium structure. The compliance constants s ij were obtained by computing the strain-strain derivatives of the total energy within DFPT [78]. We used the so-called relaxed-ion compliance tensor, which takes into account the relaxation of the atomic coordinates within the unit cell. This typically lowers the components of the elastic stiffness tensor (usually referred to as the elastic constants), as the resulting stress on the atoms is reduced by the relaxation process. Since the compliance tensor is the inverse of the elastic stiffness tensor, it typically increases the elastic compliance when compared to the clamped-ion result [79].

A. Topological phase transition in the static lattice
BiTeI is a layered material composed of alternating Bi, Te and I planes stacked along the high-symmetry crystallographicĉ axis. This small band-gap semiconductor belongs to the noncentrosymmetric trigonal space group P3m1 (no. 156, mp-22965). The crystal structure and first BZ are shown in Fig. 2.
A common signature of a TPT is an inversion of the orbital character of the band structure in the vicinity of the band gap [1]. In the case of BiTeI, a schematic analysis of the band splitting showed that the inverted bands should be of p z character [31]. Fig. 3 displays the leading orbital character of the last valence and first conduction bands, for different pressures throughout the studied range. We computed the p z character by projecting the wavefunction on the l = 1, m = 0 spherical harmonic centered around the different atoms in the unit cell. To emphasize the band inversion, Fig. 3 shows the relative projected character on each band, namely the normalized difference between the Bi-6p z and Te/I-5p z projections. In the H-A and A-L directions, which both exhibit Rashba splitting, the valence band maximum (VBM) is dominated by p z states of Te and I (dark brown), whereas the conduction band minimum (CBM) has Bi-p z (light yellow) character. It is a trivial band insulator, with topological index Z 2 =(0;000). The Z 2 topological index was computed by tracking the evolution of the hybrid Wannier charge centers, using the Z2Pack software package [80,81]. For more information about the definition of the Z 2 topological index for 3D topological insulators, we refer to Sections II-C and IV-A of the review by Hasan and Kane [1].
As pressure is increased, the band gap progressively decreases until it closes at P C1 =2.08 GPa [ Fig. 3 Fortunately, from previous studies [41,42], it is known that the gap closing occurs along symmetry line H-A (with k z = π/c), and that it will reopen along a symmetry line rotated by π/6, without symmetry constraints in the k z direction (namely, in the A-L-M mirror plane). In the following, we define the P C2 plane as the plane parallel to the H-A-L plane, with k z = 1.07698 π/c, which is the k z coordinate of the second band touching point displayed in Fig. 3(c), where the Weyl nodes' annihilation occurred. This plane corresponds to the shaded plane in Fig. 2(b). In a similar fashion, we identify the H-A-L plane as the P C1 plane.
In the vicinity of the TPT, the valence and conduction band extrema show a mixed character, with an almost equal contribution of Bi-p z and Te/I-p z states [ Fig.3  Gap energy (meV) Variation of the indirect band gap with pressure. In the trivial insulator phase (Z2 = 0), the minimal band gap is located in the H-A-L plane (diamonds), where the first band touching point will occur at PC1=2.08 GPa. The system enters the topological phase (Z2 = 1) once the pairs of Weyl points have annihilated at the second band touching point, at PC2=2.28 GPa. For higher pressures, the minimal band gap is now located in the A-L-M mirror plane, with kz = 1.07698 π/c (circles). The gray shaded area represents the WSM phase width in the static lattice approximation. See Section IV A for more details.
c)], foreshadowing the band inversion [ Fig. 3(d)] which occurs for P >2.28 GPa. The system is then a strong topological insulator, with Z 2 =(1;001). In the TI phase, the minimal band gap is shifted from the H-A-L plane (k z = π/c) to the P C2 plane (k z = 1.07698 π/c), where the Weyl nodes' annihilation occurred. The full pressure dependence of the band gap for the static lattice is shown in Fig. 4. For the static lattice at T = 0 K, we thus obtain a WSM phase width of 0.2 GPa, in good agreement with previous calculations [41][42][43]. Our value for the first gap closing P C1 =2.08 GPa agrees with many experimental estimations [7,34,38] but is slightly lower than those of other experimental [35,36] and computational works [42,43,82]. For comparison with experiments, one should note that, as the TPT of BiTeI occurs under hydrostatic pressure, direct observation of the topological surface states with angle-resolved photoemission spectroscopy (ARPES) would be challenging from an experimental point of view. As a consequence, experiments rather focused on indirect signatures of these states, for example, a broad minimum of resistivity [34,35], which do not allow a precise estimation of P C . The difference with other computational works can be explained by the use of experimental pressure-dependent lattice parameters instead of the DFT-relaxed ones [43], as well as by the choice of projector-augmented wave (PAW) pseudopotentials without semi-core states for Bi [42,82] and the use of the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional [82] to compute the static gap at ambient pressure. Indeed, the calculated critical pressure is quite sensitive to the choice of functional, as it may signifi-cantly modify the starting gap at ambient pressure. One must also recall that the DFT gap has fundamentally no physical meaning; a proper, theoretically adequate calculation of the ambient pressure gap would require using a many-body methodology like GW. Nevertheless, for the scope of this paper, we are more interested in the relative variation of the critical pressures rather than in their absolute values. These slight discrepancies with other works do not, therefore, undermine our conclusions.
Our value of P C1 corresponds to a relative volume (V /V 0 ) compression of ∼ 12%, in good agreement with other calculations relying on the PBE functional [31]. We also verified that the relative volume compression within the WSM phase was consistent with the original prediction of Liu and Vanderbilt [41].
B. Temperature dependence of the topological phase transition Within this paper, we track the temperature dependence of the critical pressures P C1 and P C2 by evaluating how the electron-phonon interaction and the thermal expansion will affect the indirect band-gap value. Our analysis is based on several assumptions. First, we suppose that the Z 2 topological index computed at T = 0 K for the static lattice remains unchanged unless the renormalization closes the gap. Second, since the system lacks inversion symmetry regardless of temperature, we assume that Liu and Vanderbilt's argument [41] is still valid, such that a gapless WSM phase still separates the two insulating phases. We finally suppose that the k z value of the P C2 plane remains unchanged by temperature, as it is unlikely that the Weyl nodes' annihilation point in the BZ changes significantly from the static result. We verified that the position of both band extrema along the A-L line in the TI phase does not change significantly over the range of temperatures and pressures studied, between the P C1 and P C2 planes. A more thorough analysis of the temperature dependence of the band-gap location within the P C1 and P C2 planes, as well as the pressuretemperature dependence of the Rashba splitting, can be found in Appendix B.
A more formal approach would also require one to redefine the topological invariant for finite temperatures, relying on the density matrix rather than on the static ground state wavefunction, since temperature modifies the occupation of the electronic eigenstates [83][84][85]. The question of the persistence of quantum topological order at finite temperatures, compared to a nonquantized "classical" topological order, is still under investigation [86].
When considering the effect of temperature on the TPT within our framework, there are thus two possible outcomes, which depend on the sign of the band-gap renormalization. In the first case, a negative correction implies that the band gap closes with increasing temperature, which can both drive a topologically trivial sys-tem towards a band inversion and stabilize an already inverted band structure, thus globally promoting the topological phase. On the other hand, a positive correction will favor the trivial phase: the band gap opens with increasing temperature, preventing the band inversion in a trivial system and destabilizing a TI band structure until the band inversion is reversed. A more detailed version of this argument can be found in our proceedings paper [87].
In the case of a pressure-induced TPT with an intermediate WSM phase, both critical pressures, P C1 and P C2 , will be modified by temperature: a negative renormalization favoring the topological phase will diminish the amount of pressure required to close the band gap, while a positive renormalization promoting the trivial phase will increase it. The WSM phase width will also be affected by the relative strength of the renormalization between the trivial and topological phases. The combined effect of the strength and sign of the renormalization will determine if the WSM phase is widened or narrowed by increasing temperature.

Electron-phonon interaction
In a typical semiconductor, with a nearly parabolic dispersion near well-defined band extrema and a sufficiently wide band gap, the Fan contribution usually dominates the self-energy. One then expects from Eq. (5) to observe a decrease of the band-gap energy with increasing temperature: the VBM is repelled by neighboring occupied states with similar but lower energy, leading to a positive renormalization. Similarly, the CBM is repelled by neighboring unoccupied states with higher energy, yielding a negative renormalization. Couplings between states with different occupation factors are disadvantaged by the large energy difference in the denominator, due to the presence of the gap. This behavior is sometimes referred to in the literature as the Varshni effect [88].
In BiTeI, the temperature dependence of the gap cannot be inferred by such a heuristic argument, since, on the one hand, the small band gap emphasizes the weight of couplings between the subsets of occupied and unoccupied bands compared to a typical semiconductor, and, on the other, the Rashba splitting creates regions in the BZ where a phonon with finite, nonvanishing wavevector q can couple electronic states with very similar eigenenergies.
The EPI contribution to the VBM (center panel), CBM (top panel) and total band gap (bottom panel) is shown in Fig. 5. In order to track any temperatureinduced change to the gap location and to accurately capture the renormalization for the minimal gap, the EPI was computed in the H-A-L plane for the trivial phase, and in the P C2 plane for the TI phase, for electronic states along the solid red and dashed blue paths shown in Fig. 2(b).
In the trivial insulator phase, both the VBM and the CBM are shifted towards higher energies. While the CBM renormalization is almost pressure independent for a given temperature, the VBM is strongly affected by pressure, its renormalization dropping by roughly a factor of 3 between 0 and 1 GPa. This behavior is most likely linked to the higher compressibility along the c axis at lower pressures [7]. The strength of this variation could, however, be amplified by the overestimated compressibility predicted by the PBE functional for BiTeI [75]. In our case, this results in a total gap renormalization that changes sign within the trivial phase, going from negative to positive shortly after 0.5 GPa. Nevertheless, the EPI contribution to the gap remains quite small compared to the bare gap energy at ambient pressure (at most 20 meV at all temperatures). As argued before, such an opening of the gap with temperature when approaching the TPT will delay the first gap closing, thus moving P C1 towards a slightly higher value.
In the TI phase, the renormalization displays a different behavior: the VBM and CBM renormalizations have opposite signs, both contributing to a decrease of the gap energy with increasing temperature. Again, this behavior is not favorable for the nontrivial case as it works to revert the already inverted gap, driving the system back to the trivial phase. One should also note that, in this phase, it is the VBM that exhibits a pressureindependent behavior, while, for the CBM, the negative renormalization increases steadily with increasing pressure. This seemingly distinct behavior is simply a signature of the band inversion: it can be understood by recalling that the orbital decomposition of the state being corrected at the VBM in the TI phase has the same character as the CBM in the trivial phase, as was emphasized in Fig. 3. Similarly, the CBM in the trivial phase has the same leading orbital character as the VBM in the TI phase. From this point of view, the band extremum dominated by Bi-p z character exhibits a quasi-pressureindependent renormalization, while the extremum with Te/I-p z character is gradually shifted towards lower energies throughout the studied pressure range. It is also rather peculiar that only the Te/I states are affected by the change of topology, which can be seen in the temperature-dependent renormalization at a given pressure going from positive to negative as the system goes through the phase transition. A more thorough explanation of this sign change would require one to probe the EPI in the gap-closing region, which was not possible in the scope of this paper. Nevertheless, if considering the gap behavior in both phases, we can infer that the EPI globally promotes the trivial phase, pushing the TPT towards higher pressures as temperature increases.

Thermal expansion
The TE contribution to the total renormalization is displayed in Fig. 6, for the CBM (top panel), VBM (center panel), and indirect band gap (bottom panel). When solely looking at the band renormalization, one finds a seemingly monotonic behavior, both bands being shifted towards lower absolute energies. The pressure-dependent renormalization at a given temperature does not seem affected by the leading orbital character of the band extrema, as was the case for EPI. In the TI phase, the renormalization rate at a given pressure is greater for the CBM than for the VBM, hence a global decrease of the gap energy with increasing temperature. In the trivial phase, the indirect band-gap renormalization shows a far more erratic behavior: while for most pressures the renormalization rate is greater for the VBM than for the CBM, hence a gap opening with increasing temperature, at 0.5 GPa we obtain the opposite trend. This is reminiscent of the sign change observed in this pressure range for the EPI contribution [see Fig. 5, lower panel], which we again attribute to an overestimated compressibility in the low-pressure regime. For a layered structure like BiTeI, a higher compressibility will mainly enhance the s 33 compliance constant entering the definition of α c (Eq. (21)), which will, in turn, modify the thermal ex-  Gap ren (meV) TE contribution to the temperaturedependent renormalization for the conduction-band minimum (CBM, top panel), valence-band maximum (VBM, center panel) and total band gap (lower panel) of BiTeI, as a function of pressure. A zero renormalization (black dashed lines) refers to the static lattice results. Black arrows indicate increasing temperature. For visual reference, the shaded gray area locates the WSM phase boundaries from the static calculation. Note that the Grüneisen formalism neglects the zero-point phonon vibrational energy, hence the identically zero correction at 0 K (gray). pansion along this axis and the resulting gap renormalization. As discussed more thoroughly in Appendix A 1, the lattice structure is reasonably well described by PBE for P ≥ 1 GPa. Thus, in the following section, we only consider the normalization trends for pressures of 1 GPa and higher to evaluate the temperature-dependent critical pressures. One should finally note that the TE contribution has the same order of magnitude as the EPI contribution, such that it cannot be neglected when investigating the temperature dependence of the gap, as it is often done for materials with light ions.
We recall here that, as mentioned in Section II B, the Grüneisen formalism neglects the zero-point phonon vibrational energy, such that the T = 0 K lattice parameters are identical to the static ones. We estimated the missing zero-point thermal expansion contribution to the band-gap renormalization using the high-temperature extrapolation method displayed in Fig. 2

of Cardona and
Thewalt [88]. In this method, the renormalized lattice parameters at T = 0 K are evaluated by extrapolating the high-temperature slope of ∆a(T )/a and ∆c(T )/c to 0 K. This led to a nearly constant contribution of ∼ 2 meV in the trivial phase and ∼ −5 meV in the topological phase. This missing contribution does not affect the qualitative trends previously described.
The renormalization trends observed when approaching the TPT in both insulating phases, a gap opening in the trivial phase and a gap closing in the topological phase, are precisely the trends observed for the EPI contribution. Consequently, just like EPI, TE delays the TPT, moving both critical pressures towards higher values with increasing temperature. Excluding P < 1 GPa from the analysis for previously stated reasons, these trends agree qualitatively with Monserrat and Vanderbilt [53], despite being smaller in absolute value. This could partly be attributed to the fact that Monserrat and Vanderbilt minimized the Gibbs free energy independently for each lattice parameter, while our methodology [see Eq. (22], correlates the expansion along a given axis with the "thermal pressure" [67] induced on the lattice by varying both lattice parameters, through γ a1 and γ c . This discrepancy could also be explained by anharmonic contributions to the phonon free energy, which our perturbative methodology does not capture. A more refined and numerically accurate calculation of the anisotropic thermal expansion coefficients would rely on the full temperature-dependent phonon-perturbed potential, taking into account all anharmonic interactions, as implemented in the Temperature Dependent Effective Potential method (TDEP) method [89].

C. Weyl semimetal phase evolution
The temperature dependence of the WSM phase width was assessed by combining the results displayed in the two previous sections to evaluate the full gap energy renormalization at each pressure and temperature. For each k point on the paths considered [see Fig. 2(b)], the EPI and TE corrections were added to the static band structure before computing the resulting temperaturedependent band gap.
Since, for technical reasons discussed in Section IV A, we could not directly probe the EPI in the gap closing region, we tracked the temperature dependence of P C1 and P C2 by extrapolating the total renormalization behavior for each band extrema from neighboring pressures towards the TPT in the trivial and topological phases. This correction was then added to the bare gap energy computed with the static lattice (black markers with denser pressure sampling), as shown by the colored dashed lines in Fig 7. The temperature-dependent critical pressures extracted from Fig. 7 are explicitly shown in the pressuretemperature phase diagram of Fig. 8 (black markers). For visual reference, the critical pressures obtained with the static lattice are shown with dotted lines. The nu- Gap energy (meV) merical values for P C1 (T ) and P C2 (T ) are detailed in Table VI of Appendix C, as well as the WSM phase boundaries obtained by considering only the contribution from EPI or TE. As discussed in the previous sections, both EPI and TE are unfavorable to the topological phase. We therefore observe without surprise that both P C1 and P C2 are shifted towards higher pressures. However, as the renormalization from both contributions is stronger in magnitude in the topological phase, the temperature-dependent variation is much weaker for P C1 than for P C2 . Thus, the WSM phase width increases with temperature. This effect is already sizable at T = 0 K, where the zero-point motion of the ions induces a 15% increase in the WSM phase width, compared to the static lattice. It has almost doubled by T = 100 K, despite the absolute WSM phase width remaining small at 0.39 GPa.
The only deviation to this behavior is the slight decrease of P C1 above room temperature. This can be explained by the fact that, for P = 1 and 1.5 GPa, EPI and TE do not shift the band extrema in the same way in the BZ. Therefore, in this pressure range, the total renormalization does not sum up to the individual corrections.
One should also take into consideration the fact that, at finite temperature, thermal excitations do not allow one to physically distinguish between an insulating and a gapless phase, if the band-gap energy is smaller than k B T . Thus, for each temperature, we defined a crossover region between topologically distinct phases by evaluating the pressure difference between a given static lattice critical pressure and the pressure at which E g = k B T in the appropriate insulating phase, using the data from Fig. 4. In Fig. 8, these regions identify the phase space where, within our methodology, the topological index cannot be properly assessed. Lastly, we also added to Fig. 8 a color intensity gradient in the high-temperature regime of the WSM phase, to emphasize the fact that our calculations did not explicitly treat the temperature dependence within this pressure region. Therefore, there remains some uncertainty about the topological nature of BiTeI in this region of the pressure-temperature phase space.
From these results, we can conclude that temperature effects do not restrain the experimental detection of the WSM phase. In the low-temperature regime, where quantum effects can more easily be observed thanks to the reduction of thermal noise, the WSM phase widening remains however small, despite being sizable when compared to the static lattice phase width. The experimental signatures of the TI phase will be found at higher pressures with increasing temperature. These qualitative temperature dependencies can assist the experimental effort when designing experiments the purpose of which is to detect signatures of a topologically nontrivial phase. Moreover, they shed some light on the subtle interplay between EPI and TE and their effect on the phase space of topological phases in noncentrosymmetric materials. Since the presence of gapless modes characterizes both the WSM and Z 2 TI phases, additional experiments focusing on observables that can physically distinguish between these two nontrivial phases, and that are feasible under isotropic hydrostatic pressure, could provide more insight about this phenomenon.

D. Discussion
The purpose of this paper was to determine the qualitative trend of the temperature-dependent renormalization of P C1 and P C2 : the global behavior of the TPT, more than the precise numerical values of the critical pressures, was our main target. From this point of view, there are several limitations to our analysis, the first of them being the q-point sampling for the EPI calculation. BiTeI is a polar material, i.e., has nonvanishing Born effective charges, with infrared-active phonon modes. At long wavelength, |q| → 0, atomic displacements along longitudinal optical (LO) modes generate a macroscopic electric field throughout the material, which can couple to the electrons through dielectric interaction. This particular type of EPI, known as the Fröhlich interaction [90], has been shown to dominate the ZPR for polar materials [91,92] and to cause an unphysical divergence of the ZPR within the adiabatic AHC framework [59], in which Eq. (5) is approximated by neglecting the phonon frequencies in the denominators.
Thus, for polar materials, evaluating the precise contribution of the long-range Fröhlich interaction to the ZPR within the nonadiabatic AHC framework requires a very dense q-point grid, which for technical reasons was not available to us at the time of this study. The methodology developed by Nery and Allen [91] to estimate the polar contribution missing from a coarse q-point grid calculation could also not be straightforwardly applied to BiTeI for two reasons. First, the band extrema are not located at the zone center and spin-orbit coupling is mandatory to describe Rashba splitting correctly. Second, the validity of the parabolic effective mass approximation near the band extrema is also questionable for this material, especially at higher pressures. Our results should consequently be understood as a lower bound to the full EPI contribution.
The extrapolation technique used to evaluate the temperature dependence of the critical pressures is also by itself limited by the assumption that the renormalization behavior we observe in both insulating phases varies monotonically towards the TPT. This assumption could only be validated by an explicit calculation of the EPI in the close vicinity of the TPT, which would require one to lift the diagonal self-energy approximation, as excitations across the gap will no longer be negligible. The addition of these nondiagonal coupling terms to the self-energy, which we conjecture will have a significant contribution, could alter the general trends observed in Fig. 7. These extra terms will also be crucial to fully characterize the effect of EPI on the Te/I-p z band extrema near the critical pressures (VBM for Z 2 = 0, CBM for Z 2 = 1; see Fig. 5 and Section IV B 1).
Lastly, our analysis supposed a perfect semiconductor crystal, while experimentally BiTeI has been shown to be naturally n-doped due to small deviations in stoichiometry [33]. The presence of such doping charge could slightly affect the phonon-perturbed crystal potential, which could, in return, modify the matrix elements entering the EPI self-energy. We do not however expect this self-doping to have a strong impact on the qualitative trends observed throughout this paper, as long as the amount of defects remains small compared to substitutional doping.

V. CONCLUSION
In the present paper, we have characterized the temperature dependence of the topological phase transition in BiTeI using first-principles methodologies based on DFPT. The electron-phonon interaction was obtained with the nonadiabatic AHC methodology, while the thermal expansion was computed using the Grüneisen parameters formalism.
The indirect band-gap renormalization induced by EPI changes sign as the system undergoes the phase transition, opening the gap in the trivial phase and closing it in the TI phase. The band extremum with leading Bi-6p z character displays a quasi pressure-independent behavior, while the Te/I-5p z extremum manifests a strong sensitivity to both pressure and topology.
The thermal expansion contribution to the band-gap renormalization has the same order of magnitude as the EPI contribution. At a given temperature, both band extrema are affected in a similar manner throughout the pressure range, regardless of their leading orbital character. The resulting band-gap renormalization nevertheless captures the change of topology: it mainly increases the indirect gap in the trivial phase and reduces it throughout the topological phase.
The combined effect of EPI and TE globally moves the TPT of BiTeI towards higher pressures. The resulting temperature-dependent renormalization is stronger in magnitude in the TI phase compared to the trivial phase, such that the intermediate WSM phase is widened by temperature. Clear signatures of the Z 2 topological insulator phase should, therefore, appear at higher pressures as temperature is increased.
Overall, our results indicate that temperature effects are not negligible for materials with heavy ions and must be accounted for when engineering devices relying on topological properties. Our findings can also aid the search for experimental evidence of the topologically nontrivial phases of BiTeI, and more generally reveal how temperature can have a substantial influence on the phase space of topological phases of noncentrosymmetric materials.
Note added. Recently, we learned of an article by Lihm and Park [106] which is relevant to this research field. Since the first prediction of the TPT in BiTeI [31], numerous studies of ground-state properties of bismuth tellurohalides (BiTeX, X=Cl, Br, I) have appeared in the first-principles literature [42,72,[93][94][95][96]. Many of them have discussed the fact that PBE-GGA does not capture the Van der Walls bonding properly along the normal axis for this class of materials, resulting in a general overestimation of the lattice parameters [42,75], up to 7% for thê c axis [42]. PBE-GGA was also criticized as yielding an unrealistic compressibility for BiTeI at low pressures [75]. It has also been shown that dispersion-corrected DFT, which accounts for Van der Waals forces through a semiempirical dispersion potential, reproduces more accurately the pressure dependence of the lattice structure for bismuth tellurohalides [44]. These studies did not, however, investigate the impact of these corrections on the phonon frequencies, nor on the electron-phonon interaction. In the following, we include the Van der Waals interaction within the DFT framework, using dispersioncorrected DFT-D3BJ, where a semi-empirical dispersion potential is added [97] in conjunction with Becke-Jonhson damping [98,99] to avoid short-range diverging behavior. DFT-D3BJ was shown to be one of the most reliable dispersion-corrected methods to describe vibrational properties [100]. Both the trivial (0, 1.5 GPa) and the topological (3.5, 5 GPa) phases are considered.

Structural properties
The lattice parameters, bare band gap energies and equilibrium distance between the Te/I and Bi planes for 0 and 5 GPa are shown in Tables I and III, with comparison to available experimental data. In good agreement with previous works [42,53,75], PBE-GGA overestimates the lattice parameters at ambient pressure more than what is typically expected from DFT calculations, at 2.3% for a and 6.7% for c, while the dispersioncorrected method reproduces more accurately the experimental structure. At 5 GPa, both methods agree within 2% of experimental value, for both a and c. The interplanar equilibrium distances at both pressures are in good agreement with experiment. Note that the positions of Te and I planes are inverted when compared to the experimental results [101]; this effect has been attributed to the fact that x-ray diffraction measurements cannot establish a clear distinction between the Te and I planes, due to their very similar atomic radii and charge [29]. This inversion is systematic throughout the first-principles literature [29,53,72,75].
Analyzing the full pressure dependence of the lattice parameters for PBE-GGA and DFT-D3BJ between 0 and 5 GPa (not included here) lead to similar results as those of Güler-Kiliç and Kiliç (see Figs. 3(g) and 3(h) of [75]), where the overestimation of a remains around 2% throughout the whole pressure range, while for c it drops rapidly under 2% after 1 GPa. Thus, for P 1 GPa, the lattice structure is similarly well described by both methods. The resulting lattice parameters and interatomic distances alongĉ obtained by the optimization procedure are presented in Table II for all pressures included in this work.
The bulk band gap of BiTeI has been measured at ambient pressure via ARPES [25], optical spectroscopy [102] and optical spectral ellipsometry [33], yielding respectively 0.38, 0.36 and 0.33 eV. Another soft x-ray ARPES experiment obtained a smaller value of 0.26 eV [103]; they however argued that their measurement probes a region between the bulk and the subsurface, where some band bending effects could lower the measured band gap compared to the bulk value. Our results with PBE, 0.30 eV for the indirect gap and 0.33 eV for the direct gap, are in pleasantly good agreement with experiment, thanks to the overestimation of the lattice parameters counterbalancing PBE's well-know gap underestimation [42]. For DFT-D3BJ, the band gap is strongly underestimated: the system is almost metallic, with a band gap of about 30 meV at ambient pressure. We are not aware of any reference for experimental gap measurements in the topological phase.
Beyond the crystal structure, the parameters characterizing the Rashba splitting can also be used to assess the level of agreement of the calculated electronic structure with experiment. The Rashba energy, E R , is defined as the energy difference between electronic states at the band extrema and at the time-reversal protected  Fig. 3]. The momentum offset between those two k-points is captured by the Rashba wavevector, k R . The Rashba parameter, is the strength of the spin-orbit interaction in the Rashba Hamiltonian. These parameters are clearly displayed in the left panel of Fig. 2a and Eq. (1) of Monserrat and Vanderbilt [53]. Our calculations agree reasonably well with available experimental data [see Table III], with a clear advantage for PBE-GGA.

Vibrational properties
At the zone center, Γ, BiTeI has C 3v point group symmetry, yielding four Raman-active modes, whose irreducible representations are two doubly-degenerate E modes and two A 1 modes [32,104]. Since BiTeI has no inversion center, all modes are both Raman and infrared active [32]. Phonon frequencies for these modes are compared to available experimental data and previous first-principles calculations in Tables IV and V. At ambient pressure, we obtain a relative mean deviation with experiment of 4.1% for PBE and 1.6% for D3BJ. At 5 GPa, these deviations drop to 3.8% for PBE and 0.8% for D3BJ. We verified that the general pressure dependence of the Raman-active frequencies agreed well with experiment in both cases. Our calculated values for LO-TO splitting are in excellent agreement with previous calculations [72]. A small feature measured at ω=118.5 cm −1 at ambient pressure [72] can be associated with the LO-TO split E 2 mode, thus validating our results.

Discussion
The purpose of this work is to analyze the temperature dependence of the TPT in BiTeI. The numerical value of the critical pressures is directly related to the bare band gap value, which will be corrected by electron-phonon interaction. From Eq. (5) and (11), we can see that the phonon frequencies and the bare electronic eigenvalues are key physical quantities in this correction. This dependence goes beyond the simple electronic dispersion: as discussed in the main text, the band gap also plays a crucial role. Indeed, the denominator captures the "unlikelyhood" of couplings between two eigenstates with a large energy difference, rooted in perturbation theory. Furthermore, considering the renormalization of both band extrema forming the fundamental band gap, the gap energy will act as a "barrier", giving a lower bound on the smallest value of the denominator for couplings with eigenstates within the subset of bands located across the gap. Therefore, a significant underestimation of the bare band gap will artificially emphasize interband couplings and bias the resulting zero-point renormalization.
While one could argue that this effect should be small at T = 0 K, especially for a material with heavy atoms like BiTeI, it should not be forgotten that it will be amplified by the phonon occupation factor at higher temperatures, and thus could lead to significantly different conclusions. Lastly, since we are dealing with a TPT, the orbital character of the valence and conduction bands close to their extrema is subject to an inversion process. This phenomenon is captured by the wavefunction, such that a wrong band character could alter the different matrix elements entering the self-energy. While DFT-D3BJ provides the most accurate descriptions of the cell geometry and Raman-active phonon frequencies, our results indicate that in our case, it is not a suitable choice to track the temperature dependence of the TPT. On the one hand, despite its accurate description of Van der Waals bonding between the atomic planes, DFT-D3BJ delivers a severe underestimation of the band gap at ambient pressure. This can be understood by considering the van der Waals correction as acting like an "effective pressure" along theĉ axis. As previously argued, such a small value of the band gap would not only flaw our bare P C predictions (by the same arguments as Rusinov et al. [42]) but could also artificially strengthen some interband couplings in the EPI self-energy.
One could easily argue that this band gap underestimation could be overturned simply by applying a scissor shift operator to the eigenvalues entering the self-energy, or be compensated by evaluating the pressure shift that would bring the gap to the experimental value. We considered both these approaches. For the latter, we estimate that one must shift towards a negative pressure P 0 −2 GPa in order to bring the DFT-D3BJ gap close to the experimental value, thus bringing the lattice parameters very close to the PBE relaxed values at ambient pressure.
About the former, we have investigated the temperature-dependent renormalization at 0, 1.5, 3.5 and 5 GPa, with and without the scissor shift operator. The results for P = 1.5 GPa shown in Fig. 9 are quite striking: PBE-GGA (red circles) and DFT-D3BJ with (green diamonds) and without (blue squares) the scissor shift disagree about the sign of the renormalization for the CBM. This sign disagreement is also reflected in the total gap renormalization.
This discrepancy brings to light a fundamental aspect of the calculation, namely the fact that the EPI renormalization is affected by the system's topology. From Z 2 topological invariant calculations, we found that at 1.5 GPa, the DFT-D3BJ wavefunction already is in the topological state (yielding a negative correction), while for PBE-GGA it is still in the trivial state (where the correction is positive). These are precisely the trends observed in Fig. 5. Moreover, if we compare the DFT-D3BJ without scissor shift (blue squares) to the PBE-GGA result at 3.5 GPa (dashed yellow line, triangle markers), which is roughly P = 1.5 GPa−P 0 , we can see that those two calculations, with both wavefunctions in the topological state and with roughly the same lattice parameters and bare band gap, deliver almost identical results. Thus, even if DFT-D3BJ delivers the most accurate lattice parameters and Raman-active phonon frequencies at 1.5 GPa when compared to experiment, it predicts a Z 2 topological index that contradicts all experimental evidences of the phase transition. Therefore, computing the renormalization with DFT-D3BJ and a scissor shift operator for a given pressure would roughly result in applying the electron-phonon contribution from a PBE-GGA calculation at another pressure, shifted by P 0 . This is especially crucial when investigating pressures neighboring the TPT, as it can lead to wrongly assign a renormalization computed with a TI wavefunction to a shifted band structure intended to describe the trivial phase. For pressures sufficiently far from the transition in the TI phase, we expect the trend of a DFT-D3BJ with scissor shift calculation to be correct, despite a discrepancy in the numerical value of the renormalization. For the trivial phase, a Z 2 calculation showed that 0 < P C1 < 0.5 GPa for DFT-D3BJ. Therefore, we cannot rely on such results to describe the ambient pressure renormalization since the resulting wavefunction does not predict the correct band character.
For these reasons, we chose to rely on PBE-GGA functional: despite a clear overestimation of the lattice parameters, it does deliver a band gap, Rashba splitting and Raman-active phonon frequencies that are in reasonable overall agreement with experiment. Moreover and most crucially, it exhibits the right band characters at ambient pressure and predicts a critical pressure P C1 in agreement with experimental evidence.
Appendix B: Pressure and temperature-dependent Rashba Splitting As discussed in the main text, one of the assumptions underlying our work is that the Weyl nodes trajectory throughout the WSM phase is not qualitatively affected by temperature. In other words, we suppose that the k z location of the P C1 and P C2 planes remains constant; hence we only investigated EPI and TE within these planes. Without having access to the pressure and temperature dependence of the band structure in the full BZ, we can nevertheless gain significant insight by tracking the pressure and temperature dependence of the different parameters characterizing Rashba splitting within these planes. The Rashba energy, E R , captures the steepness of the band's slope in the vicinity of the band gap. The Rashba momentum, k R , tracks the distance between the band extrema and the degeneracy point at (0, 0, k z ), that is, the temperature-dependent band gap location within a given k z plane. Finally, the Rashba parameter, α R , captures the strength of the Rashba interaction (Eq. (A1)), which is related to the potential gradient along theĉ axis [105]. Fig. 10 presents the results for the Rashba energy (a), Rashba momentum (b) and Rashba parameter (c), for both the valence (bottom panels) and conduction bands (upper panels). We emphasize that these parameters were computed for the minimal band gap, namely on the (H-A, P C1 ) line for the trivial phase and on the (A-L, P C2 ) for the TI phase. Hence, the numerical values for ambient pressure in this figure should not be compared to experimental data, as available ARPES measurements were done on the (A-L, P C1 ) line. The trends reported here nevertheless agree with those of Monserrat and Vanderbilt [53] for the (A-L, k z = π c ) direction. As can be seen in Fig. 10, pressure has a stronger effect than temperature on the Rashba parameters. This behavior can be qualitatively understood as increasing temperature enhances atomic vibrations around their slightly shifted equilibrium positions, while pressure modifies more significantly these equilibrium positions by reducing the distance between the atomic planes. Hence, for a given temperature, increasing pressure enhances the potential gradient along theĉ axis, thus increasing α R . On the contrary, for a given pressure, temperature increases the complex interplay between the electrons and the lattice, yielding a reduction of α R . These are the behaviors observed in the trivial phase [ Fig. 10(c), diamond mark-  ⇐ Static WSM (c) α R (eV·Å) P C1 plane P C2 plane P C1 plane P C2 plane Pressure and temperature-dependent Rashba parameters, for conduction (top panels) and valence bands (bottom panels), near the minimum band gap of the PC1 (diamond markers) and PC2 (round markers) planes, as defined in Fig. 2(b). Note that for P < PC1 these quantities were computed in the H-A direction in the PC1 plane, hence the difference with the data in Table III for the conduction band at ambient pressure, where it was computed in the A-L direction for comparison with experimental data. For P > PC2 it was computed in the A-L direction in the PC2 plane.

ers].
The TI phase displays the opposite behavior [ Fig. 10(c), round markers]: α R is reduced by pressure and increased by temperature. To understand this behavior, we must analyze more carefully both the Rashba energy and momentum, as α R captures the interplay between those two quantities. On the one hand, for all pressures, the Rashba energy [ Fig. 10(a)] almost steadily increases with temperature, except for the CB in the trivial phase, which shows a very weak decrease with temperature. On the other hand, the Rashba momentum is almost unaffected by pressure and temperature in the TI phase [ Fig. 10(b), round markers], while it increases almost steadily with temperature for pressures below the TPT.
The opposite trends of α R between the trivial and TI phase can, therefore, be understood by two observations. On the one hand, in both phases, the pressure dependence of α R is governed by the Rashba energy, which increases in the trivial phase and decreases in the TI phase. On the other, their opposite temperature dependence is caused by the Rashba momentum. In the trivial phase, both E R and k R increase with temperature, the latter being more significant than the former, hence the decrease of α R . In the TI phase, the increase of α R captures the change of the Rashba energy, as k R is scarcely affected by temperature. We cannot, however, determine if this behavior is linked to the TPT or merely a consequence of an increased bonding alongĉ induced by increased pressure. This question could be addressed by investigating the temperature dependence of the high-pressure Rashba splitting in other bismuth tellurohalides, BiTeBr and BiTeCl, which are not known to exhibit a pressureinduced TPT [42].
Lastly, we verified that the trends presented in Fig. 10 for the TI phase do not vary if the temperaturedependent band structure is evaluated in the P C1 plane rather than the P C2 plane, remaining in the A-L direction. Hence, it seems reasonable that these trends, including the quasi temperature-independence of the Rashba momentum, could be applied to any k z between those planes. We also verified that the temperaturedependent values of P C2 are not significantly altered if considering the EPI and TE corrections computed in the P C1 plane.
Appendix C: Temperature-dependent Weyl semimetal phase width In Table VI, we report the numerical values for the temperature-dependent WSM phase boundaries, namely P C1 (T ) and P C2 (T ). These values were extracted from Fig. 7 using the extrapolation procedure described in VI. Temperature-dependent WSM phase boundaries, obtained with the extrapolation procedure described in Sect. IV C. PC1 marks the phase transition from a trivial insulator to a WSM, while PC2 marks the transition from WSM to Z2 topological insulator. The EPI+TE values are plotted in Fig. 8  Sect. IV C and are explicitly plotted in Fig. 8. For comparison, we also include the temperature-dependent crit-ical pressures obtained by considering the EPI and TE corrections independently.