Theoretical investigation of twin boundaries in WO$_3$: Structure, properties and implications for superconductivity

We present a theoretical study of the structure and functionality of ferroelastic domain walls in tungsten trioxide, WO$_3$. WO$_3$ has a rich structural phase diagram, with the stability and properties of the various structural phases strongly affected both by temperature and by electron doping. The existence of superconductivity is of particular interest, with the underlying mechanism as of now not well understood. In addition, reports of enhanced superconductivity at structural domain walls are particularly intriguing. Focusing specifically on the orthorhombic $\beta$ phase, we calculate the structure and properties of the domain walls both with and without electron doping. We use two theoretical approaches: Landau-Ginzburg theory, with free energies constructed from symmetry considerations and parameters extracted from our first-principles density functional calculations, and direct calculation using large-scale, GPU-enabled density functional theory. We find that the structure of the $\beta$-phase domain walls resembles that of the bulk tetragonal {\alpha}1 phase, and that the electronic charge tends to accumulate at the walls. Motivated by this finding, we perform ab initio computations of electron-phonon coupling in the bulk $\alpha_1$ structure and extract the superconducting critical temperatures , $T_c$, within Bardeen-Cooper-Schrieffer theory. Our results provide insight into the experimentally observed unusual trend of decreasing Tc with increasing electronic charge carrier concentration.

We present a theoretical study of the structure and functionality of ferroelastic domain walls in tungsten trioxide, WO3.WO3 has a rich structural phase diagram, with the stability and properties of the various structural phases strongly affected both by temperature and by electron doping.The existence of superconductivity is of particular interest, with the underlying mechanism as of now not well understood.In addition, reports of enhanced superconductivity at structural domain walls are particularly intriguing.Focusing specifically on the orthorhombic β phase, we calculate the structure and properties of the domain walls both with and without electron doping.We use two theoretical approaches: Landau-Ginzburg theory, with free energies constructed from symmetry considerations and parameters extracted from our first-principles density functional calculations, and direct calculation using large-scale, GPU-enabled density functional theory.We find that the structure of the β-phase domain walls resembles that of the bulk tetragonal α1 phase, and that the electronic charge tends to accumulate at the walls.Motivated by this finding, we perform ab initio computations of electron-phonon coupling in the bulk α1 structure and extract the superconducting critical temperatures , Tc, within Bardeen-Cooper-Schrieffer theory.Our results provide insight into the experimentally observed unusual trend of decreasing Tc with increasing electronic charge carrier concentration.

I. STRUCTURE AND PROPERTIES OF WO3
A. Introduction Tungsten trioxide, WO 3 , is a functionally versatile material with possible applications based on electrochromism (smart windows), gasochromism (gas sensors) and photocatalysis [1][2][3] .The high-symmetry structure of WO 3 is that of a perovskite with a vacant A site (see Fig. 1), and it exhibits a series of lower symmetry phases at lower temperature.Both its structure and its properties depend on and can be tuned by doping 4 .For instance, while pure WO 3 is an insulator, it becomes metallic upon occupation of the A sites with alkali metal ions.In addition, superconductivity was reported as early as 1964 for Na-doped WO 3 , and M x WO 3 -x systems are now well established as superconductors (with M usually an alkali metal) [5][6][7][8][9][10][11][12][13][14][15][16][17] .The reported superconducting critical temperatures (T c ) for bulk M x WO 3 -x systems are generally less than 2 K 9,12 , and, interestingly, tend to decrease with increasing doping above the lowest doping level at which superconductivity is observed 9,18,19 .Moreover, different dopants result in superconductivity in different structural phases.More recently, non-bulk high-temperature superconductivity was reported on the surface of dopant-rich islands in Na x WO 3 -x 20 .Around the same time, superconductivity was discovered in reduced WO 3 -x with a superconducting critical temperature of around 3 K [21][22][23] .In this case, the bulk sample was not superconducting, but sheet superconductivity occurred along the ferroelastic domain walls of the reduced WO 3 -x crystals.

B. Structural phase transitions and domain walls
At room temperature, WO 3 shows two types of ferroelastic domain walls which correspond to two successive ferroelastic transitions at higher temperatures.Fig. 2 shows the sequence of structural phases of WO 3 as a function of temperature and doping 21,[24][25][26] .The highsymmetry cubic P m3m reference structure of WO 3 does not form under standard conditions as WO 3 sublimes be-fore reaching it 27,28 .Therefore there are no ferroelastic domain walls resulting from a cubic P m3m to tetragonal P 4/nmm (α 1 ) transition.Additionally, the subsequent α 1 − α 2 transition does not form domain walls since it does not change the point symmetry 29,30 .As temperature is further reduced, the next structural phase transition, and the first ferroelastic transition, is from the tetragonal α 2 (space group P 4/ncc) to the orthorhombic β phase (P bcn).In terms of distortions from the cubic phase, the tetragonal α 2 phase is characterized by two normal modes of the cubic perovskite structure with representations (and wave vectors) M − 2 ( 1 2 , 1 2 , 0) and R − 5 ( 1 2 , 1 2 , 1 2 ), respectively (see Fig. 2 and 3), we refer to these as the cubic M − 2 and the cubic R − 5 modes below.The cubic M − 2 mode (Fig. 3 (a)) consists of antipolar displacements of the W atoms, and the cubic R − 5 mode (Fig. 3 (b)) of out-of-phase rotations of the O octahedra with a 0 a 0 b − Glazer notation.The transition from the α 2 to the β phase then introduces in addition mainly the cubic X + 5 ( 1 2 , 0, 0) and the cubic M + 2 ( 1 2 , 1 2 , 0) modes and reorients the already present cubic M − 2 mode from along the a axis to the (a, b) diagonal spatial direction.The X + 5 mode causes further antipolar displacements of the W atoms and the M + 2 mode introduces in-phase rotations of the O octahedra (Glazer notation a 0 a 0 b + ).The primary order parameter with respect to the tetragonal α 2 phase is an M 1 ( 1 2 , 1 2 , 0) mode, which causes a doubling of the unit cell along the rotation axis, along with antipolar displacements.We refer to the domain walls that form between different orientations of the β phase at this transition as β domain walls.
Upon further reduction of temperature, the second ferroelastic transition corresponds to the transformation from the orthorhombic β phase to the monoclinic γ phase (P 2 1 /n).In addition to the monoclinic distortion of the lattice, this transition introduces additional directional components to the already present cubic X + 5 and R − 5 modes, however the additional amplitude of the X + 5 is very small.These changes require no further doubling of the unit cell, so the order parameter of this second transition is an orthorhombic Γ mode: Γ + 2 (0, 0, 0).We refer to the domain walls that form between different orientations of the γ phase at this transition as γ domain walls.
In experiments, the two types of domain walls create a pattern in which the γ domain walls (blue in Fig. 4) form in a zig-zag manner between the β domain walls (red in Fig. 4) 23,25,31 .The γ domain walls in different β domains meet at 90 degree angles, and the β and γ domain walls are oriented at 45 degrees with respect to each other.Atomic force microscopy (AFM) measurements of epitaxially grown WO 3 films have shown that the two types of domain walls correspond to crystallographic (100) (γ) and (110) (β) planes 31 .The same orientations are implied by a strain analysis of the respective ferroelastic transitions [32][33][34] : The planes of vanishing strain for a ferroelastic transition of type 4/mmmFmmm (which is the type of the α 2 − β transi-P m3m P 4/nmm (α 1 ) Symmetries of, and transitions between, the structural phases of WO 3 as a function of decreasing temperature T or increasing amount of electron doping ρ.The arrows indicate the transition modes using the notation of the structure at the arrow origin.That is blue modes denote modes of the cubic, green of the α 1 , red of the α 2 and black of the β structures.Listed in blue on the right are all modes of the cubic structure present in the respective lower-symmetry structures.The Q and R labels at the transitions represented by red arrows denote parameters of the energy expansions that will be introduced later.
tion) correspond to (110) planes and a transition of type mmF2/m (which the β − γ transition corresponds to) has (100) vanishing strain planes.

C. Questions addressed in this work
In this work, we investigate theoretically the structure and properties of the β domain walls in WO 3 .We start by constructing the Landau-Ginzburg free energy densities using parameters obtained from electronic structure calculations based on density functional theory (DFT), and use these Landau-Ginzburg expressions to calculate the structure and properties of the walls with and without doping.We benchmark our model calculations by also calculating the structure of the β domain wall directly from first-principles using a large supercell.Our calculations allow us to extract both the structural and electrostatic changes associated with domain wall formation.We find that the structure at the β domain wall resembles that of the bulk tetragonal α 1 phase, and that the electrostatic changes cause a local accumulation of electronic charge at the wall.Motivated by these results, we calculate the doping dependence of the critical temperature for electron-phonon mediated superconductivity in the bulk tetragonal α 1 structure, and find that

II. LANDAU-GINZBURG THEORY OF THE TETRAGONAL TO ORTHORHOMBIC (α2 → β) TRANSITION
We begin our treatment of the β ferroelastic domain walls in WO 3 by constructing the equations describing their free energy density according to the Landau-Ginzburg theory of the α 2 − β transition.We then determine the lowest energy domain wall profiles by numerical minimization, and partly by analytical solution, of the free energy density.The detailed analytical solution is presented in appendix section VII A.
As outlined above, the tetragonal α 2 to orthorhombic β transition, which occurs at around 1000 K, is driven by the condensation of a single mode with M 1 symmetry of the tetragonal phase (see Fig. 2 and Fig. 3) 27,29,30 .The irreducible representation of the M 1 mode is twodimensional, and so the order parameter of the transition, which we denote as Q, has two components, Q(q 1 , q 2 ).The Landau free energy density F Q of a domain wall described by Q, then depends on both q 1 and q 2 , as well as on the strain, e.It is given by the following expansion around the α 2 bulk free energy density 1 + q 2 2 ) 3 + dq 2 1 q 2 2 + e(q 4 1 q 2 2 + q 2 1 q 4 2 ) + λ 1 (q 2 1 + q 2 2 )e s + λ 2 (q 2 1 − q 2 2 )e as + λ 3 (q 2 1 + q 2 2 )e 3 + s[(∇q 1 ) 2 + (∇q 2 ) 2 ] + C ij e i e j , (1)   where we use Voigt and Einstein notations 35,36 .F Q , q 1 , q 2 and the strains are treated as continuous fields and the spatial coordinates z (for example q 1 (z), etc.) are implied.The non-symmetry breaking (e s ) and symmetrybreaking (e as ) strains can be related to the amplitudes of the in-plane eigenvectors of the tetragonal elastic tensor, e 1 and e 2 , through e s = e 1 + e 2 and e as = e 1 − e 2 , and e 3 is the strain along the tetragonal axis.Parameters a, b and c describe the Landau potential up to sixth order in q 1 and q 2 .Parameters d and e describe the additional coupling between order parameter components q 1 and q 2 up to sixth order that is not contained in b and c.The parameters λ 1 , λ 2 and λ 3 describe the separate couplings between the order parameter and the tetragonal strains.The Ginzburg parameter s accounts for the variation of the order parameter Q in the domain wall.The last term is the strain energy with C ij the elastic tensor.To circumvent the explicit calculation of the strain dependence, we incorporate the energy-minimizing strains in effective Landau parameters for q 1 and q 2 .The general free energy density for a domain wall described by Q then simplifies to Parameters a Q , b Q and c Q correspond to terms that are non-vanishing even when Q has only one component (i.e.q 2 = 0) and parameters a 2Q and b 2Q describe the bidirectional coupling between the order parameter components that only occur when both q 1 and q 2 are non-zero.
A. Extension of the Landau-Ginzburg free energy density to include the effect of additional charge Next we extend the Landau-Ginzburg free energy density to study the effect of additional charge, ρ(z), introduced by reduction or doping.
a. Effect of charge on the order parameter, Q.We begin by analyzing the effect of additional charge density on the Q order parameter describing the α 2 to β transition by extending the free energy density expression as follows.
Here, the direct effect of the charge on the free energy density in the α 2 reference structure appears explicitly as the chemical potential term µ(ρ).All additional effects of the change in chemical potential are incorporated in the ρ dependence of the Landau parameters, a Q (ρ) etc.Note that we also account for the effect of charge on the gradient parameter s(ρ).
b. Effect of charge doping on the amplitude of the cubic R − 5 mode.In addition to affecting the Q order parameter responsible for the α 2 to β transition, the addition of charge has the effect of reducing the amplitude of the cubic R − 5 mode which is present in both the α 2 and β phases 37,38 .This mode is the order parameter for the transition between the α 1 and α 2 structures (see Fig. 2).Complete suppression of the cubic R − 5 mode therefore transforms the α 2 phase to the higher symmetry α 1 phase.In order to take this into account, we extend the Landau potential for α 2 further by expanding this mode around its value in the α 2 phase, R α2 .For convenience, we define an expansion parameter R = R α2 − |R − 5 |, which is zero in the α 2 phase and increases with doping, reaching the value R α2 in the α 1 phase, and so its sign matches that of a conventional Landau theory order parameter.(|R − 5 | is the amplitude of the cubic R − 5 mode at the particular doping value of interest).The new parameter R therefore describes the reduction in the amplitude of R − 5 in the transition from α 2 to α 1 .
Including this degree of freedom in the Landau potential with this definition of R leads to the free energy density: Thus the additional energy, F R = F QR − F Q , is 0 in the α 2 phase where R is 0. Note that we included R2Q terms, in which R is present in a single and Q in two directions, up to eighth order.

III. COMPUTATIONAL DETAILS A. Choice of exchange-correlation functional
The properties of WO 3 are unusually sensitive to the choice of exchange-correlation functional, with many studies in the literature suggesting different choices.Consistently good matches of relaxed structures to experimental structures have been reported using the B1-WC hybrid functional by Hamdi et al. 39,40 , as well as by Wang et al. using HSE-06 albeit not to the same degree 38 , but use of a hybrid functional is prohibitively expensive for our calculations.We found that the generalized gradient approximation (GGA) in the PBEsol implementation grossly underestimates the amplitude of the tetragonal M 1 mode; a similar underestimation of the oxygen rotations in GGA(PBE)-relaxed monoclinic WO 3 has also been reported 41 .A more detailed comparison with published calculations is often problematic, since in many cases only the lattice parameters of relaxed WO 3 bulk structures are reported but not the internal coordinates 19,37,[42][43][44][45][46][47][48][49][50][51][52][53] .
In this work we use the local-density approximation (LDA) description of the exchange-correlation functional.Our motivation is its good description of the amplitude of the tetragonal M 1 mode, which is the order parameter Q of the α 2 − β transition that we study in detail here.The lattice constants and phonon mode amplitudes of WO 3 bulk structures that we calculate within the LDA in this work are listed in TABLES I and II.
We note, however, that our chosen LDA implementation is not suitable for describing the γ domain walls because it does not yield a pronounced and necessary decrease in energy from the β to the γ phase as for instance reported by Hamdi et al 39 .A detailed discussion of this point is provided in the appendix section VII C.

B. Calculation of Landau-Ginzburg parameters
The calculations to obtain the parameters of the Landau-Ginzburg free energies in Eqns.( 3) and (4) were performed using the Quantum Espresso (version 6.2.1) plane-wave pseudopotential DFT implementation 54,55 .The choice of DFT implementation was made to be consistent with the electron-phonon coupling calculations, which we describe later.This forced us to use normconserving pseudopotentials as these were the only available option for electron-phonon calculations when this work was started.The norm-conserving LDA pseudopotentials were generated with the ONCVPSP program and the input parameters provided by the PseudoDojo pseudopotential repository [56][57][58][59] .A high cutoff energy of 120 Ry was necessary to converge the parameters in the Landau potentials, due to the use of norm-conserving pseudopotentials and the small core of the available W pseudopotential.We used valence electron configurations of 4f 14 5s 2 5p 6 5d 4 6s 2 for the W atoms and 2s 2 2p 4 for the O atoms.The k-and q-point grid sizes were set to 12 × 12 × 12 and 4 × 4 × 4 respectively in the cubic phase, and scaled down relatively for larger unit cells, ensuring that they were always commensurate with each other as required for the electron-phonon calculations.
The Landau parameters were determined by calculating the energies of structures with different amplitudes of the Q (q 1 and/or q 2 ) and/or R distortions frozen into the reference α 2 structure in the range from 0 to 1.2 Å per unit cell.We calculated the energies of a total of 360 distinct (R, q 1 , q 2 ) points, with the size of the unit cells allowed to relax in each case to satisfy the condition of energy-minimizing strain.
Parameters to sixth order (eighth order for the R2Q) were calculated for all described terms, with higher-order terms, constrained to be small and positive, included in the fit in each case to prevent unphysical negative divergence.To calculate the change of the parameters on charge doping, we repeated the set of calculations for a total of four different amounts of additional electrons up to 0.25 electrons per f.u.The cell parameters were set to those obtained from relaxations that did not contain additional charge and they were not allowed to relax further.The changes in the parameters were then fitted up to quadratic order of the charge density ρ (see TABLE III).
We fitted the chemical potential µ(ρ) of the α 2 reference structure with a quadratic dependence on the charge density ρ in the free energy density (see TABLE III).To accurately extract a value for the quadratic term, which is small compared to the linear term, a total of 80 energies for charge densities between 0 and 0.25 electrons per f.u. were calculated.
Finally, the gradient parameters s and t were obtained with the procedure described in the appendix section VII B. The inter-atomic force constants of the α 2 phase in real space were calculated by interpolating the dynamical matrices on a q-grid.Force constant matrices were then interpolated for q-points on the q-paths (1/2, 1/2, 0) → (1/2 + δ, 1/2 + δ, 0) for the M 1 and (0, 0, 0) → (δ, δ, 0) for the Γ + 2 mode (the corresponding α 2 − β domain walls correspond to the (110) crystallographic plane).The path length δ was set to 0.04.The branches belonging to the transition modes M 1 and Γ + 2 in the force constant matrix dispersion were determined by symmetry combined with visual analysis of the respective displacements η.The corresponding gradient parameters s and t were finally obtained by performing a quadratic fit to the determined force constant branch as shown in equation (18).To describe the change of the gradient parameters with charge, calculations of s(ρ) and t(ρ) were performed in the α 2 cell with three values of ρ, and then fitted to fourth order in ρ with the third order term omitted (see TABLE III).

C. Supercell calculations of domain wall structures
In addition to our Landau-Ginzburg calculations of the β domain wall, we also performed direct calculations by explicitly relaxing the domain wall structure using DFT.We constructed a supercell containing two β domain walls corresponding to (110) crystallographic planes for subsequent relaxation (see Fig. 5).The supercell was generated as follows: First, the bulk structures were relaxed and one domain was constructed with the resulting relaxed structure.The second domain was then created by application of the point-group symmetry operations on the first domain that are lost during the transition.The supercell for the β domain wall calculation contained 512 atoms and had dimensions of approximately 84 × 11 × 8 Å.
We then relaxed the supercell structure with some atoms fixed to the bulk structure (see Fig. 5) as releasing the bulk cells causes a relaxation back to a single domain state.The relaxed structures were analyzed in terms of distortion modes.Additional charge was then introduced to the relaxed cells to determine if there was an accumulation of charge at the domain walls.
All calculations for the large domain wall supercells were performed with the GPU-accelerated VASP (version 5.4.4)DFT implementation [60][61][62][63] .The wave functions were expanded with a basis set of plane waves and their cutoff energy was set to 800 eV for the bulk structures but was decreased to 600 eV for the supercell relaxations.
The core electrons were treated with projector-augmented waves, allowing for a lower planewave cutoff energy compared to the norm-conserving pseudopotentials employed for the rest of this work 64,65 .The valence electron configurations were 5p 6 5d 4 6s 2

D. Calculations of electron-phonon coupling and superconducting critical temperature
First-principles calculations of electron-phonon coupling were performed for electron-doped 8-atom α 1 and 16-atom α 2 cells.Electron doping was achieved by adding electrons with a compensating background charge rather than explicit inclusion of point defects; this method has been shown to describe well the chargeinduced structural distortions in WO 3 37 .For reasons discussed in the next subsection III E, we only relaxed the internal coordinates and the cell parameters were manually set by linear interpolation between the calculated values for the tetragonal WO 3 and the cubic NaWO 3 cells.
We used the EPW package in conjunction with Quantum Espresso 67,68 .The cutoff energy had to be kept at an extremely large value of 120 Ry to converge the calculations, as tested by the convergence behaviour of the total electron-phonon coupling strength λ in the cubic phase.Coarse k-and q-point grids were set relative to 12 × 12 × 12 and 4 × 4 × 4 grids of the cubic phase.Additionally, fine grid sizes for the k-and q-point grids were set relative to 200,000 and 100,000 random points for the cubic phase, respectively.W 5d xy , 5d xz and 5d yz orbitals were chosen for the wannierization procedure, which was performed using the wannier90 package 69 .Superconducting critical temperatures were then extracted using the usual Eliashberg formalism of Bardeen-Cooper-Schriefer (BCS) theory, as described in Appendix section VII D. The temperature for the Fermi occupations in Eqn.(20) was set to 0.075 K and the Fermi surface energy window of considered electron states was set to 3 eV.The Coulomb pseudopotential parameter µ was set to 0.10.All results are given for smearings of 0.05 meV for the delta functions in Eqn. ( 22) and 0.25 eV for the frequency delta functions in Eqn.(26), respectively.None of the calculations were performed with the double-delta approximation.

E. Lattice relaxation of charged unit cells
As the calculated stresses in charged unit cells with a constant background charge are not well-defined in DFT implementations, we did not relax the volumes of our unit cells and supercells in the calculations for which we include additional electronic charge 70 .
To determine the validity of keeping the lattice parameters fixed in our Landau model, we compared our Landau parameters calculated using the relaxed lattice parameters of the undoped cells with calculations in which the cell parameters of the α 2 cells were interpolated to those of the cubic NaWO 3 cell.We found that the changes in the Landau curves caused by the volume change were negligible compared to those caused by the introduction of charge into the uncharged cells.Thus, the change in volume caused by doping within the doping range used here can be safely disregarded.
In the case of the electron-phonon calculations for the α 1 phase, the change in cell size on doping could not be disregarded, as the phonon dispersions depend strongly on the cell volume.For example rotational modes are artifically stabilized if the cell volume is not allowed to  increase, whereas antipolar modes are artifically destabilized.Thus, we opted for the compromise of linearly interpolating the lattice parameters between those calculated for the tetragonal WO 3 and for the cubic NaWO 3 unit cells.We found this approximation to be sufficient to describe the correct general trends of the modes upon doping in the α 1 phase.In particular, the rotational R − 5 mode became softer with decreasing charge x so that the α 2 − α 1 transition occurred at a doping value close to the experimentally observed transition value of x = 0.2 in Na x WO 3 -x 9 .Also, the amplitude of the antipolar cubic M − 2 mode decreased with increasing x consistent with the literature 37 .

IV. RESULTS OF DOMAIN WALLS CALCULATIONS A. Landau-Ginzburg domain wall profiles
Using the Landau-Ginzburg parameters obtained as described in Section III B, we calculated the profiles of the order parameters across the domain walls by numerically minimizing the total free energy density functional given in equation ( 4).The spatial grid of the order parameter fields, z, consisted of 251 points, spaced by ∆z = 0.6 Å. Self-consistent solutions were found as follows: For given q 1 (z) and q 2 (z) profiles at a specific total charge, the minimum energy charge distribution ρ(z) was calculated.For this ρ(z), the minimum energy R(z) and subsequently q 1 (z) and q 2 (z) were obtained, after which the cycle was repeated.Self-consistency was achieved  when the change in the total z-integrated energy density between steps was less than 10 −4 meV Å−2 .We checked our numerical approach for the simplified potential of eq. ( 7) by comparing with the analytical solutions of equations ( 9) and (12), and found excellent agreement between the numerical and analytical results.
The calculated evolution of the order parameters across the energetically minimized β domain walls for various doping levels is shown in Fig. 6.For clarity, we plot the magnitude |Q| = q 2 1 + q 2 2 and the angle φ = arctan q 1 /q 2 of the order parameter.Panel (a) shows the angle φ (solid lines) with respect to the z-axis along the wall.We see that, for the undoped case (blue line), the wall shows characteristic Néel-like behavior.The order parameter Q retains 80 % of its bulk amplitude across the wall and the reorientation is achieved by rotation of Q along the wall as represented by φ.The small reduction in amplitude at the wall can be explained by the bidirectional coupling of Q (that is q 1 to q 2 ) which results in an energy reduction when the amplitude of Q decreases.We extracted the wall width 2ξ by fitting Q(z) to a | tanh(z/ξ)| curve as in eq. ( 9) and obtained a value of ∼1.35 nm for the undoped case.As can be seen from equations (10) and ( 14), the widths are mostly determined by the gradient parameter s of the order parameter Q which is an order of magnitude larger than the Landau terms (see TAB. III).This value lies well within the general range of ferroelastic domain wall widths, which are generally between 0.2 and 2 nm at low temperatures 72 .
A distinct change in behavior is seen on introduction of electrons.At the most strongly doped example studied, 0.24 electrons per formula unit (yellow line), the wall is strongly Ising like, with the amplitude of Q suppressed to zero in the wall region.At this highest doping level, the domain wall width is widened by a factor of around 3.6 relative to the width in the undoped wall as measured by the fitting of the φ curves to | arctan(exp(z/ξ))| curves.
The crossover from undoped behavior to doped behavior, as well as the wall broadening, are gradual, with intermediate dopings (purple, red and orange colors) having intermediate behavior.
The origin of the evolution with doping is clear in Fig. 6 panel (b) which shows the charge density as a function of position across the wall.We see that, for all doping levels, the charge accumulates in the wall region, and no additional charge remains in the bulk of the domains.
As discussed earlier, electron doping causes a reduction in the amplitude of the Q order parameter, moving the structure towards the α 2 phase.In addition, it causes a decrease of the Z + 3 mode, parametrized by |R − 5 | as shown in Fig. 6 panel (c).As a result the structure within the domain wall approaches that of the α 1 phase.Note that we calculated Landau-Ginzburg parameters only for concentrations up to −2.5 me Å−3 .This is the origin of the forced cutoff of ρ in the yellow curve in Fig. 6 panel (b).However, we expect that |R − 5 | would decrease to 0 with increasing doping.

B. Direct calculation of domain walls using density functional theory
Motivated by our estimation of the domain wall widths of ∼14 Å from our Landau-Ginzburg model, we next performed a full density functional calculation of the domain wall structure shown in Fig. 5.
Our calculated layer-resolved order parameter angle φ and its amplitude |Q| are shown in Fig. 7 panel (a) as a function of position perpendicular to the wall plane, z.Consistent with our results from Landau-Ginzburg theory, we find the domain wall to be predominantly of the Néel type as represented by the gradual transition in φ.In Fig. 7 panel (b) we show the calculated charge density distribution obtained by adding an additional but small electronic charge of 1.6 × 10 −2 e to the supercell.We find that, as in our Landau-Ginzburg simulations, the charge accumulates at the walls.While the macroscopic planar charge density shows an alternating behaviour from site to site, there is a clear depletion of charge from the bulk towards the domain wall structure as indicated by the top and bottom envelopes of the density.We note that a definitive study would require further relaxation of the wall after the introduction of the charge.In addition, a systematic study of larger supercells would be desirable to ensure that there are no interactions between the walls and that full convergence to the bulk values is achieved in the intermediate regions.The supercells used here were barely large enough to host two domain walls as we observe no clear bulk plateau in the order parameters.
A clear difference compared to the Landau-Ginzburg model can be found in the amplitude of the Q (M 1 ) mode in the domain walls.In the DFT-calculated walls we observe a slight increase in Q, whereas it decreased in the walls obtained from Landau-Ginzburg theory.The reason for this is the limitation of the phase space used in the Landau-Ginzburg model.When comparing fully relaxed structures of the P bcn (β) (with order parameter direction Q(a, 0)) and the P 42 der parameter direction Q(a, a)) phases we observe that the latter has a lower energy and a higher total Q amplitude than the former.Thus, we expect that including additional order parameters in the Landau-Ginzburg model would also lead to an increase in Q in the center of the domain wall.However, this would lead to a highly increased dimension of the phase space, making the Landau-Ginzburg parameterization unfeasible.Due to the small amplitude of these additional distortions, it is reasonable to assume that the difference in domain wall width would be minor if they were included.Furthermore, we expect the inaccuracy in the Q displacement to be less relevant in the charged domain walls, as the charge reduces the amplitude of both Q and R − 5 at the domain wall.

C. Summary of domain wall results
In summary, we investigated the structure of the βtype domain walls that form during the phase transition from the tetragonal α 2 phase to the orthorhombic β phase in WO 3 in the framework of Landau-Ginzburg and density functional theories.Our Landau-Ginzburg calculations showed that the ferroelastic walls in the undoped case are mostly Néel-like, with the amplitude of the order parameter Q retaining ∼80% of its bulk value.We found the domain wall width 2ξ to be around 14 Å in the undoped case.Electronic doping increased our calculated domain wall width and led to an accumulation of the additional charge in the domain wall.The domain wall width of WO 3 β domain walls at very low temperatures has been reported to be around 2w = 1.2 nm in experiments 24,25 , which is very close to the value suggested by our Landau-Ginzburg model.We found that the charge accumulation at the walls caused an increasingly large Ising-type component, indicated by the drop in the order parameter amplitude |Q| across the wall.The accumulated charge also reduced the magnitude of the |R − 5 | order parameter, so that the structure approached that of the α 1 phase in the wall region.Using DFT calculations on supercells, we were able to confirm the Néel-type character of the domain walls, as well as the predicted accumulation of charge at the domain walls.

V. IMPLICATIONS FOR DOMAIN WALL SUPERCONDUCTIVITY
Motivated by our finding that the charge accumulates at the domain walls and causes a local α 1 -like structure, we next study the superconducting properties of this phase.Experimentally, the α 1 phase of Na x WO 3 -x , which also has P 4/nmm symmetry 73 , was shown to be superconducting.Similar to other superconducting tungsten bronzes, the superconductivity shows two general features.First, for each dopant type, superconductivity occurs in only one high-symmetry structure.For smaller alkali metals (Na and K) these are structures of tetragonal symmetry, while for larger alkali metals (Rb and Cs) the structures are hexagonal.At doping levels x that lie above or below the x-range of these phases, superconductivity is not found.The second feature is a decrease in T c with increasing x, within the superconducting phase.Thus, the highest T c is reported at the lowest x-value at which the superconducting phase is still retained; lower doping results in a phase transition to the lower-symmetry, non-superconducting phase.Both properties implicate the soft mode associated with the corresponding structural phase transition in the superconductivity mechanism 6,9,21,74 .
Our approach is to calculate and analyze T c as a function of doping within standard Bardeen-Cooper-Schrieffer (BCS) theory 75 for the α 1 phase of WO 3 .While BCS theory has been shown to capture some as-pects of the behavior of doped WO 3 9,22,76,77 , the absence of superconductivity in the α 2 phase, and the decrease in T c with increasing doping are not well understood (the latter has even been described as the "T c paradox"! 9,14,16,19,74,78,79), and we explore these aspects here.
We calculate the electron-phonon coupling matrix, using density functional perturbation theory (DFPT) 80 .
Here, ∆ ν q V KS is the phonon perturbation to the Kohn-Sham potential, and the matrix elements are the transition probability amplitudes for an electron in initial state ϕ nk with wave vector k and band n, scattering to final state ϕ mk+q of band m, via a phonon of wave vector q and branch ν.We then evaluate the superconducting critical temperature using the semi-empirical Allen-Dynes equation 81 with the coupling strength, λ, and the weighted phonon frequency ω , extracted from the electron-phonon matrix as described in the appendix section VII D and an empirical value of 0.1 taken for the effective Coulomb repulsion µ.
A. BCS theory applied to the bulk α1 phase We begin by calculating the superconducting T c for the α 1 phase, to see whether the measured decrease in T c with increasing doping is correctly captured within BCS theory.The α 1 structure is stable for calculated electron concentrations larger than x ≈ 0.125 (for lower concentrations, it has an unstable Z + 3 mode, indicating the transition to the lower-energy α 2 structure), which is therefore the lowest doping concentration that we consider.Our calculated electron bands, phonon bands, phonon linewidths, phonon density of states and Eliashberg spectral function α 2 F for x = 0.125 are shown in Fig. 8, where the soft Z + 3 mode is the lowest frequency Z mode in the phonon bands close to zero frequency.Interestingly, while there is some electron-phonon coupling at low frequency, it is considerably stronger at higher frequencies, with the highest values of α 2 F occuring at around 600 -800 cm −1 .This suggests that, at least in the BCS picture, the soft mode is not the most relevant in determining the superconducting T c .The subsequent changes in the phonon bands and α 2 F upon increase of x are presented in Fig. 9.We see that, as expected, the Z + 3 soft mode hardens with increasing doping, leading to a reduction of α 2 F at low frequency with increasing doping.Interestingly, the high energy phonons shift to lower frequencies as doping is increased, with corresponding shifts of the peaks in α 2 F to lower frequency.Finally, the calculated superconducting critical temperature, T c , Fig. 10 clearly illustrates that our BCS-theory calculations reproduce the experimental trend of decreasing T c with increasing doping, with the calculated maximum in T c at x = 0.125 coinciding with the calculated α 1 − α 2 transition, where the transition Z + 3 mode starts to become imaginary.In addition, our calculated T c 's are comparable to the reported values ( < ∼ 2 K in bulk samples), although we emphasize that their actual magnitudes should not be over-interpreted, since they are sensitive to the spreads in α 2 F integration in equation (26)  and the value of the screened Coulomb potential µ in the Allen-Dynes formula (6).The trend of a decreasing for the α 1 structure, for a range of added electron concentrations.Blue corresponds to the lowest (0.125 electrons per formula unit) and orange to the highest (0.24 electrons per formula unit) doping levels, with successive lines corresponding to successive points in the plots of Fig. 10.
0.12 0.14 0.16 0.18 0.20 0.22 0.24 x [e/f.u. ] T c , however, is robust to these parameters.Therefore we conclude from our calculations that conventional BCS theory captures the observed evolution of T c with doping in the α 1 phase of WO 3 . 83iven the good agreement of our computational BCS theory results with experiments, we next analyze them to rationalize the behavior.In particular, the decreasing T c with x was unexpected within a simple BCS picture, since the electron density of states at the Fermi level, n(E F ), has been reported from photoemission measurements to increase with increasing x in tetragonal Na x WO 3 -x 84,85 .Since the BCS Cooper-pair binding energy scales as exp [−1/(n(E F )V ep )] (V ep is the interelectronic attraction caused by the electron-phonon coupling), an increase in n(E F ) should in turn lead to an increase in T c , provided that the electron-phonon coupling strength remains constant with electron density.Our calculations of n(E F ) (Fig. 10) in fact indicate that, within the α 1 phase, n(E F ) (red line) at first decreases with increasing doping (as does T c ).Note that there is no inconsistency with Refs.84 and 85, which provided measured n(E F ) values only above x = 0.25.At higher dopings (above x = 0.18) n(E F ) starts to increase, while T c continues to decrease.This lack of correlation between T c and n(E F ) points to a doping dependence of the electron-phonon matrix elements.
In Fig. 11 we show our calculated doping dependence of T c and total coupling strength, λ, as defined in equation (25).First we note that, over the whole range, the coupling strength, λ, has a substantial value, consistent with the measurable superconductivity.Second, as we expected, it is clear that the value of λ decreases with increasing doping, explaining the corresponding decrease in T c according to the Allen-Dynes formula given in equation 6.In particular, the calculated T c tracks closely the calculated value of λ.
Finally, to understand the change in superconductivity across the α 1 − α 2 transition, we performed a calculation of T c in the α 2 phase where experimentally superconductivity has not been measured.We chose a value of x = 0.125 for the T c calculation in the α 2 phase, and adjusted the lattice constants as outlined in section III E for the α 1 phase, so that the resulting system in the α 2 phase was quite far from the α 1 −α 2 transition.(The amplitude of the transition mode Z + 3 was 0.51 Å, compared with the amplitude of 0.76 Å in the undoped α 2 phase).As expected, our calculated T c dropped sharply from the calculated α 1 value, to 0.018 K, consistent with a sharp drop in the calculated λ value to 0.25.The bands and linewidths for the α 2 case are shown in appendix Fig. 15.

B. Discussion of bulk superconductivity results
While our calculations indicate that the superconducting behavior of WO 3 can be reproduced within standard BCS theory, this is of course not definitive evidence that WO 3 is a BCS superconductor.In this section we discuss two other models for superconductivity -based on soft modes and bipolarons respectively -that have been discussed in the literature.
We begin with a discussion of the importance of the soft mode, whose strong change in frequency with doping was originally proposed to account for the apparently paradoxical behaviour of T c upon doping 9 , in spite of its absence in inelastic neutron scattering experiments 74,78  FIG.
12: Calculated T c as a function of the upper integration frequency limit, ω max , in equations ( 25) and (27).Blue curves show the lowest and orange curves the highest doping x (in electrons per formula unit), over the same range as in Fig. 10.
and doping dependence of α 2 F point to a small, if any, role of the soft mode at the BCS level; here we quantify its contribution.In Fig. 12, we plot the calculated T c as a function of the maximum frequency of the phonons included in the calculation, for a range of doping values within the α 1 phase.We see that the modes below 200 cm −1 , which include the soft mode, contribute negligibly to T c .As noted above, the modes above around 600 cm −1 , which had the largest phonon bandwidth and the strongest frequency shifts on doping, contribute most to T c at every doping concentration. 86While our calculations were performed for tetragonal structures, we note that Brusetti et al. have attributed the increase in T c with decreasing x in hexagonal Rb x WO 3 -x to changes in electron-phonon coupling for phonons with a frequency of more than 240 cm −179 .Before leaving the topic of soft-mode superconductivity, we point out that the superconductivity in WO 3 is somewhat reminiscent of that in SrTiO 3 , for which a model of superconductivity mediated by fluctuations associated with the ferroelectric quantum critical point has been proposed 87 .The quantum criticality model had considerable success in reproducing the measured behavior, as well as in making rather bold predictions about strain and isotope effects that were subsequently verified experimentally [88][89][90][91] .An important difference is that SrTiO 3 has a superconducting dome as a function of doping, whereas in WO 3 an analogous picture would have the left side of the dome cut off due to the absence of superconductivity in the α 2 phase.If this mechanism is relevant in WO 3 , a large and anomalous oxygen isotope effect on T c should be observed.
That the carriers in doped WO 3 form bipolarons at low temperature was first proposed by Schirmer and Salje in 1980 to explain optical absorption, conductivity and electron spin resonance data 92,93 .The idea seems to have been largely neglected until a recent measurement of sheet superconductivity in shear planes of WO 3 -x attributed the remarkably high reported T c of 80 K 94 to the W 5+ -W 5+ bipolarons observed using electron paramagnetic resonance.Bipolarons have also been speculated to be responsible for the high temperature surface superconductivity in H x WO 3 20 .A recent density functional study of self-trapped polarons in WO 3 succeeded in capturing a polaronic state, with substantial lattice distortions, in the simulations 95 , although the polaron was at higher energy than the delocalized electron.The role of electron localization and its coupling to the lattice is clearly an important area for future study 72 .

C. From bulk to sheet superconductivity
To link our bulk WO 3 results to the sheet superconductivity reported at the WO 3 β domain walls 21,22 , we revisit the β domain walls that we obtained from Landau-Ginzburg theory and our density functional calculations.We can make three main inferences.First, from both studies we see that it is lower energy for electronic charge to be at the domain than in the surrounding β bulk structure leading to local charge accumulation at the walls.Second, this local increase in charge induces a local transition to the tetragonal bulk α 1 phase in the domain walls.And third, the additional charge, combined with the presence of the α 1 phase, leads to strong enough electron-phonon coupling to enable superconductivity in the domain walls.
Many of the samples in which domain wall superconductivity was measured showed a stripe pattern of parallel ferroelastic domain walls of only one type (see for example Ref. 22).Since similar samples were characterized in detail and shown to consist of β domain walls 24 , it is likely that the superconducting samples contained only β domain walls.Whether γ domain walls are also superconducting, and if so by what mechanism, is an interesting open question for future study.

VI. CONCLUSION
In this work, we calculated the structure and properties of the ferroelastic domain walls within the β phase of WO 3 , using a combination of first-principles density functional calculations and Landau-Ginzburg theory.We showed that the ferroelastic β domain walls have mixed Néel and Ising character, and found that free electronic charge preferentially accumulates at the domain walls.We showed that this accumulation of charge leads to a broadening of the walls and an increase in their Ising character, as well as a change in the atomic structure within the domain wall structure to the α 1 (P 4/nmm) phase.This latter phase is known to be the superconducting phase in doped WO 3 , suggesting that the domain wall superconductivity is a consequence of the combined electron accumulation and local structural change at the walls.
To investigate further this possible link between domain wall and bulk superconductivity, we performed electron-phonon calculations based on DFT to calculate the T c as a function of doping in the bulk α 1 phase at the BCS-theory level.Our calculated values were comparable in magnitude to the measured values ( < ∼ 2 K) and showed the same trend of decrease in T c with increasing doping.The evolution of T c with doping correlated with a reduction in the electron-phonon coupling, with the largest contribution coming from high frequency phonons above approximately 600 cm −1 .
Our calculations suggest that the superconductivity at the domain walls in WO 3 results from the combined accumulation of charge at the walls and the structural changes at the domain walls that are induced by the presence of the carriers.
Using the Landau-Ginzburg free energy density expression, one can calculate domain wall profiles by minimiz-ing the free energy density with appropriate boundary conditions.For a 2D order parameter, in general, we can have two types of domain walls, Néel-type (order parameter rotates along the wall) and Ising-type (order parameter vanishes on the domain wall).These two limiting cases can be calculated analytically for a simple Landau theory of the form: 1. Néel wall For ease of calculation, we parametrize q 1 and q 2 with polar coordinates {Q, φ}.For a fixed amplitude Q 0 , we construct the Euler-Lagrange equation ( 2) with respect to φ(z) and obtain the following equation: Using the boundary conditions 8) is solved by a stationary Sine-Gordon equation and its solution is given by with ξ considered to be half the domain wall width.

Ising wall
The other possible domain wall for a 2D order parameter is the Ising-type wall, in which the order parameter amplitude vanishes in the middle of the wall.(Note that for 1D irreps, this is the only possibility.)For such a wall, the Euler-Lagrange equation has the form where the 6th-order Landau term has been omitted, so that we can exploit the known solutions of the 4rd-order equation.These solutions are

General domain wall profile
In reality, a structural domain wall with 2D order parameter will be a mixture of Néel-and Ising types, and its profile can be obtained by solving both ( 8) and (11)  simultaneously.This problem can likely not be solved analytically.

B. Determination of the gradient parameter
To calculate the gradient parameters we followed the procedure outlined by Artyukhin et al in Ref. 96.Consider any order parameter Q described by the eigendisplacement η Q (q) of a force constant mode corresponding to a wave vector q.Q is then given by an eigenvector of the force constant matrix The gradient energy term in q-space associated with the parameter s Q can then be written as which equals 0 if q = 0. Therefore it is possible to determine s Q for some direction of q by calculating the energies of supercells with distortions described by Eqn. 15 with various magnitudes of qs frozen in along this direction.
A more feasible approach is to determine s Q from the force constant dispersion.The Hessian of the gradient energy for all modes η q of q is equal to the force constant matrix in q-space C (q) within the harmonic approximation: Consequently, the eigenvalues of the Hessian in expression (17) are the eigenvalues of the force constant matrix and we can determine all gradient parameters of modes with wave vector q by expanding the Hessian in (17)  around q = 0: The gradient parameter corresponding to the mode η Q is then given by the eigenvalue of expression (18) corresponding to the mode η Q : C. Exchange-correlation-functional suitability for the description of the β to γ transition As mentioned in the main text, both LDA and GGA (PBEsol) exchange-correlation functionals yielded almost identical energies for the β and γ phases, even though the structural relaxations yielded distinct structures.This was the case for both Quantum Espresso and VASP calculations.Consequences of the small energy difference between the two phases were a negligible Landau parameter a, which (consistent with Eqn.14) led to unreasonably large widths for the γ domain walls, and erratic behavior in the DFT structural relaxations.
A crude estimate for the Landau parameters can be made by using the β − γ transition energy obtained with B1-WC calculations as reported by Hamdi et al. to approximate the Landau parameters a and b 39 .The condition that the mode amplitude Q 0 in (13) equals the bulk amplitude of the orthorhombic Γ + 2 mode in the γ phase (experimentally 0.531 Å) and that the Landau potential in eq. ( 2) with only the a Q and b Q terms equals this reported energy difference (F 0 P bcn -13 meV/f.u.) at Q 0 results in the values for a Q and b Q of −0.22 meV Å−5 and 0.39 meV Å−7 which are comparable to our values for the α 2 − β transition Landau potential 27,39 .However, ultimately a final estimation for the Ising domain wall width cannot be made without the Ginzburg parameter s and the calculation thereof for the orthorhombic Γ + 2 mode does not seem reasonable based upon the poor characterization of the γ phase by LDA.

D. Superconducting critical temperature
Once the electron-phonon matrix in equation ( 5) is known, the phonon linewidth γ qν resulting from electronphonon interaction can be calculated 67,68,81,[97][98][99] .Within the Migdal approximation the linewidth of a phonon with wave vector q and branch ν is given as the imaginary component of the phonon self-energy, where w k denotes the weights for the k-points, nk is the band energy, f ( ) is the associated Fermi occupancy, ω ν q is the phonon frequency and η is a smearing parameter for allowed transitions.In principle, the latter can be neglected in calculations where the k-and q-grids are dense enough.In such a limit of vanishing smearing, and additionally vanishing phonon frequencies, lim η,ωqν →0 γ qν , one arrives at the so-called double-delta approximation of the phonon linewidth γ ν q = 2πω ν q nm 1 Ω BZ BZ dkw k |g ν nm (k, q)| 2 δ( nk )δ( mk+q ), ( 21) where a smearing may be reintroduced in the the two delta functions.A similar expression as the one for the double-delta approximation (21) gives the electronphonon coupling strength λ ν q for the phonon dk|g ν nm (k, q)| 2 δ( nk )δ( mk+q ), (22)  where n(E F ) is the DOS at the Fermi energy.Consequently, the coupling strength in the double-delta approximation is given as Calculating the Brillouin-zone average of the coupling strength yields the first parameter in the McMillan formula which is the total coupling strength λ where w q now denotes the weights for the q points.In the Allen-Dynes formula, λ is evaluated as where α 2 F (ω) is the isotropic Eliashberg spectral function, which in turn is given as dqw q ω ν q λ ν q δ(ω − ω ν q ).( 26) The delta function may again be subject to a smearing for numerical calculations.The remaining characteristic phonon frequency ω , according to Allen and Dynes, is then given as ω = exp 2 λ dω ω α 2 F (ω)log(ω) .

FIG. 1 :
FIG. 1: High-symmetry cubic pseudo-perovskite P m3m structure of WO 3 .The perovskite B sites are occupied by W atoms (grey) which are encapsulated by O atom (red) octahedra.The A sites on the unit cell corners are vacant.

2 FIG. 3 :
FIG. 3: Distortions of the cubic WO 3 structure, labelled with the irreducible representations of the cubic unit cell, that lead to the tetragonal α 1 and α 2 , and orthorhombic β phases.

FIG. 4 :
FIG. 4: Idealized schematic of ferroelastic domain walls in WO 3 .The 2D coordinate system corresponds to the pseudocubic directions.Red and blue lines represent β domain walls with (110) and γ domain walls with (100) crystallographic plane orientation, respectively.Domain sizes and domain wall widths are not representative and β domain walls can also be present with 110 plane orientation.

FIG. 5 :
FIG. 5: The initial unrelaxed supercell for the β domain wall along a crystallographic (110) plane (shown with the black dashed lines).We label the direction perpendicular to the plane of the wall as z.The orthorhombic β unit cell is indicated by the blue box and the ions that were kept fixed during the relaxation by red boxes (with black lines depicting their center).The cubic phonon modes and directions that are already present in the α 2 phase are denoted by black vectors and those introduced by the α 2 − β transition by green vectors.

FIG. 6 :
FIG. 6: Profiles of (a) the angle φ (solid lines) and amplitude |Q| (dotted lines) of the mode Q, (b) the charge density ρ (in milli-electronic charges per cubic angstrom) and (c) the total amplitude |R − 5 | of the cubic R − 5 mode in (c) across domain walls, calculated with the Landau-Ginzburg model.The total amount of additional charge corresponds to the doping level x (in electrons per formula unit) depicted in the color bar.

FIG. 7 :
FIG. 7: (a) Rotation angle φ (blue) and amplitude |Q| (red) of the Q mode along the [110] direction (z) of the relaxed β domain wall supercell.(b) Change in macroscopic planar charge density ∆ρ (in milli-electronic charges per length(in Å)) upon addition of charge in the relaxed β domain wall structure.The red and green lines represent the upper and lower envelopes of the charge densities of red boxes as depicted in Fig. 5.The dashed and solid lines in both panels indicate the positions of the fixed bulk unit cells and domain wall centers as in Fig. 5.

FIG. 8 : 2 F
FIG. 8: Calculated α 1 electron bands (a), phonon bands and phonon linewidths (b), phonon density of states and α 2 F (c) for the case of x = 0.125 additional electrons per formula unit.The phonon linewidths are shown as vertical red bars in the phonon band structure plot and they are scaled by a factor of 10 for visibility.Occupied states in the electron band structure are shown in red, unoccupied states in blue, and the horizontal line at zero eV is the Fermi level.

FIG. 9 :
FIG.9: Calculated phonon bands (left) and α 2 F (right)for the α 1 structure, for a range of added electron concentrations.Blue corresponds to the lowest (0.125 electrons per formula unit) and orange to the highest (0.24 electrons per formula unit) doping levels, with successive lines corresponding to successive points in the plots of Fig.10.

FIG. 10 :
FIG. 10: Calculated superconducting critical temperature T c and electron density of states at the Fermi level n(E F ) in the α 1 phase, as a function of added electrons per formula unit, x.

5 FIG. 13 :
FIG. 13: Initial unrelaxed supercell for the γ domain wall corresponding to a crystallographic (100) plane.The monoclinic γ unit cell is indicated by the blue box and the ions that were kept fixed during the relaxation by red boxes.The black vectors indicate the cubic phonon modes and directions that are already present in the β phase; those introduced by the β − γ transition are shown by green vectors.The orthorhombic Γ + 2 mode causes a slight monoclinic tilt in one of the directions of the cubic M − 2 and R − 5 modes each of which is illustrated by circular arrows.

FIG. 15 : 2 F
FIG. 15: Calculated electronic bands (a), phonon bands and phonon linewidths (b), phonon density of states and α 2 F (c) for x = 0.125 additional electrons per f.u. in the α 2 phase.The phonon linewidths are shown as vertical red bars in the phonon plots and they are scaled by a factor of 10 for visibility.Occupied states in the electron band structure are shown in red, unoccupied states in blue.Horizontal line at zero denotes the Fermi level.

TABLE I :
Lattice parameters (in Å and degrees) of relaxed structures obtained in this work with LDA-VASP (upper rows) and LDA-Quantum Espresso (lower rows) compared with experimental values from the literature.

TABLE II :
36,66 amplitudes of the main cubic phonon modes (in Å) in relaxed structures obtained with LDA-VASP (upper rows) and LDA-Quantum Espresso (lower rows) compared with experimental values from the literature.Phonon mode amplitudes were obtained with ISODISTORT36,66.The amplitudes correspond to the summed atomic displacements normalized to the cubic cell relaxed with LDA-VASP.
. As mentioned above, our calculated phonon linewidths FIG.11: Calculated superconducting temperature T c and total coupling strength λ as a function of added electrons per formula unit, x, in the α 1 phase.

TABLE III :
Landau-Ginzburg parameters and their dependence on doping calculated to sixth order using DFT for a doping range of -0.25 to 0 me Å−3 .