# First-principles study of lattice dynamical properties of the room-
temperature $P2_{1}/n$ and ground-state $P2_{1}/c$ phases of WO3

Hamideh Hassani  Physique Théorique des Matériaux, QMAT, CESAM, Université de
Lìege, B-4000 Sart-Tilman, Belgium  Department of Physics, University of
Antwerp, Groenenborgerlaan 171, 2020 Antwerp, Belgium  Bart Partoens
Department of Physics, University of Antwerp, Groenenborgerlaan 171, 2020
Antwerp, Belgium  Eric Bousquet  Physique Théorique des Matériaux, QMAT,
CESAM, Université de Lìege, B-4000 Sart-Tilman, Belgium  Philippe Ghosez
Physique Théorique des Matériaux, QMAT, CESAM, Université de Lìege, B-4000
Sart-Tilman, Belgium

###### Abstract

Using first-principles density functional theory, we investigate the dynamical
properties of the room-temperature $P2_{1}/n$ and ground-state $P2_{1}/c$
phases of WO3. As a preliminary step, we assess the validity of various
standard and hybrid functionals, concluding that the best description is
achieved with the B1-WC hybrid functional while a reliable description can
also be provided using the standard LDA functional. We also carefully
rediscuss the structure and energetics of all experimentally observed and few
hypothetical metastable phases in order to provide deeper insight into the
unusual phase diagram of WO3. Then, we provide a comprehensive theoretical
study of the lattice dynamical properties of the $P2_{1}/n$ and $P2_{1}/c$
phases, reporting zone-center phonons, infrared and Raman spectra as well as
the full phonon dispersion curves, which attest for the dynamical stability of
both phases. We carefully discuss the spectra, explaining the physical origin
of their main features and evolution from one phase to another. We reveal a
systematic connection between the dynamical and structural properties of WO3,
highlighting that the number of peaks in the high-frequency range of the Raman
spectrum appears as a fingerprint of the number of antipolar distortions that
are present in the structure and a practical way to discriminate between the
different phases.

##  I Introduction

Tungsten trioxide, WO3, is a wide bandgap semiconductor transition metal oxide
with peculiar electronic, chromic and optical properties that makes it
appealing for photocatalytic, electrochromic and optoelectrical devices
applications [1, 2]. WO3 can be classified as a perovskite-like ABO3 compound
but with an empty $A$ site. As such, it possesses an aristotype $Pm\bar{3}m$
cubic structure, in which W atoms sit at the center of oxygen octahedra
sharing a corner and forming a continuous network. However, WO3 is typically
not observed in that reference cubic structure and experimental measurements
at various temperatures reveal a very complex phase diagram, not fully
elucidated yet [3]. At about 1300 K, WO3 is reported to crystallize in a
$P4/nmm$ tetragonal phase [3, 4, 5]. Then, on cooling, it shows consecutive
phase transitions to (i) a $P4/ncc$ tetragonal phase at about 1150 K [3, 4, 6,
5], (ii) an eventual $P2_{1}/c$ monoclinic phase at about 1050 K [3], (ii) a
$Pbcn$ orthorhombic phase at about 1000 K [6, 3, 4, 7], (iii) a $P2_{1}/n$
monoclinic phase at about 600 K [3, 4, 6, 8], (iv) a $P\overline{1}$ triclinic
phase at 273 K [9, 4, 10] and, ultimately, (v) a $Pc$ monoclinic polar phase
at around 200 K [11, 4, 12, 13] that was recently proposed to be instead a
related non-polar $P2_{1}/c$ monoclinic phase [14]. All these phases appear as
small distortions of the cubic artistotype structure but including distinct
rotations and tilts of the oxygen octahedra and different shifts of W atoms
within these octahedra.

Experimentally, WO3 appears also in practice as a sub-stoichiometric (WO3-x)
compound. This, on the one hand, confers it some unique properties exploited
in a variety of technological applications [15] but, on the other hand, makes
its study significantly more complex. In fact, the off-stoichiometry makes it
a $n$-doped semiconductor and the excess electrons interact with the
polarizable crystal lattice in order to form polarons and bipolarons [16] that
strongly influence its properties.

Although extensive studies have been reported concerning the electronic and
optical properties of WO3 [17, 18, 19, 15], there is a much more limited
amount of works dedicated to the investigation of its lattice dynamical
properties. The phonon dispersion curves of cubic WO3 and its elastic
properties have been studied by Fan et al. [20] or together with the electron-
phonon coupling by Mascello et al. [21]. Some lattice dynamical properties of
distinct WO3 phases have also been reported by Yang et al. [22]. Furthermore,
Hamdi et al. [14] studied the phonon instabilities of the cubic phase together
with the various potential metastable phases originating from these
instabilities. This latter study has shown that the B1 Wu-Cohen (B1-WC) hybrid
exchange-correlation functional [23] can accurately reproduce the experimental
results regarding both the electronic and structural properties of WO3.

In the present work, we first carefully re-discuss the phase diagram of WO3,
assessing the validity of distinct functionals. We then provide a
comprehensive theoretical study of the lattice dynamical properties of the
$P2_{1}/n$ room temperature phase and $P2_{1}/c$ expected ground-state phase,
reporting infrared and Raman spectra together with the full phonon dispersion
curves. On the one hand, we carefully discuss the spectra, explaining the
physical origin of their main features and evolution from one phase to another
and providing meaningful benchmark results for the interpretation of
experimental measurements. On the other hand, polaron formation is directly
linked to electron-phonon interaction so that the present understanding of the
phonon properties appears as a useful step toward the understanding of
polarons in these two phases.

The paper is organized as follows. In the first part, we compare the
structural parameters and energetics of various phases of WO3, as obtained
using standard (LDA and GGA) and more advanced (B1-WC, HSE06) hybrid
functionals. We also take this opportunity to re-discuss briefly the phase
diagram of WO3. In the second part, we carefully discuss the lattice dynamical
properties of the $P2_{1}/n$ room temperature and $P2_{1}/c$ ground state
structures, comparing our theoretical results with experimental data when
available and interpreting the shape of the spectra. From this analysis, Raman
measurements emerges as a promising tool to identify the various phases of
WO3.

##  II Computational details

Our calculations are performed in the framework of density functional theory
(DFT) as implemented in CRYSTAL17 [24] based on a linear combination of
localized basis functions and ABINIT [25, 26] making use a plane-wave basis
set. Calculations are reported using distinct approximations of the exchange-
correlation energy including the local density approximation (LDA), the
generalized gradient approximation (GGA) and different hybrid schemes.
Integrations over the Brillouin zone are approximated by sums over a mesh of
$8\times 8\times 8$ $k$-points for the cubic phase or meshes providing a
comparable sampling for other phases (meshes of $6\times 6\times 8$, $6\times
6\times 4$ and $4\times 4\times 4$ are used for P4/nmm, P4/ncc and Pbcn phases
respectively).

ABINIT is used with the LDA with Perdew-Wang’s parametrization [27] and the
GGA with both the PBEsol [28] and Wu and Cohen (WC) functionals [29]. These
calculations make use of norm–conserving pseudopotentials from the Pseudo Dojo
Table v0.4 [30], considering the $2s$ and $2p$ orbitals of O and the $5s$,
$5p$, $5d$ and $6s$ orbitals of W as valence states. The electronic wave
functions are expanded in plane-waves up to an energy cutoff of 60 Ha. Full
relaxations of the lattice parameters and atomic positions are performed with
each functional using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm
[31] until the maximal forces and stresses are less than $10^{-6}$ Ha/Bohr and
$10^{-4}$ GPa, respectively. The electronic self-consistent calculations were
converged until the difference of the total energy is smaller than $10^{-9}$
Ha. Dynamical matrices, Born effective charges and dielectric tensors are
obtained within a variational approach [32] to density functional perturbation
theory (DFPT) [33].

CRYSTAL is used for the hybrid functional calculations with both the B1-WC
functional [34] and the HSE06 functional [35]. Benchmark calculations in LDA
have also been performed for comparison with ABINIT. We use all-electron
double-$\zeta$ basis sets for oxygen whereas, for the heavy tungsten atom, we
use the effective core pseudopotential (ECP) technique as implemented in the
small core Detlev-Figgen pseudopotentials associated with correlation
consistent polarized Valence Triple-$\zeta$ Basis Set (cc-pVDZ) [36] for
explicit treatment of valence electrons. The truncation thresholds in the
evaluation of Coulomb and exchange series appearing in the SCF equation for
periodic systems are adjusted to ${10}^{-7}$ for Coulomb overlap tolerance,
${10}^{-7}$ for Coulomb penetration tolerance, ${10}^{-7}$ for exchange
overlap tolerance, ${10}^{-7}$ for exchange pseudo-overlap in the direct
space, and ${10}^{-14}$ for exchange pseudo-overlap in the reciprocal space.
The tolerance on change in total energy for the SCF convergence is adjusted to
${10}^{-10}$ Ha. Full structural relaxations are carried out by using a quasi-
Newton algorithm with a BFGS Hessian updating scheme so that a convergence of
less than $5\times 10^{-5}$ Ha/Bohr and $10^{-3}$ Bohr is reached respectively
on the root mean square of the gradient and displacements. The phonon
frequencies are computed using the frozen phonon numerical differences [37,
38] (for detailed information check Fig. 1(a) in the appendix). The Born
effective charges and the optical dielectric tensors are evaluated via the
Berry phase technique [39] and the Coupled-Perturbed Kohn-Sham approach [40,
41], respectively.

##  III Structural properties

Figure 1:  The phonon dispersion curve of the cubic phase of WO3, computed
with the B1-WC hybrid functional and proper treatment of the dipole-dipole
interactions. Lines are colored according to the involvement of W (red) and O
(blues) atoms in each mode.

We begin our study by considering the high-symmetric parent cubic structure
($Pm\overline{3}m$) of WO3. Fig. 1 shows its phonon dispersion curves
calculated using the B1-WC hybrid functional. The high-frequency branches show
slightly distinct dispersion than those previously reported in [14]. The
previous results were unfortunately affected by a bug in the implementation of
the Fourier interpolations technique in the CRYSTAL package while the new
curves reported here rely on a proper treatment of the dipole-dipole
interactions as implemented in ABINIT (for further information we refer to
Fig. 1(a) in the appendix).

Table 1: Experimental and calculated lattice parameters (Å) of different distorted phases of WO3, calculated by LDA, GGA-WC, PBEsol, HSE06 and B1-WC hybrid functionals. Also the PBE results of Yang et al. [22] are included. Data in parentheses represent relative deviations from the experimental values in percent.  Phase | Exp | HSE06 $(\Delta\%)$ | B1-WC $(\Delta\%)$ | GGA-WC $(\Delta\%)$ | PBEsol $(\Delta\%)$ | PBE [22] $(\Delta\%)$ | LDA $(\Delta\%)$
---|---|---|---|---|---|---|---
| 5.29111Ref [3] | 5.29 (+0.0) | 5.29 (+0.0) | 5.31 (+0.3) | 5.30 (+0.1) | 5.37 (+1.5) | 5.27 (-0.3)
P4/nmm | 3.92 | 3.96 (+1.0) | 3.92 (+0.0) | 3.92 (+0.0) | 3.92 (+0.0) | 4.05 (+3.3) | 3.87 (-1.1)
| 5.27a | 5.20 (-1.3) | 5.16 (-2.0) | 5.26 (-0.1) | 5.24 (-0.5) | - | 5.17 (-1.8)
P4/ncc | 7.84 | 7.95 (+1.4) | 7.87 (+0.3) | 7.85 (+0.1) | 7.85 (+0.1) | - | 7.74 (-1.2)
| 5.27a | 5.28 (+0.1) | 5.26 (-0.1) | 5.33 (+1.1) | 5.27 (+0.0) | - | 5.23 (-0.7)
$P2_{1}/c$-HT | 5.26 | 5.21 (-0.9) | 5.15 (-2.0) | 5.31 (+0.9) | 5.27 (+0.1) | - | 5.16 (-1.9)
| 7.83 | 7.70 (-1.6) | 7.61 (-2.8) | 7.72 (-1.4) | 7.72 (-1.4) | - | 7.62 (-2.7)
| 7.33222Ref [6] | 7.46 (+1.7) | 7.28 (-0.6) | 7.45 (+1.6) | 7.42 (+1.2) | 7.46 (+1.7) | 7.36 (+0.4)
Pbcn | 7.57 | 7.50 (-0.9) | 7.52 (-0.6) | 7.62 (+0.6) | 7.59 (+0.2) | 7.61 (+0.5) | 7.39 (-2.3)
| 7.74 | 7.69 (-0.6) | 7.68 (-0.7) | 7.76 (+0.2) | 7.76 (+0.2) | 7.98 (+3.1) | 7.62 (-1.5)
| 7.30a | 7.43 (+1.7) | 7.37 (+0.9) | 7.52 (+3.0) | 7.51 (+2.8) | 7.49 (+2.6) | 7.37 (+0.9)
$P2_{1}/n$ | 7.53 | 7.55 (+0.2) | 7.47 (-0.7) | 7.61 (+1.0) | 7.59 (+0.7) | 7.54 (+0.1) | 7.39 (-1.8)
| 7.69 | 7.62 (-0.9) | 7.54 (-1.9) | 7.67 (-0.2) | 7.66 (-0.3) | 8.02 (+4.2) | 7.59 (-1.3)
| 7.30333Ref [42] | 7.36 (+0.8) | 7.33 (+0.4) | 7.48 (+2.4) | 7.44 (+1.9) | 7.50 (+2.7) | 7.32 (+0.2)
$P\overline{1}$ | 7.52 | 7.54 (+0.2) | 7.44 (-1.0) | 7.62 (+1.3) | 7.59 (+0.9) | 7.60 (+1.0) | 7.43 (-1.1)
| 7.67 | 7.70 (+0.3) | 7.61 (-0.7) | 7.72 (+0.6) | 7.72 (+0.6) | 7.98 (+3.8) | 7.63 (-0.5)
| 5.27 444Ref [11] | 5.28 (+0.1) | 5.26 (-0.1) | 5.33 (+1.1) | 5.27 (+0.0) | - | 5.23 (-0.7)
$P2_{1}/c$-LT | 5.15 | 5.21 (+1.1) | 5.15 (+0.0) | 5.31 (+3.1) | 5.27 (+2.3) | - | 5.16 (+0.1)
| 7.66 | 7.70 (+0.5) | 7.61 (-0.6) | 7.72 (+0.7) | 7.72 (+0.7) | - | 7.62 (-0.5)
Table 2: Lattice parameters (Å) of metastable distorted phases of WO3 not observed experimentally, calculated using the LDA, GGA-WC, PBEsol, HSE06 and B1-WC hybrid functionals. Dash symbols correspond to cases for which the phase is not metastable and the system relaxes back to a higher symmetry. Phase | HSE06 | B1-WC | GGA-WC | PBEsol | LDA
---|---|---|---|---|---
$Pm\overline{3}m$ | 3.78 | 3.78 | 3.79 | 3.79 | 3.76
| 5.30 | 5.26 | - | 5.35 | 5.25
P4/mbm | 3.79 | 3.78 | - | 3.79 | 3.76
| 3.88 | 3.85 | 3.85 | 3.85 | 3.81
$P2_{1}/m$ | 7.43 | 7.44 | 7.46 | 7.45 | 7.41
| 3.88 | 3.85 | 3.85 | 3.85 | 3.81
$R\overline{3}c$ | 5.31 | 5.35 | 5.36 | 5.35 | 5.29
R3m | 3.82 | 3.81 | 3.81 | 3.81 | 3.78
R3c | 5.39 | 5.35 | - | - | 5.32
| 5.32 | 5.32 | - | - | 5.30
Pnma | 7.44 | 7.42 | - | - | 7.39
| 5.23 | 5.22 | - | - | 5.20

Imaginary frequencies appear as negative numbers in Fig. 1 and are associated
to phonon structural instabilities of the parent cubic phase. In fact, the
distinct phases of WO3, observed at different temperatures, arise from the
individual or combined condensation of the cooperative atomic motions related
to these unstable phonon modes. The latter include (i) polar modes (PM) at
$\Gamma$ (corresponding to irreducible representation $\Gamma^{-}_{4}$), which
are related to the motion of W against the O atoms, (ii) anti-polar modes (AP)
at the M and X points (M${}^{-}_{3}$, X${}^{-}_{5}$), which correspond to
opposite displacements W against the O atoms from unit cell to unit cell along
the [110] and [100] directions, respectively, leading to the displacement of
the W atoms from the center of the octahedral and forming chains with
alternating short and long W–O bond lengths, and finally (iii)
antiferrodistortive modes (AFD) at the M and R points (M${}^{+}_{3}$ and
R${}^{+}_{4}$), which consist in rotation and tilts of the oxygen cages about
the central W atom.

Phonon dispersion curves obtained with different functionals are compared in
the appendix. Although they appear very similar, we point out that the
M${}^{+}_{3}$ AFD mode is slightly unstable using hybrid and LDA functionals
while it is slightly stable using GGA functionals.

In what follows, we will consider all the experimentally observed phases as
well as six additional metastable phases arising all from distinct
combinations of the unstable modes to gain more insight into the energy
landscape of WO3. Full structural optimizations of the lattice parameters and
atomic positions have then been performed for each phase using different
exchange-correlation functionals.

###  III.1 Lattice parameters

In Table 1, we compare to experiment the calculated lattice parameters
obtained with the HSE06 and B1-WC hybrid functionals, the WC and PBEsol GGA
functionals, the LDA functional as well as the results relying on the PBE GGA
functional as reported recently by Yang et al. [22]. In line with what was
previously reported by Hamdi et al. [14], starting from the experimental
structure with the $Pc$ space group and performing full structural
relaxations, the system relaxes back to a $P2_{1}/c$ phase of higher symmetry
using all considered functionals. Although a $P2_{1}/c$ phase has only been
experimentally observed at about 1000 K [3], its structural parameters are
very similar to those of the $Pc$ phase observed at about 200 K [11] (the only
difference is a small polar distortion as further discussed in the next
Section). In Table 1, the lattice parameters of the relaxed $P2_{1}/c$ phase
are therefore compared both to those of the high-temperature $P2_{1}/c$ phase
($P2_{1}/c$-HT) and low-temperature $P_{c}$ phase to which we also refer to as
$P2_{1}/c$-LT hereafter.

For the high-temperature $P4/nmm$ phase, all functionals provide a rather fair
agreement with experimental data although PBE results of Yang et al. show a
significant overestimation (also for other phases). Both GGA-WC and PBEsol
provide also a good prediction of the lattice constants for the $P4/ncc$
phase. It should be noticed however that the experiments on these high-T
phases are performed above 1000 K [3] while the calculations correspond to 0
K. It is therefore possible that this good agreement arises from a
cancellation of errors coming, on the one hand, from neglecting the thermal
expansion and, on the other hand, from a natural propensity of these GGA
functionals to overestimate the lattice parameters of WO3. In line with that
suggestion, the GGA-WC and PBEsol are indeed in better agreement with the
$P2_{1}/c$-HT phase and their overestimate of lattice parameters tends to
increase for the low-temperature phases. On the contrary, the LDA tends to
systematically underestimate the lattice constants as expected and provides a
good agreement for the $P2_{1}/c$-LT phase. HSE06 provides good estimates
although with a tendency to slightly overestimate the volume of the low-
temperature phases. B1-WC is confirmed to be the one providing the best
overall estimates for all phases, and gives good description of the
$P2_{1}/c$-LT and $P2_{1}/n$ phases under investigation here. The computed
lattice parameters of few additional metastable phases, not observed
experimentally, are reported in Table 2.

Figure 2: Symmetry adapted mode decomposition of distorted WO3 phases that are
observed experimentally. Comparison between B1-WC, HSE06, GGA-WC, PBEsol, LDA
and experimental results [3, 6, 11]. The modes with an amplitude lower than
0.04 Angstrom are not shown. Figure 3: The symmetry adapted mode decomposition
of distorted WO3 phases which are not observed experimentally. Comparison of
B1-WC, HSE06, GGA-WC, PBEsol and LDA. The modes with an amplitude lower than
0.04 Angstrom are not shown. The black lines correspond to the cases where the
phases are not condensed.

###  III.2 Atomic distortion

In this section, we analyze further the structural distortions of each phase
with respect to the reference cubic phase by projecting the internal atomic
displacements on symmetry adapted modes of the parent cubic phase using the
AMPLIMODE software package [43]. The results are presented in Fig. 2 for the
experimentally observed phases and in Fig. 3 for the few other phases. For the
latter, missing phases in GGAs correspond to cases for which the phase is not
metastable and the system relaxes back to a higher symmetry.

As highlighted in Fig. 2, the $P4/nmm$ phase arises from the condensation of
the only anti-polar $M^{-}_{3}$ mode along the $z$ direction. This mode will
remain present in all experimentally observed phases with a rather constant
amplitude. The $P4/ncc$ phase includes the additional condensation of the AFD
R${}^{+}_{4}$ mode along the same direction yielding a tilt pattern
($a^{0}a^{0}c^{-}$) in Glazer’s notations. The $P2_{1}/c$ phase arises from
the addition of AFD R${}^{+}_{4}$ contributions, yielding a tilt pattern
($a^{-}a^{-}c^{-}$). A secondary antipolar $X_{5}^{-}$ mode in the $xy$-plane
also appears thanks to an anharmonic coupling between the R${}^{+}_{4}$ and
M${}^{-}_{3}$ modes [44, 45, 46]. In the $Pbcn$ phase, a M${}^{+}_{3}$ mode
emerges giving rise to a distinct tilt pattern $a^{0}b^{+}c^{-}$. The
secondary antipolar $X_{5}^{-}$ mode in this phase is only along the $y$
direction. Also, another secondary antipolar $X_{5}^{+}$ mode along $y$
appears through a coupling between the R${}^{+}_{4}$ and M${}^{+}_{3}$ modes.
Besides, the antipolar $M^{-}_{3}$ mode is along the $z$ direction but with a
very slight tilt toward the $x$ direction. This small $x$ component of the
$M^{-}_{3}$ mode appears through a coupling between the secondary $X_{5}^{+}$
and $X_{5}^{-}$ modes and the primary R${}^{+}_{4}$ mode. The $P2_{1}/n$ phase
involves the same modes but combined differently within a tilt pattern
$a^{-}b^{+}c^{-}$. Again, the antipolar $M^{-}_{3}$ mode is along $z$ with a
very small $x$ component. Also, the antipolar $X_{5}^{-}$ and $X_{5}^{+}$
modes are along $y$ direction but with a very small $z$ component. The
$P\bar{1}$ phase shows an out-of-phase octahedral tilting around the second
axis yielding to tilt pattern $a^{-}b^{-}c^{-}$. This phase has the same
antipolar $M_{3}^{-}$ and $X_{5}^{+}$ modes as in the $P2_{1}/n$ phase, but
the antipolar $X_{5}^{-}$ mode has components in all directions. Finally, the
experimental $Pc$ phase appears similar to the $P2_{1}/c$ phase except for the
additional appearance of a polar $\Gamma_{4}^{-}$ mode. As previously
discussed, our calculations in $Pc$ symmetry always relax back to $P2_{1}/c$
to which the $Pc$ phase is therefore compared. We insist that all the phases
except $Pc$ are non-polar. Yang et al. [22] report a net dipole moment for the
$P2_{1}/n$ phase which is however inconsistent with the symmetry assignment.
Fig. 3 also compares the results for few additional metastable phases not
observed experimentally.

From Fig. 2 and Fig. 3, it appears that, on the one hand, all functionals
predict rather consistently the appearance and amplitudes of polar and
antipolar distortions ($M_{3}^{-}$, $X_{5}^{-}$, $X_{5}^{+}$,
$\Gamma_{4}^{-}$) while, on the other hand, GGA functionals strongly
underestimate the amplitudes of rotations and tilts ($R_{4}^{+}$,
$M_{3}^{+}$), as further discussed in the next Section. In line with that,
GGAs do not reproduce some metastable phases and also produce the worst
description of the $P2_{1}/n$ and and $P2_{1}/c$ phases. Hybrid functionals
typically provide a very good description of all phases while LDA gives also
reasonable results, significantly better than GGAs.

Figure 4: Calculated energy gain with respect to the cubic phase of different
phases of WO3. Comparison between B1-WC, HSE06, GGA-WC, PBEsol and LDA.

###  III.3 Energetics

Fig. 4 compares the energies of the various phases under study in comparison
to that of the parent cubic phase taken as zero energy reference. The results
are shown for distinct DFT functionals. The energy is reported in meV per
four-atoms formula unit (meV/f.u.).

A first observation is that, globally, the gains of energy are slightly larger
using the hybrid functionals than the LDA/GGA functionals. This is in line
with the slightly larger amplitudes of phonon unstabilities predicted with the
hybrid functionals (see the appendix). It is also consistent with the slightly
larger amplitudes of distortions obtained with the hybrid functionals,
although LDA shows energies more comparable to GGA but amplitudes of
distortion closer to hybrid functionals.

Comparing the lowerings of energy produced by the individual condensation of
distinct kinds of instabilities (AFD, PM, AP), we observe that the smallest
ones are related to the AFD motions ($P4/mbm$, $Pnma$, $R\bar{3}c$ phases).
This is compatible with the fact that AFD modes are less unstable than FE or
AP modes. The energy lowerings produced by AFD modes are also significantly
smaller with GGA functionals that also show slightly smaller AFD instabilities
(the $M_{3}^{+}$ mode is even stable in GGA-WC). Consequently, some AFD phases
are not even appearing as metastable in GGA. We further notice that the energy
lowerings produced by AFD modes are also smaller with HSE06 than B1-WC
although in this case the phonon dispersion curves and the amplitudes of
distortion are relatively close. Only LDA shows an energetics of AFD
comparable to B1-WC.

Regarding the individual condensation of PM and AP modes, although the
strongest instability is observed in the dispersion curves for the
$\Gamma^{-}_{4}$ mode with all functionals, the largest gain of energy arises
from the condensation of the AP X${}^{-}_{5}$ ($P2_{1}/m$ phase) except with
the B1-WC for which the $R3m$ phase resulting from the condensation of
$\Gamma^{-}_{4}$ is lower in energy. We notice also that while condensation of
the AP M${}_{3}^{-}$ mode ($P4/nmm$ phase) produces a lowering of energy
smaller than the polar $\Gamma^{-}_{4}$ ($R3m$ phase) with hybrid functionals,
the opposite is true for LDA and GGA functionals.

As previously discussed by Hamdi et al. [14], the lowest energy phases of WO3
result from the combination of distinct instabilities and anharmonic couplings
play an important role in the formation of $Pbcn$, $P2_{1}/n$, $P\overline{1}$
and $P2_{1}/c$ phases. For this reason it is not straightforward to clarify
the origin of differences in the energy landscape obtained with different
functionals. As illustrated in Fig. 4, all functionals locate consistently
these four phases very close in energy (except for B1-WC that locates the
$Pbcn$ phase at slightly higher energy) but nevertheless in a different order
for each of them. We notice that only B1-WC and LDA predict the $P2_{1}/c$
phase as the ground state. Moreover, although not a requirement for internal
energies, the energetics of the different phases in B1-WC is in line with the
sequence of phase transitions observed experimentally.

All this highlights that WO3 exhibits a very delicate energy landscape.
Choosing the most appropriate functional remains a delicate issue that goes
beyond trying to best reproduce cell parameters and atomic distortions. Our
study confirms the B1-WC hybrid functional as the most appropriate choice and
points the LDA as a reasonable alternative.

Figure 5: The experimental symmetry adapted mode decomposition of observed
phases presented in the order they appear with decreasing temperatures.

###  III.4 Phase diagram

In line with its complex energy landscape, WO3 presents a complicated phase
diagram which, as mentioned before, is still partly debated. One question is
the potential polar nature of its ground state phase: as previously discussed
and confirmed by the present study, the $P2_{1}/c$ phase of stoichiometric WO3
is predicted as the theoretical ground state and does not show any tendency to
polar instability so that the experimental $P_{c}$ assignment, deviating only
marginally from the $P2_{1}/c$ prediction, might be related to extrinsic
effects [16]. Then, particularly puzzling is the intermediate $P2_{1}/c$ phase
appearing in a narrow temperature range at high-temperature and re-emerging as
the ground state at low temperature.

In Fig. 5, we compare the symmetry-mode analysis of all experimental phases
from high to low temperatures. The AP M${}_{3}^{-}$ is common to all phases
and remains roughly constant in amplitude. Then, on cooling, the AFD
R${}^{+}_{4}$ mode appears from the $P4/ncc$ phase and is preserved down to
the ground state. We notice that the Glazer tilt pattern associated to the
different phases (i.e. number and orientation of R${}^{+}_{4}$ modes) is
evolving from one phase to another. In the $Pbcn$ and $P2_{1}/n$ phases there
is the emergence of an additional AFD M${}_{3}^{+}$ distortion which however
compete with the R${}^{+}_{4}$ one [14] and disappears in the $P\bar{1}$ and
$P2_{1}/c$-LT phases.

The main difference between the high- and low-temperature $P2_{1}/c$ phases is
the amplitude AFD R${}^{+}_{4}$ mode. It might happen that at high-temperature
at which the R${}^{+}_{4}$ distortion are reduced, the M${}_{3}^{+}$
distortion can appear to produce intermediate phases while, due to mode
competition, the latter disappear when the R${}^{+}_{4}$ distortion becomes
larger at low temperatures.

Our calculations show that the $Pbcn$, $P2_{1}/n$, $P\overline{1}$ and
$P2_{1}/c$-LT phases are very close in energy. This might contribute to
explain the origin of the very complex and unusual phase diagram of WO3. It
suggests also that distinct phases could coexist and/or could be stabilized
under slight variation of the experimental conditions.

##  IV Born effective charges and dielectric constant

We consider now the Born effective charges ($Z^{*}$) and dielectric tensor
($\epsilon^{\infty}$). All results in the main text are obtained using the LDA
functional. Since the optimized LDA lattice parameters slightly underestimate
the experimental values, for better comparison with experiment, we decided to
work with phases for which (i) lattice parameters are fixed at the
experimental values for $P2_{1}/n$ phase and at the optimized B1-WC values for
$P2_{1}/c$ phase and (ii) atomic positions are fully relaxed under that
constraint. In practice, the deviations of atomic positions remain negligible
compared to the structures fully relaxed with B1-WC (see Tab. A1). Further
comparison of the $Z^{*}$ and $\epsilon^{\infty}$ computed in LDA and B1-WC
functionals are also provided in the appendix (Tables A2 and A3), showing
close agreement.

In Tab. 3 the Born effective charge tensors ($Z^{*}$) of WO3 are presented for
the eight nonequivalent atoms in the $P2_{1}/n$ phase, and the four
nonequivalent atoms in the $P2_{1}/c$ phase. Some $Z^{*}$ are anomalously
large in comparison to the nominal atomic charges ($+6e$, $-2e$ for W and O
atoms, respectively) and these values can be linked to the dynamical changes
in the hybridization between the $2p$ orbitals of the oxygen and the $5d$
orbitals of the W atoms. The main values of the symmetric part of these
tensors have also been reported for easier comparison of the different phases
and atoms.

In the cubic structure of WO3, the WO6-octahedra are undistorted and no
symmetry breaking is observed. This implies the existence of diagonal,
isotropic Born effective charge tensors: $+13.47e$ for W, $-9.79e$ for O||,
and $-1.8e$ for O⟂ consistent with those previously reported in Ref. [47]. By
distorting the cubic environment, the off-diagonal values of $Z^{*}$ increase
and the main values of $Z^{*}$ decrease with respect to the cubic phase, as
previously observed in perovskites [48, 49]. In both the $P2_{1}/c$ and
$P2_{1}/n$ phases the off-diagonal terms of $Z^{*}$ are small, with the only
exceptions the O1 and O2 atoms in the basal plane of the $P2_{1}/c$ crystal.
This is due to the fact that $Z^{*}$ is presented in the cartesian basis while
the unit cell basis of $P2_{1}/c$ is rotated by 45 degrees around the $z$
axis.

We believe that the significant differences between the diagonal components of
$Z^{*}$($W_{1}$) (and $Z^{*}$($W_{2}$)) for the $P2_{1}/n$ phase can be
associated with the distance between two consecutive tungsten atoms in W-O-W
chains, which are 3.70, 3.80 and 3.84 Å along the $x$, $y$ and $z$ directions,
respectively. Therefore, by decreasing the interatomic distance the diagonal
components of $Z^{*}$($W_{1}$) and $Z^{*}$($W_{2}$) will increase.
Furthermore, the same trend exists when the long-short bond alternation of W-O
is decreased. Although the same pattern is observed for the $P2_{1}/c$ phase,
the differences between the diagonal components of $Z^{*}$($W_{1}$) are
smaller than for the $P2_{1}/n$ phase.

Table 3:  Calculated (LDA) Born effective charge tensors ($Z^{*}$) and their main values (between parentheses) of the atoms of the asymmetric units (see notations in Tab. A1 in the appendix) of the $P2_{1}/n$ and $P2_{1}/c$ phases of WO3.  $P2_{1}/n$ |  |  |
---|---|---|---
$Z^{*}$($W_{1}$) |  $Z^{*}$($O_{1}$) |  $Z^{*}$($O_{3}$) |  $Z^{*}$($O_{5}$)
$\begin{pmatrix}10.79&-1.15&0.64\\\ 0.75&8.55&-0.37\\\ -0.49&-0.30&7.79\end{pmatrix}$ | $\begin{pmatrix}-8.00&0.93&-0.94\\\ 0.34&-1.25&0.08\\\ -0.44&0.10&-1.12\end{pmatrix}$ | $\begin{pmatrix}-1.43&0.09&-0.10\\\ 0.16&-6.22&-0.69\\\ 0.05&-0.28&-1.02\end{pmatrix}$ | $\begin{pmatrix}-1.35&0.03&-0.05\\\ -0.01&-1.05&-0.45\\\ -0.05&-0.14&-5.65\end{pmatrix}$
$\begin{pmatrix}10.81&8.65&7.66\\\ \end{pmatrix}$ | $\begin{pmatrix}-8.12&-1.19&-1.04\\\ \end{pmatrix}$ | $\begin{pmatrix}-6.26&-1.42&-0.97\\\ \end{pmatrix}$ | $\begin{pmatrix}-5.66&-1.34&-1.03\\\ \end{pmatrix}$
$Z^{*}$($W_{2}$) |  $Z^{*}$($O_{2}$) |  $Z^{*}$($O_{4}$) |  $Z^{*}$($O_{6}$)
$\begin{pmatrix}10.80&1.17&-0.78\\\ -0.78&8.62&-0.66\\\ 0.36&-0.04&7.82\end{pmatrix}$ | $\begin{pmatrix}-8.01&0.93&1.00\\\ 0.35&-1.25&-0.08\\\ 0.39&-0.10&-1.14\end{pmatrix}$ | $\begin{pmatrix}-1.43&-0.10&0.08\\\ -0.16&-6.28&-0.75\\\ -0.05&-0.27&-1.03\end{pmatrix}$ | $\begin{pmatrix}-1.35&0.01&0.13\\\ -0.03&-1.07&0.47\\\ 0.23&0.19&-5.63\end{pmatrix}$
$\begin{pmatrix}10.83&8.72&7.68\\\ \end{pmatrix}$ | $\begin{pmatrix}-8.14&-1.19&-1.06\\\ \end{pmatrix}$ | $\begin{pmatrix}-6.33&-1.39&-1.01\\\ \end{pmatrix}$ | $\begin{pmatrix}-5.59&-1.22&-1.18\\\ \end{pmatrix}$
$P2_{1}/c$ |  |  |
$Z^{*}$($W_{1}$) |  $Z^{*}$($O_{1}$) |  $Z^{*}$($O_{2}$) |  $Z^{*}$($O_{3}$)
$\begin{pmatrix}10.38&-1.07&-0.15\\\ 0.99&9.40&0.95\\\ -0.17&-0.89&7.74\end{pmatrix}$ | $\begin{pmatrix}-4.57&3.05&-0.66\\\ 2.92&-4.13&0.54\\\ -0.24&0.23&-1.06\end{pmatrix}$ | $\begin{pmatrix}-4.62&2.90&0.65\\\ 2.89&-4.01&-0.57\\\ 0.23&-0.23&-1.11\end{pmatrix}$ | $\begin{pmatrix}-1.19&0.00&0.16\\\ -0.03&-1.25&0.21\\\ 0.18&0.45&-5.56\end{pmatrix}$
$\begin{pmatrix}10.39&9.39&7.72\\\ \end{pmatrix}$ | $\begin{pmatrix}-7.46&-1.36&-1.01\\\ \end{pmatrix}$ | $\begin{pmatrix}-7.17&-1.40&-1.16\\\ \end{pmatrix}$ | $\begin{pmatrix}-5.66&-1.34&-1.04\\\ \end{pmatrix}$
Table 4: Optical dielectric tensor (${\varepsilon}^{\infty}_{ij}$) of the $P2_{1}/n$ and $P2_{1}/c$ phases of WO3 within LDA.  $P2_{1}/n$ | $P2_{1}/c$
---|---
$\begin{pmatrix}7.05&0&-0.02\\\ 0&5.83&0\\\ -0.02&0&5.35\end{pmatrix}$ | $\begin{pmatrix}6.95&0&-0.08\\\ 0&6.50&0\\\ -0.08&0&5.44\end{pmatrix}$
Table 5: Calculated (LDA) and experimental refractive indices ($n_{i}$) of the $P2_{1}/n$ and $P2_{1}/c$ phases of WO3 Phase | $n_{1}$ | $n_{2}$ | $n_{3}$
---|---|---|---
$P2_{1}/n$ (EXP) [50] | 2.70$\pm$0.035 | 2.37$\pm$0.035 | 2.28$\pm$0.035
$P2_{1}/n$ (LDA) | 2.65 | 2.41 | 2.31
$P2_{1}/c$ (LDA) | 2.63 | 2.54 | 2.33

The optical dielectric tensor (${\varepsilon}^{\infty}_{ij}$), and related
refractive indices ($n_{i}$) are reported in Tab. 4 and Tab. 5. The
experimental values [50] of the refractive indices measured at room
temperature are also reported for comparison. The results for the room
temperature $P2_{1}/n$ phase show a good agreement with the experimental data.
Our calculations indicate a high index of refraction for both the $P2_{1}/n$
and $P2_{1}/c$ phases. This is due to the large local field in these two
phases of WO3. Also, the birefringence ($\Delta n$) is larger in the
$P2_{1}/n$ phase (0.23, 0.34 and 0.1 for $n_{1}$-$n_{2}$, $n_{1}$-$n_{3}$ and
$n_{2}$-$n_{3}$, respectively) than in the $P2_{1}/c$ phase [51].

For comparison, we also computed the optical dielectric tensor
${\varepsilon}^{\infty}_{ij}$, and Born effective charges $Z^{*}$ of the
ground state $P2_{1}/c$ phase using the B1-WC and HSE06 hybrid functional. The
results are presented in the appendix Tab. A2, A3. We can conclude that LDA
results are in reasonable agreement with the B1-WC results. However, LDA
usually overestimates the absolute value of ${\varepsilon}^{\infty}$, which is
due to the underestimation of the electronic bandgap in line with the lack of
polarization dependence of this local exchange-correlation functional [52].

##  V Dynamical properties

We focus now on the dynamical properties of the room-temperature $P2_{1}/n$
phase and ground-state $P2_{1}/c$ phase of WO3. In each case, we first discuss
the phonons at the zone center ($\Gamma$ point), reporting both Raman and
infra-red spectra, and also the full phonon dispersion curves. As in the
previous Section, all results in the main text have been obtained in LDA.
Further comparison of LDA and B1-WC frequencies are reported in the appendix
(Table A4), showing very close agreement regarding both phonon frequencies and
eigenvectors.

###  V.1 The room temperature $P2_{1}/n$ phase

####  V.1.1 Irreducible representations at $\Gamma$

For the non-polar $P2_{1}/n$ structure, which belongs to the $C^{5}_{2h}$
point group, there are 32 atoms in the primitive cell and so 96 phonon modes
at each $k$-point. The zone-center phonons can be classified according to the
irreducible representations of this point group as

$\Gamma_{phonons}=24A_{g}\oplus 24A_{u}\oplus 24B_{g}\oplus 24B_{u}.$

These include: (i) three acoustic modes ($\Gamma_{acoustic}=A_{u}+2B_{u}$),
(ii) 48 Raman active modes ($A_{g}$ and $B_{g}$) and (iii) 45 infrared (IR)
active modes ($A_{u}$ and $B_{u}$).

####  V.1.2 Infrared active modes

Table 6:  The TO phonon frequencies of the infrared active modes of the $P2_{1}/n$ room temperature phase within LDA in cm-1. The values in brackets correspond to the experimental measurements reported in [53] Irrep | LDA [Exp] | Irrep | LDA [Exp]
---|---|---|---
$A_{u}$(1) | 46 |  $B_{u}$(1) | 126
$A_{u}$(2) | 57 |  $B_{u}$(2) | 194
$A_{u}$(3) | 74 |  $B_{u}$(3) | 203
$A_{u}$(4) | 108 |  $B_{u}$(4) | 221
$A_{u}$(5) | 134 |  $B_{u}$(5) | 235
$A_{u}$(6) | 206 |  $B_{u}$(6) | 265
$A_{u}$(7) | 226 [230] |  $B_{u}$(7) | 280
$A_{u}$(8) | 257 |  $B_{u}$(8) | 283 [285]
$A_{u}$(9) | 266 |  $B_{u}$(9) | 303
$A_{u}$(10) | 272 |  $B_{u}$(10) | 310 [310]
$A_{u}$(11) | 308 |  $B_{u}$(11) | 315
$A_{u}$(12) | 326 |  $B_{u}$(12) | 332 [335]
$A_{u}$(13) | 333 |  $B_{u}$(13) | 344
$A_{u}$(14) | 346 |  $B_{u}$(14) | 362
$A_{u}$(15) | 362 [370] |  $B_{u}$(15) | 379
$A_{u}$(16) | 425 |  $B_{u}$(16) | 407
$A_{u}$(17) | 446 |  $B_{u}$(17) | 625 [665]
$A_{u}$(18) | 609 |  $B_{u}$(18) | 756
$A_{u}$(19) | 695 |  $B_{u}$(19) | 776
$A_{u}$(20) | 699 |  $B_{u}$(20) | 847 [825]
$A_{u}$(21) | 766 |  $B_{u}$(21) | 1042
$A_{u}$(22) | 1010 [920] |  $B_{u}$(22) | 1048
$A_{u}$(23) | 1015 | - | -

TO modes – The calculated TO frequencies of the IR active phonons are listed
in Tab. 6, besides the measured experimental data (in brackets) from Ref. [53]
and their possible assignments. In fact, the assignment of the IR active modes
remains still experimentally unexplored. Thus, we assigned each experimental
frequency to the mode with largest mode effective charge in the same frequency
range (see Tab. A6 in the appendix). It can be observed that the experimental
values and our calculated frequencies are overall in good agreement.

LO modes –The IR active LO modes take different frequencies when approaching
$\Gamma$ from different directions due to the low symmetry of the $P2_{1}/n$
structure. In Tab. A6 in the appendix, the LO frequencies of the IR active
modes along different directions are reported together with related TO
frequencies and mode effective charges. Since there is no one-to-one
correspondence, the mapping between LO and TO modes is done by comparing the
eigenvectors and looking for the maximum overlap.

The polarities of the $A_{u}$ phonon modes are along the $y$ direction, thus,
they do not show any LO-TO splitting along the $\Gamma-Y$ and $\Gamma-B$
directions. On the other hand, some $A_{u}$ modes experience large shifts of
their LO frequencies along the $[010]$ direction (with a width of 45 cm-1 for
$A_{u}$(8), 48 cm-1 for $A_{u}$(15) and 262 cm-1 for $A_{u}$(19)). Along the
$[011]$ direction, although there is a component parallel to the polarity, no
significant shift of LO frequencies is found for $A_{u}$ modes. The
$A_{u}$(15) and $A_{u}$(19) modes experience a strong LO-TO shifting along
$[110]$ (44 cm-1 and 263 cm-1, respectively). Along the $[111]$ direction the
$A_{u}$(7) mode and again the $A_{u}$(19) mode show giant shifts of 171 cm-1
and 268 cm-1, respectively. In fact, all these giant shifts correspond to
modes with very large mode effective charges that are associated with the
opposite motion of W and O atoms along the $y$ direction.

Table 7:  Phonon frequencies of Raman active modes of the $P2_{1}/n$ room temperature phase within LDA in cm-1. Irrep | LDA | Irrep | LDA
---|---|---|---
$A_{g}$(1) | 48 |  $B_{g}$(1) | 46
$A_{g}$(2) | 53 |  $B_{g}$(2) | 75
$A_{g}$(3) | 59 |  $B_{g}$(3) | 81
$A_{g}$(4) | 72 |  $B_{g}$(4) | 178
$A_{g}$(5) | 79 |  $B_{g}$(5) | 186
$A_{g}$(6) | 117 |  $B_{g}$(6) | 202
$A_{g}$(7) | 159 |  $B_{g}$(7) | 218
$A_{g}$(8) | 206 |  $B_{g}$(8) | 249
$A_{g}$(9) | 214 |  $B_{g}$(9) | 255
$A_{g}$(10) | 269 |  $B_{g}$(10) | 291
$A_{g}$(11) | 271 |  $B_{g}$(11) | 312
$A_{g}$(12) | 304 |  $B_{g}$(12) | 336
$A_{g}$(13) | 321 |  $B_{g}$(13) | 354
$A_{g}$(14) | 329 |  $B_{g}$(14) | 359
$A_{g}$(15) | 342 |  $B_{g}$(15) | 363
$A_{g}$(16) | 419 |  $B_{g}$(16) | 393
$A_{g}$(17) | 434 |  $B_{g}$(17) | 430
$A_{g}$(18) | 443 |  $B_{g}$(18) | 433
$A_{g}$(19) | 589 |  $B_{g}$(19) | 608
$A_{g}$(20) | 689 |  $B_{g}$(20) | 738
$A_{g}$(21) | 704 |  $B_{g}$(21) | 756
$A_{g}$(22) | 766 |  $B_{g}$(22) | 842
$A_{g}$(23) | 809 |  $B_{g}$(23) | 1003
$A_{g}$(24) | 1051 |  $B_{g}$(24) | 1079

As illustrated in Table A6, the $B_{u}$ modes have their polarity in the
$xz$-plane. The two largest LO-TO splittings of 334 and 275 cm-1 are observed
for the $B_{u}$(17) and $B_{u}$(1) modes respectively and are related to the
opposite motion of W and O atoms along $x$. The mode effective charges of the
$B_{u}$(17) and $B_{u}$(1) modes are 18.1$e$ and 19.5$e$, respectively, along
the $x$ direction, and almost zero along $z$, confirming that these modes are
mostly polarized along the $x$ direction. Note that this point is in line with
the considerably large (xx) components of the $Z^{*}$($W_{1}$),
$Z^{*}$($W_{2}$), $Z^{*}$($O_{1}$) and $Z^{*}$($O_{2}$) tensors of the
$P2_{1}/n$ phase.

Figure 6: The calculated infrared reflectivity of the $P2_{1}/n$ phase: (a) $B_{u}$ with the polarity along $x$, (b) $A_{u}$, (c) $B_{u}$ with the polarity along $z$, also the experimental reflectivity spectrum of the $[001]$ crystal surface of WO3 measured at room temperature in the spectral range from 50 to 1200 cm-1, taken from Ref. [54]. Table 8:  The comparison between calculated Raman modes of WO3 $P2_{1}/n$ phase within the LDA and available experimental data. The bold fonts represent the frequencies of Raman lines in the calculated Raman spectrum of Fig. 8. Irrep | LDA | Expt.[55] | Expt.[56] | Expt.[57] | Expt.[53] | Expt.[51]
---|---|---|---|---|---|---
$A_{g}$(1) | 48 | 34 | - | - | 33-$A_{g}$ | -
$A_{g}$(3) | 59 | 61 | - | 63 | 60 | -
$A_{g}$(4) | 72 | 71 | - | 71 | 73 | -
$A_{g}$(5) | 79 | 93 | - | - | 93 | -
$A_{g}$(6) | 117 | 134 | 137 | 136 | 133-$A_{g}$ | -
$A_{g}$(7) | 159 | 187 | 187 | 187 | - | -
$A_{g}$(9) | 214 | 218 | - | 221 | - | -
$A_{g}$(11) | 271 | 273 | 275 | 273 | 275 | 272
$A_{g}$(13) | 321 | 327 | 329 | 328 | 330 | 322
$B_{g}$(16) | 393 | - | - | 376 | - | -
$A_{g}$(17) | 434 | 434 | - | 436 | - | 430
$A_{g}$(18) | 443 | - | - | - | - | 447
$A_{g}$(20) | 689 | 715 | 718 | 717 | 719-$A_{g}$ | 709
$A_{g}$(22) | 766 | 807 | 808 | 808 | 808-$A_{g}$ | 807
Figure 7: The polarized Raman spectra of the $P2_{1}/n$ phase, calculated
within the LDA (for visualization purposes the low frequency range between
100-550 cm-1 is magnified).

IR reflectivity spectra – The infrared reflectivity spectra are also
calculated at normal incidence for all three Cartesian directions, showing
contributions of $B_{u}$ (along $x$), $A_{u}$ (along $y$) and of $B_{u}$
(along $z$) modes respectively. The experimental reflectivity spectrum of the
$[001]$ crystal surface of WO3 measured by Gabrusenoksis et al. [54] is also
presented in panel $(c)$ of Fig. 6. In our calculations, the reflectivity
spectra saturate to unity because our formalism neglects damping. In spite of
that, a reasonably good agreement between the theoretical and experimental
reflectivity spectrum is observed. However it seems that the small feature at
850 cm-1 corresponding to the $B_{u}(20)$ mode is missed in the theoretical
spectrum. This mode is correctly predicted by the calculation but is predicted
with a negligible mode charge.

####  V.1.3 Raman active modes

Raman modes – Tab. 7 shows the calculated frequencies of Raman active phonons
classified according to their symmetry. To the best of our knowledge, the only
theoretical study of the Raman spectrum of the $P2_{1}/n$ phase of WO3 so far
is done by Yang et al. [22]. In contrast to both experimental data and our
calculations, they report three dominant Raman peaks in the high-frequency
range, located around 610, 650, 850 cm-1. The comparison between calculated
Raman active mode frequencies, and available experimental data is given in
Tab. 8 where a reasonable agreement is observed. In Raman measurements of Ref
[53], it has been mentioned that the crystal samples were so small that the
Raman tensors could not be determined except for the ‘c’ element corresponding
to $A_{g}$ symmetry and these modes are in line with the assignments from our
calculations. The peak located at 393 cm-1 has only been experimentally
measured by Garcia-Sanchez et al. [57] at a frequency of 376 cm-1. The reason
that this peak is absent in most of the other experimental studies is that it
corresponds to irreducible representation $B_{g}$ and the Raman polarizability
tensor for $B_{g}$ has only non-zero off-diagonal components (see next
Section), resulting in a very low Raman intensity (the concentration of the
active species is very low in off-diagonal directions) [58, 59, 60, 61, 62].

Raman spectra – Based on group theory, the Raman susceptibility tensors for
the $A_{g}$ and the $B_{g}$ modes are given by: $A_{g}$=
$\begin{pmatrix}a&d&0\\\ d&b&0\\\ 0&0&c\end{pmatrix}$ , $B_{g}$=
$\begin{pmatrix}0&0&e\\\ 0&0&f\\\ e&f&0\end{pmatrix}$ so that tensor elements
could be determined by the polarized Raman scattering measurement on single
crystals in different experimental configurations. As illustrated in Fig. 7,
the calculated polarized Raman spectra show different intensities depending on
the geometrical orientation. Unfortunately, no experimental polarized Raman
spectra have been reported so far for WO3 single crystals so that our
calculations can be used as benchmark results for the interpretation of future
measurements. We also determined the powder spectrum by averaging over
different orientations of the polarized spectrum. The experimental Raman
spectrum measured on nanoparticle powder of WO3 reported in a recent work by
Thummavichai et al. [56], and the calculated powder Raman spectrum are
compared in Fig. 8. It appears there that our calculations provide an overall
good description for the powder Raman spectrum. In the high-frequency range,
the spectrum has two peaks at the correct frequency position and with a
relatively good prediction of the intensities. In the low-frequency range
peaks are at the right frequency position, but the intensities are strongly
underestimated.

Figure 8: The unpolarized Raman spectra of the $P2_{1}/n$ phase, the room
temperature structure of WO3, the calculated within the LDA (for the
visualization purposes the low frequency range (between 100-550 cm-1) are
magnified), and the experimental Raman spectrum measured on nanoparticle
powder of WO3. The spectrum has been acquired using a 532 nm laser at room
temperature [56].

The two intense peaks in the high-frequency range are dominated by antipolar
motions of the oxygen atoms (heavy W atoms are not moving significantly in
this frequency range) along $z$ and $y$ orientations (see III.2). The
$A_{g}$(22) mode at 766 cm-1 with the maximum intensity is dominated by the
similar antipolar mode at M in the reference cubic phase (antipolar
M${}^{-}_{3}$ mode related to oxygen motions along $z$). The second high-
intensity peak at 689 cm-1 is associated to the $A_{g}$(20) mode which is
dominated by the similar antipolar mode at X in the reference cubic phase
(antipolar X${}^{-}_{5}$ mode related to oxygen motions along $y$). Also, in
the low-frequency range, we observe peaks associated to modes with dominant
antipolar motions, but involving this time both W and O atoms ($A_{g}$(6) and
$A_{g}$(11)). Other peaks are mostly associated with W-O-W bending modes of
the bridging oxygen.

####  V.1.4 Phonon dispersion curves

Fig. 9 illustrates the phonon dispersion curves and the phonon density of
states of the $P2_{1}/n$ phase. Red and blue colors distinguish the
involvement of W and O atoms. There are no imaginary frequencies in the phonon
bands, indicating that this phase, although not the ground state, is
dynamically stable. The phonon bands are spread into three regions separated
by two gaps of 100-150 cm-1. The low-frequency region is associated to modes
involving motions of W and O atoms while the medium-frequency and high-
frequency regions concern nearly pure oxygen motions.

Figure 9: The calculated full phonon dispersion curve and phonon density of
states of the room temperature structure $P21/n$ of WO3 within LDA. The blue
and red colors distinguish the contributions of oxygen and tungsten atoms,
respectively.

###  V.2 The ground state $P2_{1}/c$ phase

In this section, we report the dynamical properties of the $P2_{1}/c$ phase
with the LDA functional and further assess the dynamical stability of the
$P2_{1}/c$ phase with respect to a potential $Pc$ ground state (discussed in
Sec. III.1).

####  V.2.1 Irreducible representations at $\Gamma$

The non-polar $P2_{1}/c$ structure is similar to the $P2_{1}/n$ phase and it
belongs to the $C^{5}_{2h}$ point group. However, the primitive unit cell is
rotated by 45∘ and only contains 16 atoms so that there are 48 $\Gamma$-phonon
modes at each $k$-point. The zone-center phonons can be classified according
to the irreducible representations of this point group as

$\Gamma_{phonons}=12A_{g}\oplus 12A_{u}\oplus 12B_{g}\oplus 12B_{u}.$

These include: (i) three acoustic modes ($\Gamma_{acoustic}=A_{u}+2B_{u}$),
(ii) 24 Raman active modes ($A_{g}$ and $B_{g}$) and (iii) 21 infrared (IR)
active modes ($A_{u}$ and $B_{u}$).

####  V.2.2 Infrared active modes

TO modes – In the first part of Tab. 9 we report the calculated TO frequencies
of the infrared active modes. According to this Table, there is no unstable
mode at $\Gamma$: the lowest polar mode has a frequency of 139 cm-1 so that it
is far from being unstable and the result appears robust (in agreement with
Hamdi et al. [14]).

LO modes – The investigation of the LO modes along the different high-symmetry
directions of the Brillouin zone, their corresponding irreducible
representations, and mode effective charges $\overline{Z}^{*}$ of the infrared
active modes along the $x$, $y$ and $z$ direction are presented in Tab. A5 of
the appendix. This Table shows that along all the high-symmetry directions of
the Brillouin zone, there are only three large LO-TO splittings in the same LO
frequency ranges: at 260 $\pm$ 20 cm-1, 415 $\pm$ 10 cm-1, and 960 $\pm$ 20
cm-1. In addition, the strongest LO-TO splitting is located at 962 cm-1 with a
width of 322 cm-1 and is associated to the $B_{u}$(8) mode in the $[101]$
direction. This mode is associated with the opposite motions of W and O along
the $x$ direction.

IR reflectivity spectra – The infrared reflectivity spectra are also
calculated at normal incidence for all the Cartesian directions. Comparing to
the reflectivity spectra of the $P2_{1}/n$ phase, here, between 750 to 800
cm-1 two peaks of $A_{u}$ and $B_{u}(z)$ are absent. This point can be
exploited for distinguishing these two phases in future infrared measurements.

Figure 10: The calculated infrared reflectivity of the $P2_{1}/c$ phase: (a) $B_{u}$ with the polarity along $x$, (b) $A_{u}$, (c) $B_{u}$ with the polarity along $z$. Figure 11: The polarized Raman spectra of the $P2_{1}/c$ phase of WO3 calculated within the LDA (for visualization purposes the low frequency range between 450-150 cm-1 is magnified). Figure 12: The unpolarized Raman spectra of the $P2_{1}/c$ phase, the ground state of WO3, calculated within LDA (for the visualization purposes the low frequency range (between 450-150 cm-1) is magnified). Table 9: The TO phonon frequencies of the ground state $P2_{1}/c$ phase at $\Gamma$ within LDA. Irrep | Freq(cm-1) | Irrep | Freq(cm-1)
---|---|---|---
$A_{u}$ | 139 | $B_{u}$ | 166
$A_{u}$ | 174 | $B_{u}$ | 241
$A_{u}$ | 221 | $B_{u}$ | 272
$A_{u}$ | 255 | $B_{u}$ | 318
$A_{u}$ | 284 | $B_{u}$ | 345
$A_{u}$ | 302 | $B_{u}$ | 360
$A_{u}$ | 339 | $B_{u}$ | 453
$A_{u}$ | 365 | $B_{u}$ | 640
$A_{u}$ | 659 | $B_{u}$ | 766
$A_{u}$ | 737 | $B_{u}$ | 1041
$A_{u}$ | 1020 |  |
$A_{g}$ | 50 | $B_{g}$ | 53
$A_{g}$ | 79 | $B_{g}$ | 76
$A_{g}$ | 136 | $B_{g}$ | 178
$A_{g}$ | 179 | $B_{g}$ | 196
$A_{g}$ | 195 | $B_{g}$ | 272
$A_{g}$ | 252 | $B_{g}$ | 327
$A_{g}$ | 294 | $B_{g}$ | 410
$A_{g}$ | 372 | $B_{g}$ | 423
$A_{g}$ | 427 | $B_{g}$ | 450
$A_{g}$ | 644 | $B_{g}$ | 628
$A_{g}$ | 716 | $B_{g}$ | 839
$A_{g}$ | 765 | $B_{g}$ | 1076
Figure 13: The calculated full phonon dispersion curves and phonon density of
states of the ground state $P2_{1}/c$ phase of WO3 within the LDA. The blue
and red colors represent the contribution of oxygen and tungsten atoms,
respectively.

####  V.2.3 Raman active modes

Raman modes – In the second part of Tab. 9 we report the calculated TO
frequencies of the Raman active modes classified according to their
irreducible representations. We are not aware of any experimental data nor
other theoretical values for comparison.

Raman spectra – The calculated polarized, and unpolarized Raman spectra of the
$P2_{1}/c$ phase are presented in Fig. 11 and Fig. 12, respectively. The Raman
susceptibility tensors for the $A_{g}$ and $B_{g}$ modes are the same as in
the $P2_{1}/n$ phase. Although the number of Raman active modes of the
$P2_{1}/c$ phase are half the number of Raman active modes of the $P2_{1}/n$
phase, more peaks are observed in the Raman spectrum of the $P2_{1}/c$ phase,
because the Raman intensities are higher in the $P2_{1}/c$ phase. This might
be because in this phase the density of identical bond characters are more
pronounced. In the high-frequency range, the powder Raman spectrum of the
$P2_{1}/n$ phase shows only two intense peaks, whereas there are three intense
and one small peak in $P2_{1}/c$ phase spectrum. The small peak at 716 cm-1
corresponds mainly to the Jahn–Teller distortion. The mode with the highest
intensity at 766 cm-1 is dominated by the antipolar oxygen motions along the
$z$ direction (related to M${}^{-}_{3}$ from the reference cubic phase, see
III.2). Other two peaks with noticeable intensities are also dominated by the
antipolar oxygen motions but this time in the $xy$-plane (X${}^{-}_{5}$
antipolar mode from the reference cubic phase).

As a result, the observation of the high-frequency range of the spectrum can
provide a clear way to distinguish the $P2_{1}/n$ and $P2_{1}/c$ phases. More
generally, it is worth noticing that the number of intense peaks in the high-
frequency range reflects the number of antipolar distortions that have been
condensed in the crystal structure. Indeed, condensing an anti-polar
distortion transforms the related anti-polar mode originally at the zone
boundary into a Raman active mode at $\Gamma$. In the $P2_{1}/n$ (resp.
$P2_{1}/c$) phase, we have antipolar distortions along $z$ and $y$ directions
(resp. $z$, $y$ and $x$ directions) and two (resp. three) intense peaks. In
order to further validate this statement, we calculated also the Raman
spectrum of $P4/ncc$ and $P4/nmm$ phases, both containing only one antipolar
distortion along $z$ (M${}^{-}_{3}$) and, in both cases, only one intense peak
appears in the high-frequency range.

Comparing the Raman spectrum of the $P2_{1}/n$ and $P2_{1}/c$ phases, we
further observe a small shift of the peaks. The larger the alternation of
long-short $O-W-O$ bonds, the higher is the frequency of the peaks. The 766
cm-1 peak is an exact match. This is expected because $P2_{1}/c$ and
$P2_{1}/n$ phases have the same long-short $O-W-O$ bonds length equal to
1.76-2.09 Å along $z$. The 689 cm-1 peak in the Raman spectrum of the
$P2_{1}/n$ is attributed to the long-short $O-W-O$ bonds equal to 1.78-2.04 Å
along $y$, while in the $P2_{1}/c$ phase, the second peak is shifted to the
lower frequency (644 cm-1) due to the smaller alternation of long-short
$O-W-O$ bonds (1.81-1.95 Å along $y$).

####  V.2.4 Phonon dispersion curves

The phonon dispersion curves and phonon projected density of states of the
$P2_{1}/c$ phase are shown in Fig. 13. Similarly to the $P2_{1}/n$ phase,
there is no imaginary frequency, confirming the dynamical stability of this
phase. Again, the frequencies are spread into 3 well-separated regions. The
main differences between the two phases appear in the intermediate region,
involving the stretching vibration of O-W-O. Also the DOS of the low-energy
region is shifted to slightly higher frequencies in the $P2_{1}/c$ phase.

##  VI Conclusions

In this work we have verified the ability of various DFT functionals (B1-WC,
HSE06, LDA, GGA-WC, PBEsol) to describe the structural energetics of all the
experimentally observed phases as well as six additional metastable phases of
WO3. We have analyzed the atomic distortions of these phases in terms of
symmetry adapted modes of the hypothetical cubic phase and quantified the
amplitude of the distortions. We have analyzed their energetics based on the
instabilities presented in the phonon dispersion curve of the parent cubic
phase. It can be concluded that (i) the B1-WC hybrid functional is the most
appropriate to describe WO3 while (ii) among the investigated standard
functionals, LDA is the most suitable to reproduce the structural distortions
of the different phases of WO3 and to predict the correct ground state. We
have also re-discussed briefly the unusual phase diagram of WO3.

We have exhaustively described the dielectric and lattice dynamical properties
of the ground state $P2_{1}c$ phase and room temperature $P2_{1}/n$ phase. The
Born effective charges, dielectric tensors and the refractive indices were
obtained. The phonon dispersion curves were also reported; no imaginary
frequency was observed, indicating that both phases are dynamically stable.
The $\Gamma$-phonon modes of both phases were carefully studied and compared
to available data. The infrared reflectivity and Raman spectra were also
calculated and reasonable agreement between experimental and calculated
spectra was found for the room temperature phase. Also, we have carefully
discussed the Raman spectra, explaining the physical origin of their main
features and evolution from one phase to another. We pointed out that the
alternation of long-short bonds along $O-W-O$ chains can lead to important
consequences on the dynamical properties of different phases of WO3. We also
revealed that the number of peaks in the high-frequency range of the Raman
spectrum appears as a fingerprint of the number of antipolar distortions that
are present in the structure so that measurement of the Raman spectrum of WO3
appears as an efficient way to distinguish between its different phases.

Our calculations appear as benchmark results for the interpretation of future
experimental measurements. Furthermore, many physical properties of WO3 are
affected by polarons, which arise from the coupling of excess charge carriers
with specific phonon modes. Our work provides all the ingredients for future
more systematic analysis of which phonon modes are involved in polaron
formation and key to localize the charge.

## Acknowledgements

This work has been funded by the Communauté Française de Belgique (ARC AIMED
G.A. 15/19-09) and a Methusalem project of the University of Antwerp . EB
thanks the FRS-FNRS for support. The authors acknowledge the CECI
supercomputer facilities funded by the F.R.S-FNRS (Grant No. 2.5020.1), the
Tier-1 supercomputer of the Fédération Wallonie-Bruxelles funded by the
Walloon Region (Grant No. 1117545), and the computing facilities of the
Flemish Supercomputer Center. We acknowledge that the results of this research
have been achieved using the DECI resource BEM based in Poland at Wrocław with
support from the PRACE OFFSPRING project.

##  Appendix A Phonon dispersion calculations

In the CRYSTAL17 package, the phonon frequencies are computed using the frozen
phonon numerical differences approach. Fig. 1(a) shows the phonon bands of the
parent cubic phase of WO3 generated by CRYSTAL17. The unphysical oscillations
and intersections observed in the high frequency range of the phonon bands are
due to the improper implementation of the mixed-space approach [63] (the
cancellation of the artificial macroscopic electric field is not applied)
which causes problems in Fourier interpolations at arbitrary, non high
symmetry q-points. As a result, the correct phonon bands has been obtained
just after a proper use of the Fourier interpolation technique (shown in Fig.
1(b)).

(a)

(b)

Figure A1: The phonon dispersion curves of the cubic phase of WO3, obtained by
(a) the CRYSTAL17 code, and the improper implementation of the mixed-space
approach, and (b) the ABINIT code, and the proper employment of Fourier
interpolation technique.

(a)

(b)

(c)

(d)

Figure A2: The comparison of phonon dispersion curves of the cubic phase obtained by B1-WC and (a) HSE06, (b) LDA, (c) PBEsol, (d) GGA-WC. Table A1: Atomic positions of $P2_{1}/n$ and $P2_{1}/c$ phases of $WO_{3}$ in reduced coordinates. Only atomic positions are relaxed by LDA with the fixed lattice constants of the experimental values, and calculated within B1-WC for $P2_{1}/n$ (7.30, 7.53, 7.69 Å), and $P2_{1}/c$ (5.26, 5.15, 7.61 Å), respectively. Phase | atom | Wyckoff | Position
---|---|---|---
$P2_{1}/n$ | $W_{1}$ | 4e | 0.2487 0.0274 0.2203
| $W_{2}$ | 4e | 0.2507 0.0286 0.7211
| $O_{1}$ | 4e | 0.5011 0.0338 0.2713
| $O_{2}$ | 4e | 0.4993 0.4657 0.2707
| $O_{3}$ | 4e | 0.2158 0.2611 0.2380
| $O_{4}$ | 4e | 0.2851 0.2611 0.7509
| $O_{5}$ | 4e | 0.2269 0.0199 0.4927
| $O_{6}$ | 4e | 0.2281 0.5067 0.5067
$P2_{1}/c$ | $W_{1}$ | 4e | 0.2541 -0.2338 -0.2180
| $O_{1}$ | 4e | 0.0343 0.0383 -0.2209
| $O_{2}$ | 4e | 0.4580 0.4637 -0.2794
| $O_{3}$ | 4e | 0.2521 -0.3077 0.0078
Table A2:  The Born effective charge tensors ($Z^{*}$) of the $P2_{1}/c$ phase of WO3, calculated by both LDA and B1-WC approaches (see notations of Tab. A1).  B1-WC | LDA
---|---
$Z^{*}$($W_{1}$) |  $Z^{*}$($W_{1}$)
$\begin{pmatrix}9.97&-1.07&-0.16\\\ 0.9&8.78&1.05\\\ -0.19&-1.09&7.44\end{pmatrix}$ | $\begin{pmatrix}10.38&-1.07&-0.15\\\ 0.99&9.40&0.95\\\ -0.17&-0.89&7.74\end{pmatrix}$
$Z^{*}$($O_{1}$) |  $Z^{*}$($O_{1}$)
$\begin{pmatrix}-4.4&3.04&-0.63\\\ 2.75&-3.87&0.47\\\ -0.28&0.26&1\end{pmatrix}$ | $\begin{pmatrix}-4.57&3.05&-0.66\\\ 2.92&-4.13&0.54\\\ -0.24&0.23&-1.06\end{pmatrix}$
$Z^{*}$($O_{2}$) |  $Z^{*}$($O_{2}$)
$\begin{pmatrix}-4.52&2.83&0.6\\\ 2.71&-3.73&-0.54\\\ 0.29&0.3&-1.04\end{pmatrix}$ | $\begin{pmatrix}-4.62&2.90&0.65\\\ 2.89&-4.01&-0.57\\\ 0.23&-0.23&-1.11\end{pmatrix}$
$Z^{*}$($O_{3}$) |  $Z^{*}$($O_{3}$)
$\begin{pmatrix}-1.05&-0.006&0.19\\\ -0.04&-1.2&0.34\\\ 0.17&0.6&-5.42\end{pmatrix}$ | $\begin{pmatrix}-1.19&0.00&0.16\\\ -0.03&-1.25&0.21\\\ 0.18&0.45&-5.56\end{pmatrix}$
Table A3: Comparison between optical dielectric tensor (${\varepsilon}^{\infty}_{ij}$) of the $P2_{1}/c$ phase of WO3, calculated by HSE06, B1-WC hybrid functional and LDA.  HSE06 | B1-WC | LDA
---|---|---
$\begin{pmatrix}5.34&0&-0.09\\\ 0&4.94&0\\\ -0.09&0&4.40\end{pmatrix}$ | $\begin{pmatrix}6.06&0&-0.09\\\ 0&5.53&0\\\ -0.09&0&4.87\end{pmatrix}$ | $\begin{pmatrix}6.95&0&-0.08\\\ 0&6.50&0\\\ -0.08&0&5.44\end{pmatrix}$

##  Appendix B Comparison of LDA and B1-WC for lattice dynamical properties of
$P2_{1}/c$ phase

In order to be able to accurately match the $\Gamma$\- phonon modes obtained
by LDA and B1-WC, we have evaluated the overlap matrix between the phonon
eigenmodes computed from the diagonalization of the dynamical matrix by
ABINIT/LDA (${\eta}_{A}$) and those computed by CRYSTAL/B1-WC (${\eta}_{C}$),
namely
$\left\langle{\eta}_{C}\mathrel{\left|\vphantom{{\eta}_{C}M{\eta}_{A}}\right.\kern-1.2pt}M\mathrel{\left|\vphantom{{\eta}_{C}M{\eta}_{A}}\right.\kern-1.2pt}{\eta}_{A}\right\rangle$,
where $M$ is the mass matrix. In case of two corresponding eigenmodes, this
overlap
$\left\langle{\eta}_{C}\mathrel{\left|\vphantom{{\eta}_{C}M{\eta}_{A}}\right.\kern-1.2pt}M\mathrel{\left|\vphantom{{\eta}_{C}M{\eta}_{A}}\right.\kern-1.2pt}{\eta}_{A}\right\rangle$
should be approximately 1. We could realize such a one-to-one correspondence
between the phonon eigenmodes, and thus also a match between the corresponding
phonon frequencies, which are shown in Tab. A4. The comparison of the
calculated phonon frequencies by both LDA and B1-WC indicates that they are in
good agreement.

Table A4: The TO phonon frequencies of the ground state $P2_{1}/c$ phase at $\Gamma$ within the B1-WC and LDA approaches, in cm-1, and the overlap matrix between the phonon eigenmodes computed from the diagonalization of the dynamical matrix by ABINIT/LDA (${\eta}_{A}$) and those computed by CRYSTAL/B1-WC (${\eta}_{C}$).  Irrep | B1-WC | LDA | $\left\langle{\eta}_{C}\mathrel{\left|\vphantom{{\eta}_{C}M{\eta}_{A}}\right.\kern-1.2pt}M\mathrel{\left|\vphantom{{\eta}_{C}M{\eta}_{A}}\right.\kern-1.2pt}{\eta}_{A}\right\rangle$
---|---|---|---
$A_{g}$ | 59 | 50 | 1.00
$A_{g}$ | 93 | 79 | 1.00
$A_{g}$ | 143 | 136 | 0.97
$A_{g}$ | 183 | 179 | 0.92
$A_{g}$ | 215 | 195 | 0.93
$A_{g}$ | 269 | 252 | 0.99
$A_{g}$ | 310 | 294 | 0.98
$A_{g}$ | 386 | 372 | 0.98
$A_{g}$ | 429 | 427 | 1.00
$A_{g}$ | 681 | 644 | 1.00
$A_{g}$ | 750 | 716 | 1.00
$A_{g}$ | 805 | 765 | 1.00
$B_{g}$ | 64 | 53 | 0.98
$B_{g}$ | 91 | 76 | 0.98
$B_{g}$ | 204 | 178 | 0.92
$B_{g}$ | 174 | 196 | 0.91
$B_{g}$ | 281 | 272 | 0.96
$B_{g}$ | 332 | 327 | 0.97
$B_{g}$ | 429 | 410 | 0.96
$B_{g}$ | 418 | 423 | 0.99
$B_{g}$ | 457 | 450 | 0.96
$B_{g}$ | 640 | 628 | 1.00
$B_{g}$ | 866 | 839 | 1.00
$B_{g}$ | 1083 | 1076 | 1.00
$A_{u}$ | 144 | 139 | 1.00
$A_{u}$ | 184 | 174 | 0.98
$A_{u}$ | 216 | 221 | 0.98
$A_{u}$ | 262 | 255 | 0.96
$A_{u}$ | 287 | 284 | 0.95
$A_{u}$ | 308 | 302 | 0.96
$A_{u}$ | 359 | 339 | 0.98
$A_{u}$ | 383 | 365 | 0.98
$A_{u}$ | 698 | 659 | 1.00
$A_{u}$ | 765 | 737 | 1.00
$A_{u}$ | 1038 | 1020 | 1.00
$B_{u}$ | 188 | 166 | 0.99
$B_{u}$ | 234 | 241 | 1.00
$B_{u}$ | 276 | 272 | 1.00
$B_{u}$ | 327 | 318 | 0.99
$B_{u}$ | 350 | 345 | 0.99
$B_{u}$ | 383 | 360 | 0.99
$B_{u}$ | 456 | 453 | 1.00
$B_{u}$ | 651 | 640 | 1.00
$B_{u}$ | 796 | 766 | 1.00
$B_{u}$ | 1052 | 1041 | 1.00

##  Appendix C LO-TO modes

In Tab. A5 and Tab. A6, the LO frequencies of the IR active modes for the
$P2_{1}/c$ and $P2_{1}/n$ phase, respectively, are reported along different
directions together with related TO frequencies and mode effective charges.

Table A5:  The LO frequencies (in cm-1) of the infrared active modes along the different high-symmetry directions of the Brillouin zone, their corresponding irreducible representations, and mode effective charges ($\overline{Z}^{*}$) of the TO infrared active modes along the $x$, $y$ and $z$ is presented for the $P2_{1}/c$ phase calculated by LDA. The frequencies shown by light gray in the LO columns are equal to the TO frequencies. The frequencies shown in bold font present the strong shifts in LO-TO splitting. The coordinates of the high-symmetry points are as follows: Y(1/2,0,0), Z(0,1/2,0), B(0,0,1/2), C(1/2,1/2,0), D(0,1/2,1/2), A(1/2,0,1/2) and E(1/2,1/2,1/2). Irrep | TO Modes | $\overline{Z}^{*}$ | LO Modes
---|---|---|---
|  | x | y | z | Y | Z | B | C | D | A | E
$A_{u}$(1) | 139 | 0 | 2.8 | 0 | 139 | 143 | 139 | 141 | 143 | 139 | 141
$A_{u}$(2) | 174 | 0 | 9.9 | 0 | 174 | 239 | 174 | 171 | 226 | 174 | 171
$A_{u}$(3) | 221 | 0 | -2.2 | 0 | 221 | 217 | 221 | 221 | 211 | 221 | 221
$A_{u}$(4) | 255 | 0 | -2.6 | 0 | 255 | 267 | 255 | 252 | 261 | 255 | 245
$A_{u}$(5) | 284 | 0 | -4.6 | 0 | 284 | 412 | 284 | 274 | 306 | 284 | 303
$A_{u}$(6) | 302 | 0 | -0.3 | 0 | 302 | 302 | 302 | 302 | 302 | 302 | 302
$A_{u}$(7) | 339 | 0 | -3.6 | 0 | 339 | 317 | 339 | 329 | 347 | 339 | 342
$A_{u}$(8) | 365 | 0 | -2.3 | 0 | 365 | 358 | 365 | 365 | 365 | 365 | 365
$A_{u}$(9) | 659 | 0 | 11.8 | 0 | 659 | 941 | 659 | 650 | 962 | 659 | 956
$A_{u}$(10) | 737 | 0 | -2.1 | 0 | 737 | 734 | 737 | 737 | 737 | 737 | 725
$A_{u}$(11) | 1020 | 0 | 0.4 | 0 | 1020 | 1020 | 1020 | 1020 | 1020 | 1020 | 1020
$B_{u}$(1) | 166 | 12.1 | 0 | 0.0 | 413 | 166 | 166 | 412 | 166 | 289 | 262
$B_{u}$(2) | 241 | 3.9 | 0 | -2.9 | 232 | 241 | 248 | 234 | 246 | 241 | 241
$B_{u}$(3) | 272 | -0.1 | 0 | 3.1 | 272 | 272 | 290 | 272 | 277 | 251 | 277
$B_{u}$(4) | 318 | 2.9 | 0 | -1.5 | 301 | 318 | 321 | 308 | 321 | 318 | 318
$B_{u}$(5) | 345 | -0.8 | 0 | -1.0 | 341 | 345 | 345 | 345 | 345 | 345 | 345
$B_{u}$(6) | 360 | 3.3 | 0 | 5.1 | 349 | 360 | 410 | 353 | 411 | 426 | 421
$B_{u}$(7) | 453 | -0.4 | 0 | 0.3 | 453 | 453 | 453 | 453 | 453 | 453 | 453
$B_{u}$(8) | 640 | 12.5 | 0 | -0.1 | 952 | 640 | 640 | 947 | 640 | 962 | 649
$B_{u}$(9) | 766 | 0.3 | 0 | -11.1 | 766 | 766 | 978 | 766 | 709 | 706 | 739
$B_{u}$(10) | 1041 | -0.2 | 0 | -0.5 | 1041 | 1041 | 1043 | 1041 | 1041 | 1043 | 1041
Table A6:  The LO frequencies (in cm-1) of the infrared active modes along the different high-symmetry directions of the Brillouin zone, their corresponding irreducible representations, and mode effective charges ($\overline{Z}^{*}$) of the TO infrared active modes along the $x$, $y$ and $z$ is presented for the $P2_{1}/n$ phase calculated by LDA. The frequencies shown by light gray in the LO columns are equal to the TO frequencies. The frequencies shown in bold font present the strong shifts in LO-TO splitting. The coordinates of the high-symmetry points are as follows: Y(1/2,0,0), Z(0,1/2,0), B(0,0,1/2), C(1/2,1/2,0), D(0,1/2,1/2), A(1/2,0,1/2) and E(1/2,1/2,1/2). Irrep | TO Modes | $\overline{Z}^{*}$ | LO Modes
---|---|---|---
|  | x | y | z | Y | Z | B | C | D | A | E
$A_{u}$(1) | 46 | 0 | 0.2 | 0 | 46 | 46 | 46 | 46 | 46 | 46 | 46
$A_{u}$(2) | 57 | 0 | -0.3 | 0 | 57 | 57 | 57 | 57 | 57 | 57 | 57
$A_{u}$(3) | 74 | 0 | -0.6 | 0 | 74 | 74 | 74 | 74 | 74 | 74 | 74
$A_{u}$(4) | 108 | 0 | -0.04 | 0 | 108 | 108 | 108 | 108 | 108 | 108 | 108
$A_{u}$(5) | 134 | 0 | -0.09 | 0 | 134 | 134 | 134 | 134 | 134 | 134 | 134
$A_{u}$(6) | 206 | 0 | -0.3 | 0 | 206 | 206 | 206 | 206 | 206 | 206 | 212
$A_{u}$(7) | 226 | 0 | -6.9 | 0 | 226 | 236 | 226 | 215 | 233 | 226 | 397
$A_{u}$(8) | 257 | 0 | 7.2 | 0 | 257 | 302 | 257 | 246 | 264 | 257 | 242
$A_{u}$(9) | 266 | 0 | 1.3 | 0 | 266 | 266 | 266 | 266 | 269 |  266 | 266
$A_{u}$(10) | 272 | 0 | -2.8 | 0 | 272 | 272 | 272 | 272 | 275 | 272 | 269
$A_{u}$(11) | 308 | 0 | 1.3 | 0 | 308 | 314 | 308 | 308 | 308 | 308 | 308
$A_{u}$(12) | 326 | 0 | 0.4 | 0 | 326 | 326 | 326 | 326 | 326 | 326 | 326
$A_{u}$(13) | 333 | 0 | 1.4 | 0 | 333 | 335 | 333 | 329 | 333 |  333 | 333
$A_{u}$(14) | 346 | 0 | -2.2 | 0 | 346 | 349 | 346 | 350 | 351 | 346 | 343
$A_{u}$(15) | 362 | 0 | 7.3 | 0 | 362 | 410 | 362 | 407 | 397 | 362 | 352
$A_{u}$(16) | 425 | 0 | 0.03 | 0 | 425 | 425 | 425 | 425 | 425 | 425 | 425
$A_{u}$(17) | 446 | 0 | 0.06 | 0 | 446 | 446 | 446 | 446 | 446 | 446 | 446
$A_{u}$(18) | 609 | 0 | 0.3 | 0 | 609 | 609 | 609 | 609 | 609 | 609 | 609
$A_{u}$(19) | 695 | 0 | -13.9 | 0 | 695 | 957 | 695 | 958 | 698 | 695 | 963
$A_{u}$(20) | 699 | 0 | -7.0 | 0 | 699 | 699 | 699 | 699 | 722 | 699 | 699
$A_{u}$(21) | 766 | 0 | 6.3 | 0 | 766 | 754 | 766 | 756 | 766 | 766 | 766
$A_{u}$(22) | 1010 | 0 | 0.8 | 0 | 1010 | 1013 | 1010 | 1012 | 1012 | 1010 | 1010
$A_{u}$(23) | 1015 | 0 | -0.01 | 0 | 1015 | 1015 | 1015 | 1015 | 1015 | 1015 | 1015
$B_{u}$(1) | 126 | 19.5 | 0 | -0.06 | 401 | 126 | 126 | 337 | 126 | 390 | 183
$B_{u}$(2) | 194 | 0.1 | 0 | -0.4 | 194 | 194 | 194 | 194 | 194 | 194 | 194
$B_{u}$(3) | 203 | 8.3 | 0 | -0.5 | 193 | 203 | 203 | 186 | 203 | 191 | 206
$B_{u}$(4) | 221 | -2.0 | 0 | -1.2 | 221 | 221 | 221 | 221 | 221 | 218 | 221
$B_{u}$(5) | 235 | 5.5 | 0 | -1.1 | 231 |  235 | 235 | 232 | 235 | 231 | 233
$B_{u}$(6) | 265 | -0.7 | 0 | 1.4 | 265 | 265 | 265 | 265 | 265 | 265 | 265
$B_{u}$(7) | 280 | 2.8 | 0 | -2.3 | 275 | 280 | 280 | 278 | 280 | 280 | 280
$B_{u}$(8) | 283 | -1.8 | 0 | -4.1 | 283 | 283 | 291 | 283 | 293 | 257 | 275
$B_{u}$(9) | 303 | 0.01 | 0 | 0.04 | 303 | 303 | 303 | 303 | 303 | 297 | 303
$B_{u}$(10) | 310 | 0.6 | 0 | 5.0 | 310 | 310 | 314 | 310 | 314 | 303 | 300
$B_{u}$(11) | 315 | 1.8 | 0 | -1.1 | 315 | 315 | 322 | 315 | 320 | 315 | 315
$B_{u}$(12) | 332 | 0.3 | 0 | 4.4 | 332 | 332 | 384 | 332 | 341 | 324 | 324
$B_{u}$(13) | 344 | -0.7 | 0 | 0.8 | 344 | 344 | 344 | 344 | 344 | 344 | 344
$B_{u}$(14) | 362 | -1.5 | 0 | -0.3 | 362 | 362 | 362 | 362 | 362 | 362 | 362
$B_{u}$(15) | 379 | 0.2 | 0 | 0.6 | 379 | 379 | 376 | 379 | 379 | 379 | 379
$B_{u}$(16) | 407 | -0.04 | 0 | 1.8 | 407 | 407 | 413 |  407 | 413 | 411 | 411
$B_{u}$(17) | 625 | 18.1 | 0 | 0.6 | 959 | 625 | 625 | 664 | 625 | 700 | 660
$B_{u}$(18) | 756 | 0.6 | 0 | -14.6 | 756 | 756 | 975 | 758 | 967 | 966 | 729
$B_{u}$(19) | 776 | 0.3 | 0 | -6.5 |  776 | 776 | 773 | 776 | 776 | 773 | 776
$B_{u}$(20) | 847 | 0.5 | 0 | 0.4 | 847 | 847 | 847 | 847 | 847 | 847 | 847
$B_{u}$(21) | 1042 | 0.02 | 0 | 0.06 | 1042 | 1042 | 1042 | 1042 | 1042 | 1042 | 1042
$B_{u}$(22) | 1048 | -0.03 | 0 | -0.9 | 1048 |  1048 | 1050 | 1048 | 1050 | 1050 | 1050

## References

  * Niklasson and Granqvist [2007] G. A. Niklasson and C. G. Granqvist, Journal of Materials Chemistry 17, 127 (2007).
  * Deb [2008] S. K. Deb, Solar Energy Materials and Solar Cells 92, 245 (2008).
  * Howard  _et al._ [2001] C. J. Howard, V. Luca, and K. S. Knight, Journal of Physics: Condensed Matter 14, 377 (2001).
  * Locherer  _et al._ [1999] K. R. Locherer, I. P. Swainson, and E. K. H. Salje, Journal of Physics: Condensed Matter 11, 6737 (1999).
  * Kehl  _et al._ [1952] W. L. Kehl, R. G. Hay, and D. Wahl, Journal of Applied Physics 23, 212 (1952).
  * Vogt  _et al._ [1999] T. Vogt, P. M. Woodward, and B. A. Hunter, Journal of Solid State Chemistry 144, 209 (1999).
  * Salje [1977] E. Salje, Acta Crystallographica Section B 33, 574 (1977).
  * Loopstra and Rietveld [1969] B. O. Loopstra and H. M. Rietveld, Acta Crystallographica Section B 25, 1420 (1969).
  * Diehl  _et al._ [1978a] R. Diehl, G. Brandt, and E. Salje, Acta Crystallographica Section B 34, 1105 (1978a).
  * Woodward  _et al._ [1995] P. M. Woodward, A. W. Sleight, and T. Vogt, Journal of Physics and Chemistry of Solids 56, 1305 (1995).
  * Salje  _et al._ [1999] E. Salje, S. Rehmann, F. Pobell, D. Morris, K. Knight, T. Herrmannsdörfer, and M. Dove, Journal of Physics: Condensed Matter 9, 6563 (1999).
  * WOODWARD [1997] WOODWARD, Journal of Solid State Chemistry 131, 9 (1997).
  * Salje [1976] E. Salje, Ferroelectrics 12, 215 (1976).
  * Hamdi  _et al._ [2016] H. Hamdi, E. K. Salje, P. Ghosez, and E. Bousquet, Physical Review B 94, 245124 (2016).
  * Migas  _et al._ [2010] D. B. Migas, V. L. Shaposhnikov, V. N. Rodin, and V. E. Borisenko, Journal of Applied Physics 108, 093713 (2010).
  * Bousquet  _et al._ [2020] E. Bousquet, H. Hamdi, P. Aguado-Puente, E. K. H. Salje, E. Artacho, and P. Ghosez, Physical Review Research 2, 10.1103/physrevresearch.2.012052 (2020).
  * Johansson  _et al._ [2013] M. B. Johansson, G. Baldissera, I. Valyukh, C. Persson, H. Arwin, G. A. Niklasson, and L. Österlund, Journal of Physics Condensed Matter 25, 10.1088/0953-8984/25/20/205502 (2013).
  * González-Borrero  _et al._ [2010] P. P. González-Borrero, F. Sato, A. N. Medina, M. L. Baesso, A. C. Bento, G. Baldissera, C. Persson, G. A. Niklasson, C. G. Granqvist, and A. Ferreira Da Silva, Applied Physics Letters 96, 061909 (2010).
  * Ping  _et al._ [2013] Y. Ping, D. Rocca, and G. Galli, Physical Review B - Condensed Matter and Materials Physics 87, 1 (2013).
  * Liu and Fan [2018] X. Liu and H.-Q. Fan, Royal Society Open Science 5, 171921 (2018).
  * Mascello  _et al._ [2020] N. Mascello, N. A. Spaldin, A. Narayan, and Q. N. Meier, Phys. Rev. Research 2, 033460 (2020).
  * Yang  _et al._ [2019] H. Yang, H. Sun, Q. Li, P. Li, K. Song, B. Song, and L. Wang, Vacuum 164, 411 (2019).
  * Bilc  _et al._ [2008a] D. I. Bilc, R. Orlando, R. Shaltaf, G. M. Rignanese, J. Íñiguez, and P. Ghosez, Physical Review B - Condensed Matter and Materials Physics 77, 165107 (2008a).
  * Causà  _et al._ [1987] M. Causà, R. Dovesi, C. Pisani, R. Colle, and A. Fortunelli, Physical Review B 36, 891 (1987).
  * Gonze and et.al [2009] X. Gonze and P. G. et.al, Computer Physics Communications 180, 2582 (2009).
  * Gonze  _et al._ [2002] X. Gonze, J. M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G. M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, P. Ghosez, J. Y. Raty, and D. C. Allan, Computational Materials Science 25, 478 (2002).
  * Perdew and Wang [1992] J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100, 10.1103/physrevlett.100.136406 (2008).
  * Wu and Cohen [2006] Z. Wu and R. E. Cohen, Phys. Rev. B 73, 10.1103/physrevb.73.235116 (2006).
  * van Setten  _et al._ [2018] M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G.-M. Rignanese, Comput. Phys. Commun. 226, 39 (2018).
  * BROYDEN [1970] C. G. BROYDEN, IMA Journal of Applied Mathematics 6, 222 (1970).
  * Gonze  _et al._ [1992] X. Gonze, D. C. Allan, and M. P. Teter, Phys. Rev. Lett. 68, 3603 (1992).
  * Baroni  _et al._ [1987] S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 58, 1861 (1987).
  * Bilc  _et al._ [2008b] D. I. Bilc, R. Orlando, R. Shaltaf, G.-M. Rignanese, J. Íñiguez, and P. Ghosez, Phys. Rev. B 77, 165107 (2008b).
  * Heyd  _et al._ [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 118, 8207 (2003).
  * Figgen  _et al._ [2009] D. Figgen, K. A. Peterson, M. Dolg, and H. Stoll, Journal of Chemical Physics 130, 10.1063/1.3119665 (2009).
  * Pascale  _et al._ [2004] F. Pascale, C. M. Zicovich-Wilson, F. López Gejo, B. Civalleri, R. Orlando, and R. Dovesi, Journal of Computational Chemistry 25, 888 (2004).
  * Zicovich-Wilson  _et al._ [2004] C. M. Zicovich-Wilson, F. Pascale, C. Roetti, V. R. Saunders, R. Orlando, and R. Dovesi, Journal of Computational Chemistry 25, 1873 (2004).
  * King-Smith and Vanderbilt [1993] R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993).
  * Ferrero  _et al._ [2009] M. Ferrero, B. Civalleri, M. Rérat, R. Orlando, and R. Dovesi, The Journal of Chemical Physics 131, 214704 (2009).
  * Ferrero  _et al._ [2008] M. Ferrero, M. Rérat, B. Kirtman, and R. Dovesi, The Journal of Chemical Physics 129, 244110 (2008).
  * Diehl  _et al._ [1978b] R. Diehl, G. Brandt, and E. Salje, Acta Crystallographica Section B 34, 1105 (1978b).
  * Orobengoa  _et al._ [2009] D. Orobengoa, C. Capillas, M. I. Aroyo, and J. M. Perez-Mato, Journal of Applied Crystallography 42, 820 (2009).
  * Bousquet  _et al._ [2008] E. Bousquet, M. Dawber, N. Stucki, C. Lichtensteiger, P. Hermet, S. Gariglio, J.-M. Triscone, and P. Ghosez, Nature 452, 732 (2008).
  * Young  _et al._ [2015] J. Young, A. Stroppa, S. Picozzi, and J. M. Rondinelli, Journal of Physics Condensed Matter 27, 283202 (2015).
  * Benedek  _et al._ [2015] N. A. Benedek, J. M. Rondinelli, H. Djani, P. Ghosez, and P. Lightfoot, Dalton Trans. 44, 10543 (2015).
  * Detraux  _et al._ [1997] F. Detraux, P. Ghosez, and X. Gonze, Phys. Rev. B 56, 983 (1997).
  * Ghosez  _et al._ [1995] P. Ghosez, X. Gonze, P. Lambin, and J.-P. Michenaud, Phys. Rev. B 51, 6765 (1995).
  * Ghosez  _et al._ [1998] P. Ghosez, J.-P. Michenaud, and X. Gonze, Phys. Rev. B 58, 6224 (1998).
  * Sawada and Danielson [1959] S. Sawada and G. C. Danielson, Physical Review 113, 1008 (1959).
  * Santato  _et al._ [2001] C. Santato, M. Odziemkowski, M. Ulmann, and J. Augustynski, Journal of the American Chemical Society 123, 10639 (2001).
  * Gonze  _et al._ [1995] X. Gonze, P. Ghosez, and R. W. Godby, Phys. Rev. Lett. 74, 4035 (1995).
  * Salje [1975] E. Salje, Acta Crystallographica Section A 31, 360 (1975).
  * Gabrusenoks  _et al._ [2001] J. Gabrusenoks, A. Veispals, A. von Czarnowski, and K.-H. Meiwes-Broer, Electrochimica Acta 46, 2229 (2001).
  * Daniel  _et al._ [1987] M. F. Daniel, B. Desbat, J. C. Lassegues, B. Gerand, and M. Figlarz, Journal of Solid State Chemistry France 67, 235 (1987).
  * Thummavichai  _et al._ [2018] K. Thummavichai, N. Wang, F. Xu, G. Rance, Y. Xia, and Y. Zhu, Royal Society Open Science 5, 10.1098/rsos.171932 (2018).
  * Garcia-Sanchez  _et al._ [2013] R. F. Garcia-Sanchez, T. Ahmido, D. Casimir, S. Baliga, and P. Misra, Journal of Physical Chemistry A 117, 13825 (2013).
  * Chen [2011] X. B. Chen, Journal of the Korean Physical Society 58, 100 (2011).
  * zobeir [2019] zobeir, Journal of Physical Chemistry C 123, 23236 (2019).
  * Berg  _et al._ [2006] R. W. Berg, I. M. Ferré, and S. J. Schäffer, Vibrational Spectroscopy 42, 346 (2006).
  * Hoang  _et al._ [2010] L. H. Hoang, N. T. Hien, W. S. Choi, Y. S. Lee, K. Taniguchi, T. Arima, S. Yoon, X. B. Chena, and I. S. Yang, Journal of Raman Spectroscopy 41, 1005 (2010).
  * Maczka [2009] Maczka, Journal of Raman Spectroscopy 40, 2099 (2009).
  * Wang  _et al._ [2016] Y. Wang, S. L. Shang, H. Fang, Z. K. Liu, and L. Q. Chen, npj Computational Materials 2, 10.1038/npjcompumats.2016.6 (2016).