Electron Density Response to Phonon Dynamics in MgB 2 : An Indicator of Superconducting Properties

Electron density differences resulting from atom displacement patterns aligned with phonon modes in MgB2 have been calculated using density functional theory (DFT). The extent of phonon anomalies, identified as indicators of the superconducting transition temperature, Tc, under a range of conditions in AlB2-type structures, reduce as boron atoms are displaced from their equilibrium positions along E2g mode directions. The Fermi energy for displacements along the directions of the E2g phonon mode accounts for changes in the covalent B−B bond electronic charge density. We applied differential atom displacements to show that the shifted σ band structure associated with the light effective mass became tangential to the Fermi level and that the Fermi surface undergoes a topological transition at a critical relative displacement of ~0.6% of the boron atoms from equilibrium. The difference in Fermi energies at this critical displacement and at the equilibrium position correspond to the superconducting energy gap. The net volume between tubular σ surfaces in reciprocal space correlated with the depth of the phonon anomaly and, by inference, it is a key to an understanding of superconductivity. This ab initio approach offers a phenomenological understanding of the factors that determine Tc based on knowledge of the crystal structure.


Introduction
Density Functional Theory (DFT) offers an effective alternative to solution of a the ground state electronic properties of a crystalline material [1].For many materials including superconductors, DFT provides facile resolution of fine-scale structural variations (i.e.picometer) because electron density is the key variable in real, three-dimensional coordinate space for a particular crystallography [1].This fundamental attribute of DFT arises because charge distribution calculations are an explicit feature of these methods [2].It follows then, that detailed analysis of calculated charge distributions could provide insight into the underlying mechanisms of superconductivity and other electron dependent phenomena.We evaluate charge distributions in the boron-boron (B−B) plane of a conventional superconductor, MgB 2 , that complement earlier work on phonon distributions in the type structure.There are many theoretical and experimental studies on MgB 2 extant in the literature.This abundance of information combined with a simple structure containing two atom types provides an excellent reference for detailed DFT modelling of electron distributions in a conventional BCS superconductor.
The superconducting transition temperature, T c , not only determines the practical applicability of a superconducting material, but it also helps identify mechanisms from which superconductivity emerges [3].Current approaches to determine the T c of a material from knowledge of its crystal structure still rely on adjustment of Coulomb pseudo-potentials and effective electron phonon coupling constants; thus, significant ab initio value and prediction capacity are lost.In light of excellent correlations found between temperatures derived from anomalies in phonon dispersions and T c for AlB 2 type compounds [4] [5], under a wide range of external conditions and without any corrections or adjusting parameters [6] [7], this work explores electron density responses to the phonon dynamics.This ab initio approach provides a basis for design of new materials with a priori predictable superconducting transition temperatures.
A valuable contribution to superconductivity models is the use of DFT to calculate the deformation potentials [8] of MgB 2 , a superconductor at 39 K [9].This approach determines the amount of energy that breaks the degeneracy of σ electronic bands per unit length of B−B inter-bond distances.These distances are modified according to specific vibration patterns and then frozen in a model calculation [2].Using this approach, An and Pickett [2] observe that the displacement patterns of B 1g , A 2u , and E 1u modes produce no visible effect in the electronic bands of MgB 2 .They estimate that the deformation potentials of these modes are at least a factor of 25 smaller than those for displacement patterns associated with the E 2g mode [2].This approach clearly identifies the E 2g vibration mode as a key factor in the superconducting mechanism for MgB 2 as emphasised by many researchers [10] [11] [12].However, the variation of these modes, or of the total population of these modes, with temperature and the interplay between inter-bond distances with temperature are less evident even though these features may be directly relevant to superconductivity mechanisms.In this work, Modeling and Numerical Simulation of Material Science we explore these variations of modes with inter-bond distance which we show corresponds to effects of temperature.
Recently, deformation potentials have been used in McMillan equations to explain the effects of structural parameters on the T c of MgB 2 and substituted forms [7] [8].The McMillan equation determines reasonably well the T c of known superconductors based on normal state parameters [8] [10] [13].This fact offers hope that standard DFT may provide all the necessary characteristics of a superconducting material based only on knowledge of the crystal structure [14].However, the temperature dependencies of T c estimates using McMillan equations are often not identified in preference to comparisons of the electron density of states for different related structures [7] [8].For example, recent attempts to identify thermal effects on superconductivity from lattice vibrations invoke zero point lattice fluctuations [15] without consideration of increased amplitude of vibration for temperatures above 0 K.We now extend these concepts and, in taking a lead from the early work of An and Pickett [2], use DFT calculations to examine the deformation of specific atomic bonds or their charge distribution relevant to particular vibration modes within a structure.We also make an explicit correlation of the increased amplitude of vibration expected as a function of temperature for different phonon modes.In this work, we show that such an approach using ab initio DFT methods is an effective way to visualise the transition to superconductivity in MgB 2 .
In addition, we use this approach to re-examine the effects of displaced atom positions that align with key phonon vibration modes on the electronic properties of MgB 2 .We utilize this technique to identify specific response(s) of electron density to these atom displacements and relate these effects to changes in the Fermi energy and Fermi surface of MgB 2 .These changes, explicitly determined using DFT calculations, explain variations of the depth of phonon dispersion (PD) anomalies with strong correlation to T c [5] [6], directly influence superconductivity and establish a direct connection to variations in temperature.
While these correlations are exemplified for MgB 2 in this article, we submit that this phenomenological approach is applicable across a broad range of compounds.

Computational Methods
DFT calculations are undertaken using the CASTEP module of Materials Studio 8.0, using the linear response method from this software suite [16]  Unless specified, calculations have been customized to spin-unpolarized metal, cut-off energy 990 eV, a non-relativistic treatment and convergence criteria as follows: energy at 5 × 10 −6 eV per atom; maximum force at 0.01 eV Å −1 ; maximum stress at 0.02 GPa and maximum atom displacement at 5 × 10 −4 Å.
For unit cell calculations on MgB 2 with P6/mmm symmetry, lattice parameters a = b = 3.085A and c = 3.523A are used as input values [9] [18] [19].For electronic band structure calculations, the computation time is significantly shorter than that required for PD calculations [5].Therefore, higher resolution grids (i.e. for k-grid separation Δk < 0.02 A −1 ) are utilized for band structure calculations.We also define relative atom displacement values, D x , in terms of relative unit cell coordinates for the B−B plane.Thus, a value for D x = 0.008 defines a movement, or displacement of the boron atoms, along a real space direction relative to the a lattice parameter of 0.8% relative where the cell parameter is 1.0.Using the unit cell dimensions for MgB 2 , this example represents a displacement of 2.468 pm.A positive or negative sign for the D x value designates an extension or compression of the atom-atom displacement, respectively.This atom displacement is annotated with respect to the B-B plane due to its significant contribution to superconductivity in MgB 2 [2] [11] [13].
Materials Studio 8.0 allows for the determination of electron density differences and electron localization functions by selecting appropriate options in the set-up for each calculation [17].Materials Studio provides the option of generating the electron density difference in two different ways: 1) to calculate electron density differences with respect to a linear combination of the atomic densities or 2) to calculate the electron density differences with respect to a linear combination of the densities of sets of atoms contained in the model.Alternatively, the program provides for both approaches to be selected simultaneously.
For this work, electron density differences are calculated with respect to a linear combination of the atomic densities.
The electron localization function (ELF) provides a useful qualitative description of chemical bonding.This function, retrievable from Materials Studio 8.0 through set-up conditions, follows the approach introduced by Becke and Edgecombe [20] and the revised interpretation in more recent literature [21].The outputs from electron density difference and electron localization function calculations are very sensitive to choice of the iso-value for the corresponding iso-surfaces.Within a comparison series, iso-values are selected in order to enhance the contrast of properties but are maintained at constant value within a series of computations.This approach highlights the influence of the displacement parameter and changes made to this parameter in a computational series.Electron density differences are determined for symmetric arrangements of atoms and for corresponding modified positions that are displaced along orientations aligned with particular phonon modes.This computational approach is similar to calculations used by An and Pickett [2] to determine deformation potentials.Deformation potentials in conjunction with McMillan's formulation Modeling and Numerical Simulation of Material Science [22], are shown to successfully estimate T c for MgB 2 [8].This approach is visually demonstrated on two web sites [23] [24] which display the dynamic behaviour of σ electronic bands that are strongly linked to E 2g phonon modes in MgB 2 .
For the initial structural model, geometry optimization is undertaken as reported previously [4] [5].However, for subsequent calculations in which an atom displacement position is assumed, further geometry optimization is not undertaken.In this way, a frozen picture of the electron distribution at an intermediate position of the particular phonon mode is created.We have evaluated displacement patterns corresponding to the A 2u , E 1u and both E 2g modes, as described by Kunc et al. [10].For the E 2g modes, these displacement patterns predominantly involve boron atoms in the a-b plane.
With displaced atom positions frozen in place, electronic band structures are calculated using the selected configuration with appropriate consideration of structural symmetries invoked by the displaced atom positions.By constraining these atom positions to a new symmetry condition, the computation is simplified and the calculation time is reduced.We have demonstrated in earlier work [5] [6] [7] that comparison of models using the LDA and GGA functionals provides an excellent measure of internal consistency for AlB 2 -type structures and a source of error estimates for a generic model.
In order to calculate the phonon dispersions (PD's) for displaced positions of boron atoms, geometry optimization is required, and may be limited to atom projections without changing lattice parameters.To prevent geometry optimization from relocating displaced atoms to their equilibrium positions, a primitive symmetry and very high force convergence tolerances are used.For the MgB 2 structure, force tolerances of 0.35 eV/Å and 0.5 eV/Å are required for bond distance changes up to ±0.006 (= ±0.6%) and for ±0.008 (= ±0.8%) in terms of relative coordinates.Thus, large restoring forces are maintained and accordingly, computational models achieve convergence for PD calculations.Because these PD calculations are computationally intense (using 196 q-points), we have limited this type of study to one E 2g PD mode and the corresponding σ band Fermi surface using the LDA functional.Outcomes from this calculation are shown in Section 3.3 for selected displaced atom positions with reduced symmetry approximating the G−A, G−K and G−M directions.

Modelling Results
We present modelling outcomes in the following sections that analyse the impact of atom displacement on phonon modes and on band structure.Section 3.1 identifies the phonon modes that produce electron density differences as a result of an induced "frozen" displacement pattern.The modes that induce electron density differences are described in detail in Section 3.2 through effects on the electronic bands with increased atom displacement from an equilibrium position.We briefly describe the influence of these atom deformations on Fermi energies and their surfaces.Finally, Section 3.3 confirms that at an appropriate Modeling and Numerical Simulation of Material Science displacement from equilibrium, the phonon anomalies collapse with clear implications for the Fermi surfaces and, by inference, for superconductivity of MgB 2 .

Electron Density
As noted, the choice of iso-value greatly influences the visual representation of electron density differences and the electron localization function.Furthermore, to illustrate the influence of atom displacement, we show in Figure 1 the difference in electron density (in blue) for the symmetric MgB 2 structure and for two situations with perpendicular in-plane displacements of boron atoms.These in-plane displacement directions correspond to two different E 2g modes in the MgB 2 structure, with an exaggerated relative displacement of +0.05 (in relative coordinates) using the same iso-value of 0.15.This relative displacement value may be extreme compared to what is physically reasonable.However, the displacement demonstrates that substantial changes to the electron density occur within the a-b plane of the MgB 2 structure when boron atoms are shifted along the dominant optical phonon mode direction.
For each calculation represented in Figure 1, the iso-value is kept constant in order to avoid introducing secondary effects on these results.Figure 1    shows the in-plane symmetry of MgB 2 is similar to that in Figure 1(a).
Figure 3 shows the electron density difference of the E 1u mode for a similar level of relative displacement with (+0.05 in relative coordinates) along the direction of the B−B bond. Figure 3 shows atom positions viewed along the c axis in which the Mg atoms are shifted in the direction of the arrows or approximately along <110>.In Figure 3, Mg atoms are closer to the B atoms of the adjacent hexagonal ring rather than symmetrically at the centre of the (001) projection in an undistorted unit cell.Nevertheless, a change in symmetry of the electron density difference on the boron plane is not detected with this imposed relative displacement of atoms along E 1u modes.
To illustrate the effect of an incremental increase or decrease in B-B bond length due to incremental displacement of atoms along an E 2g phonon mode, Figure 4 shows the resulting electron density differences in a series of single unit cells with superimposed electron density differences (in light blue) and localization functions (in dark blue) using the same iso-value (e.g.0.15) as in Figure 1.
Scrutiny of the series of relative displacements in Figure 4 shows that there is an asymmetric arrangement of the electron density differences and the electron   The data in Figure 4 show an asymmetry between the elongated and the compressed directions of the inter-bond distance.For example, in this configuration a complete B−B bond is located in the centre of the unit cell and four half-bonds (i.e. two net bonds) point towards the unit-cell edges.This results in an imbalanced response (by a factor of 2) to compression or extension.In compression, an increased electronic charge is established between the boron atoms within the same unit cell.When these bonds are stretched or extended, adjacent B−B bonds are compressed with a resulting build-up of equivalent increased electronic charge at the edges of the unit cell, but distributed to twice the number of bonds.
It is clear from the above analyses that the E 2g modes are the only modes to generate observable modulated electron density differences within the boron plane.This modelled outcome is consistent with previous experimental data and calculations that show the E 2g mode is strongly electron-phonon coupled and a key enabler of superconductivity in MgB 2 .These electron density differences produce bond charge modulations which have super-lattice periodicities with respect to the original MgB 2 crystal.These changes reflect the effects of mobile electronic charge response along the bond direction(s) due to the changes in interatomic attractions and/or repulsions from atomic displacements induced by the E 2g phonons (see Discussion).

Band Structures
We have evaluated the influence of atom displacements-using similar parameters as in Section 3.1-on the electronic band structure of MgB 2 .Figure 5 shows the electronic band structures for a series of incremental changes in relative displacement from the equilibrium position that follow the movement of an E 2g mode along the B−B inter-bond direction.These band structures are calculated using the LDA functional with k-grid separation Δk = 0.01 A −1 .
With increased distance between displaced boron atoms, the primary influence on the band structures is removal of band degeneracy along the G−A (or G−Z) reciprocal direction adjacent to the Fermi level (Figure 5; circled).In addition, a corresponding separation of the inverted parabolas (nearly free-hole like) occurs at G. Figure 5(e) also shows that an electronic σ band (the inverted parabolas at the G point near the Fermi level) becomes nearly tangential to the Fermi level for the relative atom displacement D x = ±0.006.

Fermi Energy
We determine the calculated Fermi energy for each model of MgB 2 with a specific deformation potential using the LDA and the GGA functionals.As is well known, the Fermi level (and energy) represents the most energetic valence electrons at 0 K which, for bands of metallic behaviour, may participate in electron transport.Only these electrons with energies within a range of k B T, where k B is Modeling and Numerical Simulation of Material Science Boltzmann's constant and T is temperature, participate in electronic transport properties at the Fermi energy.Use of DFT with a consistent approach and software suite can provide accurate values of the Fermi energy for MgB 2 and similar structures with reference to the total number of valence electrons.Note that some earlier computational studies choose to denote a reference point for the Fermi energy at different values to that used in this and our earlier studies lengths and Fermi energy determined with the LDA functional for a range of atom displacements and for a lattice expansion of 1.006x in Table 1.This factor of 1.006x is an approximation to an isotropic expansion of the entire lattice (or unit cell) equivalent to ~40 K.In contrast, an expansion factor of 1.000x represents lattice parameters optimised for absolute zero temperature.
Figure 6 shows that the variation of calculated Fermi energies using the LDA and GGA approximations are similar in form and appearance but with a difference in energy of ~0.38 eV.In all other respects, the curves for both calculations are identical.In both cases, the Fermi energy displays local minima at relative displacements of about D x = ±0.006and a local maxima at the equilibrium position (D x = 0.000).The depth of the energy minimum at the elongated position (D x = +0.006) is approximately two times the value of the energy minimum at the contracted position (D x = −0.006).The difference in Fermi energies between the elongated minimum and the equilibrium position is ~17.5 meV and 15.9 meV, for lattice expansion factors 1.000x and 1.006x, respectively.These values, particularly the difference for the lattice expanded by 1.006x, are comparable to two times the value of the MgB 2 superconducting gap [25] [26].The differences in Fermi energies for the contracted and elongated displacements appear related to the number of contracted bonds identified in Figure 4.

Fermi Surface
The widely accepted topology of the MgB 2 Fermi surface consists of coaxial, nearly parallel "tubular" sections for σ bands in the c axis direction and three dimensional interconnected bands [27].These tubular sections resemble an hourglass shape and are also referred to as warped cylinders [28].Fermi surfaces provide a three dimensional, iso-value representation in reciprocal space of the Fermi energy.We have calculated the changes in Fermi surface for MgB 2 with Modeling and Numerical Simulation of Material Science incremental increase in atom displacement from the equilibrium position.These Fermi surfaces are shown for calculations using the LDA functional with Δk = 0.01 A −1 in Figure 7.
At small atom displacements, the Fermi surface shows a familiar coaxial warped tubular section consistent with previous calculations Parallel warped tubes coaxial with the z-direction shown in Figure 7(a) correspond to the heavy (outer tube) and light (inner tube) effective mass σ bands.The light effective mass σ band displays a topological transition to separate paraboloid sections facing in opposite directions (Figure 7(d)) as the B-B distance is gradually increased (or decreased) about 0.6% in relative coordinates.Other Fermi surfaces near the surrounding Brillouin zone boundaries do not change significantly and correspond to the π bands.

Phonon Dispersions
We have utilized the phonon dispersion (PD) in previous computational analyses to illustrate the strong link between optical phonon behaviour and band structure in MgB 2 [5].We have also shown that the depth of the PD anomaly provides a temperature which corresponds to the superconducting transition temperature for a wide range of external conditions, accurate within experimental and calculation errors [5] [6] [7] [14].Consistent with this approach, we also highlight in Figure 8 the influence of incremental atom displacements on specific phonon modes.As noted in earlier work, PD calculations are computationally intensive [5] [7].Hence, this evaluation of atom displacements is performed using the LDA functional with fewer sampling points and limited reciprocal space directions (e.g.along G−M, G−K and G−A only).Figure 8 shows the changes in PD modes for G−M and G−K reciprocal directions with particular emphasis on the E 2g phonon with increase in relative atom displacement from the equilibrium position.
Figure 8 shows significant change in the form of the E 2g phonon anomaly with increased atom displacement from equilibrium and an increase in energy of one of the degenerate modes at G. At large atom displacement (e.g.D x = 0.008; Figure 8(d)), the E 2g phonon anomaly (below the B 2g mode) is absent and the order of modes (i.e.relative energy or frequency) is reversed between the B 2g and E 2g phonons.Evaluation of the band structures in Figure 5 shows a corresponding change in form of a σ electronic band branch which shifts below the Fermi level.This feature also corresponds with changes in the Fermi surface as shown in Figure 7 where, at relative displacement D x = 0.008, the surface becomes two separate paraboloids.These computational outcomes further confirm the strong link between the electronic bands, Fermi surface nesting and the depth of the PD anomaly.These inter-relationships are clearly associated with the superconducting transition temperature for AlB 2 -type structures [4]

Discussion
In earlier publications, we utilize ab initio DFT to recognize the shape and extent of the E 2g phonon anomaly in MgB 2 [4] [5] as well as in metal-substituted MgB 2 [5] [7] and MgB 2 under the influence of hydrostatic pressure [6].In each of these studies, we compare computational outcomes with experimental data to show that the phonon dispersion (PD) is an effective tool for recognition of a phonon anomaly and that the thermal energy of this anomaly is correlated with the experimentally determined value of T c for the superconducting transition of MgB 2 compounds.This approach-based primarily on calculated PDs-is extended to other AlB 2 -type compositions [5] and is also the basis for prediction of new compounds that are likely to be superconductors [5] [7].
As is well-known, a strong crystallographic link exists between band structure calculations and PD models based on DFT.This link between phonon behaviour and electronic band structure requires fine k-grid resolution and thus high computational cost to elucidate PD data [5] [7].However, MgB 2 offers an ideal structural format to highlight this link between electron and phonon behaviour consistent with the well understood view of electron-phonon coupling as key to conventional superconductivity [30].A connection between the E 2g phonon and the Fermi surface as well as components of the E 2g phonon anomaly with the band structure at or near G is described for MgB 2 [5] [6] and more recently, elucidated by experiment [31] using Inelastic X-ray Scattering (IXS).In this work, we elaborate further on the interaction between the electronic band structure and phonon behaviour of MgB 2 , in particular, as revealed by electron density difference distributions in response to deformation potentials.

Fermi Energy
Born and Oppenheimer introduced the adiabatic approximation, in which electrons are expected to move much faster than heavier nuclei, allowing for a very rapid redistribution of electronic charge, every time a nuclei undergoes movement [32] [33].This assumption allows for substantial separation of nuclear and electronic coordinates while solving the Schrodinger equation [33].As pointed out by Boeri et al. [34], the Born-Oppenheimer approximation is implicitly assumed in DFT, whereby the electron density distribution can be calculated for each ionic configuration.According to Uchiyama et al. [35] the close agreement between DFT calculated electronic bands and Angle Resolved Photoemission (ARPES) data up to binding energies as high as 2.5 eV suggests that electronic correlations in MgB 2 are unimportant.As a result, they contend that this material is well described by conventional band theory [35].
Boeri et al. [34] argue against the validity of the adiabatic approximation for MgB 2 , based on an abnormally small (~0.5 eV) Fermi energy (E F ) calculated for the σ bands, assumed to be of the same order as the highest phonon energies.However, the Fermi energy assigned by Boeri et al. [34] appears related to the distance to the vertex of the parabolic bands and does not correspond to valence Modeling and Numerical Simulation of Material Science electron densities or charge carriers determined experimentally from resistivity measurements and from Drude approximations.For example, measurements of the Hall coefficient for MgB 2 , indicate that the charge carriers are holes with a density of 1.7 -2.8 × 10 23 holes cm −3 at 300 K [27] [36] [37].In contrast, calculations by Suzuki [38], referred to the vertex of the parabolic bands, estimate a carrier concentration of ~10 22 holes cm −3 an order of magnitude different to experiment.
The disparity between calculated Fermi energy and experimental measurements-albeit indirect-appears in many publications that refer the energies to the top of the valence bands and reset the reference Fermi level to zero energy.
We consider this disconnect between the absolute Fermi energy values calculated by DFT (from the total number of available valence electrons and states) and the traditional "zeroing" of the reference level to cause loss of valuable insight and to be a source of potential confusion in the literature.A critical factor for this work, and others requiring similar detail and accuracy, is that the topologies of Fermi surfaces are determined by the absolute calculated Fermi energy with respect to Brillouin zone dimensions, where valence states are essentially filled with the total number of valence electrons (with a two times division to account for spin directions in spin unpolarized calculations).This change in reference point for the Fermi energy, combined with aggregation of information in reduced Brillouin zone schema in reciprocal space lead to complications with interpretation as outlined in the next section.
The calculated Fermi energies for a range of configurations in MgB 2 , presented in Table 1 are consistent with earlier DFT models for metal substituted or pressure dependent forms [5] [6] [7].DFT calculations for MgB 2 using CASTEP in Materials Studio result in a Fermi energy of 8.3998 eV and 8.1165 eV using the LDA and GGA approximations, respectively.The results from CASTEP are more consistent with the high carrier densities determined from Hall effect measurements [27] [36] [37].Thus, we consider that the absolute value of Fermi energy ranges between 8.1 eV and 8.4 eV as consistently demonstrated in a number of DFT computations undertaken previously [5] [6].The insight we learn from ab initio DFT is that choice of appropriate functionals and proper account of valence electron numbers help elucidate the strong correlation between electronic, vibrational and superconducting properties for MgB 2 .Furthermore, we argue that the adiabatic approximation is a justifiable assessment for MgB 2 for a wider range of conditions than previously considered.

Fermi Surfaces, Phonons and Gap Values
The Fermi surface separates the filled electronic states from the empty electronic states in three-dimensional reciprocal space at absolute zero temperature [39].
The widely accepted topology of the MgB 2 Fermi surface consists of coaxial, nearly parallel "tubular" sections for σ bands in the c axis direction [27].The outer and inner tubular sections for MgB 2 correspond to the heavy and light ef-Modeling and Numerical Simulation of Material Science fective masses for the σ bands, respectively [2] [28] [40].The hole nature of the tube implies that at absolute zero temperature the inside of the tube is empty and the outside of the tube is filled.This requires that in an extended zone scheme, tubular Fermi surfaces should be located away from the origin, possibly at the corners or outside the first Brillouin zone.Referral of the Fermi energy to the vertex of the approximate parabolic σ bands may unnecessarily confuse interpretation of electronic parameters and construction of Fermi surface topologies in this structure and others.Apparent confusion in the literature arises because in graphic representations, boundaries for Fermi surfaces are displayed without a clear delineation of where electrons and empty states, respectively, reside.
Furthermore, two parallel tubular Fermi surfaces also create a dilemma when the two tubes are represented in a coaxial fashion.This is because at zero temperature, in the volume between these two surfaces, electrons occupy energy levels outside the inner light effective mass tubular section, while empty states with identical magnitude of momentum and energy exist inside the outer, heavy effective mass tubular section [27] [30].We suggest that this dilemma is, in practical terms, resolved via the schematic shown in Figure 9 for the coaxial σ Fermi surfaces of MgB 2 .In this schematic, we neglect the warping of the σ Fermi surfaces as previously described [6].
In the reduced zone scheme of the Brillouin zone, the inter-tube volume is an important region where levels close to the Fermi energy with identical k-vectors plus or minus a reciprocal space G vector (and similar energy) show filled character relative to the light effective mass, but an empty character relative to the heavy effective mass [6].We assume that a reduced zone representation, in essence, aggregates tubes of different diameter that occur at different places in an extended Brillouin scheme, connected by reciprocal space G vectors.Therefore, we infer that an extended zone scheme explicitly manifests a resonating behaviour between the tubes at different reciprocal space locations (shown in Figure 9) [6] and may be a more appropriate representation of these states compared with a reduced zone schematic.
We have previously shown that the apparent inter-tubular volume controls the depth of the phonon dispersion anomalies in AlB 2 -type structures [5] [7].
We have also demonstrated a correlation between a temperature associated with the PD anomaly (i.e.T δ ) and T c , and, as an extension of this phenomenom, propose that this inter-tubular volume is a key region in reciprocal space for superconductivity.As pointed out by Geballe et al. [3], T c itself and its response to controlled changes in external parameters, often help reveal the responsible mechanisms for superconductivity.Given the excellent correlation between the depth of PD anomalies, Fermi surfaces and T c , it appears that useful insights on superconducting mechanisms are hidden within these inter-tubular Fermi surface constructs.

Dynamic Fermi Surfaces and Bond Charge Modulations
Figure 8 shows a clear transition of Fermi surfaces with increased atom displacement from the equilibrium position for MgB 2 .This transition involves the collapse of one σ Fermi surface as the B−B distances are displaced relative to the equilibrium value.This transition also corresponds to a disruption of the coherent interconnection, or nesting, between the two tubular sections.When these tubes are parallel, the hole/electron states are at exactly the same k-point plus or minus a reciprocal space vector G and similar energy for the separate light and heavy mass surfaces.These filled/vacant electron states can be nested by phonons of the exact same magnitude.However, once one surface collapses, this perfect nesting is disrupted and the model suggests that this is when the superconducting T c is breached.This modelling also suggests that the condition at which these tubular sections cannot be connected by the same magnitude phonon may signal the beginning of incoherent scattering in MgB 2 .
The form of parallel Fermi surfaces shown in Figure 9 (and, more precisely, in Figure 7) are similar to the ideas and figures discussed by Bardeen (Figure 1 in ref [41]).In Bardeen's article, a shell structure around the Fermi surface is identified as the region of electron states that are important for superconductivity.The Fermi surface of MgB 2 for significantly displaced boron atoms is also similar to the effect of Al substitution in Mg 1-x Al x B 2 [42] (compare Figure 7 of this work with Figure 4 of [42]).Both results show an evolution to a topology with two separate cusps, with corresponding degradation of superconducting properties.However, Fermi surfaces described by de la Peña et al. [42] consist of two Modeling and Numerical Simulation of Material Science separate sets of cusps that retain similar parallel shapes.In contrast, and as shown in Figure 7, after the critical displacement or separation of boron atoms is achieved (at relative displacement D x = ±0.006),one surface completely collapses and the remaining Fermi surface cusps consist of single surfaces.Alternatively, we may consider that as the temperature is lowered towards the superconducting transition temperature, the mixed topology of parabolic cusps and tubular sections is interacting more strongly and evolving into the set of approximately parallel tubular sections.This interpretation would then be more in-line with important concepts proposed earlier on topological transitions that determine or enhance T c [14] [15] [43].The E 2g frequency of the appropriate branch increases as the B-B distance is displaced from the equilibrium position.The simultaneous drop in the Fermi energy seems to account for the energy that is passed onto the generated phonons, as a result of an electron jumping to alternative bonds responsible for the modified electron density differences (see Figure 4 and Figure 6).In the PD, this manifests as a higher E 2g frequency with an associated reduced PD anomaly.Once the displacement exceeds the "gap" limit, it will generate such energetic phonons that other modes become involved and the system cannot maintain coherency.A mechanical analogy to this condition is as if a natural spring is overstretched and then released, causing a disturbed movement back to the equilibrium position.
A schematic of the bond charge modulation pattern in the boron plane of MgB 2 is shown in Figure 10.The charge modulation has a complex pattern that only appears as the atoms deviate from their equilibrium positions.As the restoring forces bring atoms back to their equilibrium positions, the modulated charge pattern disappears.The dynamic origin of the bond charge modulation observed in MgB 2 may be different to that of a standard charge density wave (CDW) because of the specific crystallography of the structure.The restoring forces from phonons or lattice vibrations are poised to re-equilibrate geometrical and electronic structures to no modulation.For motion in the reverse direction, a new electron density difference is generated and then disappears in a resonant fashion.In a CDW, charge modulation is expected to be more persistent or permanent, since its manifestation is more prevalent as determined via a wide range of experimental characterization techniques.Nevertheless, this work may shed light on geometrical relationships between nesting vectors (in reciprocal space) and bond charge modulation (in real space) [44].

Temperature Dependent Mechanism
Two of the most direct manifestations of temperature effects on the geometrical structure that can be measured are the thermal expansion and the temperature dependence of lattice parameters.Both thermal expansion and changes in lattice parameters with temperature have been investigated for MgB 2 [45] [46] [47] [48] [49].Thermal expansion is modest or even negative before and soon after the superconducting transition temperature, but develops a much sharper and linear increase at higher temperatures.The change in lattice parameters with temperature follows a similar pattern.Detailed neutron refinement investigations by Jorgensen et al. [45], provide Debye Waller factors and mean square atomic displacements over a wide temperature range for MgB 2 .
Our atom displacement calculations impose an increasing change in the B-B inter-bond distance from the equilibrium distance.This approach indirectly mimics the effect of increased temperature on this inter-bond distance.For example, the mean square displacement for boron at 37 K [45] is measured as U 11 = 0.00333 (11) Å 2 , U 33 = 0.00455 (18) Å 2 and U 12 = 0.00166 (6) Å 2 .The square roots of these values are 0.0577 (10) Å, 0.0674 (13) Å and 0.0407 (7) Å, respec-Modeling and Numerical Simulation of Material Science tively.If we compare these values to the difference in B-B bond distance at equilibrium and at a relative displacement D x = ±0.006at the transition temperature (~39.4K), we obtain 0.0631 Å (= 1.7551 -1.6920) and 0.064 Å (= 1.7782 -1.7142) for the LDA and GGA functionals.This relative displacement is for an expansion factor of 1.000x in all three real space directions.These calculated values display good correspondence with the square root of the mean square displacement for boron at 37 K based on precise experimental data [45].Therefore, our assumption that a relative displacement of D x = ±0.006corresponds approximately to the superconducting transition temperature (~39.4K) is reasonable given the relative errors in calculations and in experimental measurements.The estimated inter-bond displacement required for collapse of superconductivity is consistent with experimental measurements of the square root of the mean square displacements by Jorgensen et al. [45].This consistency between theory and experiment further supports our suggestion that these calculated effects are directly linked to temperature.
A critical, and perhaps unique, aspect of superconductivity in MgB 2 is that electron density differences and phonon modes are largely confined to the a-b plane.This configuration allows a direct correlation of DFT calculations to key properties in the structure.This computational approach highlights that the electron density difference responds to increased atom-atom repulsion as the distance between atoms is reduced.As a result, an increased electronic bond charge is required to minimize the overall energy.These bond charges are localized along the highly covalent inter B−B bond distances, and along with a phonon mode that is oriented along the same direction, seem to moderate the correlation.The proportion of transferred charge density directly affects the Coulomb repulsions and can translate into a net force change.This net force is related to a change in the magnitude of the phonon frequency.
We note that other, more complex structures may not show phonon modes with movement aligned to a specific inter-bond distance.In such a case, more than one phonon may influence electron density differences for a particular bond direction with resulting complex patterns of electron density differences and phonon displacements.In addition, magnetic interactions may also add to the force constants with consequent influence on interatomic distances and displacement patterns.Not all of the structures will exhibit a situation equivalent to the important interplay of approximately parallel σ bands and maintain coherency as electron density differences minimize the energy in response to the atom movements along specific phonon modes.However, the existence of parallel tubular Fermi surfaces, coaxial in a reduced Brillouin zone scheme, is ubiquitous to all layered superconductors.This phenomenom may hide important information about superconductivity which may require disaggregation to improve our understanding for effective materials design.Given the accurate validation of these calculations for a range of external conditions, we propose that this approach may be valid for many other examples of superconductivity.

Conclusions
The Fermi energy of the valence electrons in MgB 2 is sufficient to justify the adiabatic or Born-Oppenheimer approximation in quantum mechanics calculations for MgB 2 .In this work, we present a model that identifies the distribution of electrons and empty states on the MgB 2 Fermi surface and, through this, a mechanism associated with the inter-tubular volume between the σ Fermi surfaces that tracks the onset (or demise) of superconductivity.
The net volume enclosed between the parallel heavy and light effective mass σ Fermi surfaces when represented in the reduced Brillouin zone is an important region of reciprocal space for superconductivity.This region, when represented in an extended zone (e.g. Figure 10; [6]) is suitable for a resonant interaction between separate heavy and light hole Fermi surfaces.This inter-surface region is also ideal for nesting to occur between different sections of the Fermi surface, by a phonon vector of constant magnitude, within a thermal energy range k B T from the Fermi energy.Analysis of frozen, displaced boron atoms along the B−B inter-bond distance confirms that the phonon anomaly collapses as the inter-bond distance increases.The variation of Fermi energy for displacements along the E 2g phonon mode that stretches or contracts along the inter-bond distance is smooth and accounts for changes in the electronic charge density along predominantly covalent B−B bonds, which reduces the overall energy by neutralizing increased atomic core repulsions.
The difference between the Fermi energy at a critical displacement of ~0.6% (D x = ±0.006) in relative coordinates and the Fermi energy at the equilibrium position displays a strong correlation to the superconducting energy gap, or multiples of the gap, depending on the number of bonds that are involved in electron density transfer.Our analysis of the shift in E 2g and B 2g PD bands with atom displacement suggests that the interplay between potential energy and kinetic (or thermal) energy stored in the inter-bond distance is a key determinant of T c in AlB 2 -type materials.This analysis of displaced atomic positions also highlights that the Fermi surface is time dependent and dynamic in character.
The dynamic character of the Fermi surface is not only important for an understanding of superconductivity, but also for many other electronic properties.
Earlier studies using McMillan equations, invoked normal state parameters to explain or predict T c , and, consistent with this, DFT is a tool that provides all necessary information to understand superconductivity in AlB 2 -type structures.
Phonon dispersion (PD) calculations shown in this work linked to band structure calculations via DFT, is an effective, reliable and consistent method to accurately determine a temperature strongly correlated with the superconducting T c of an AlB 2 -type material from knowledge of the crystal structure.
How to cite this paper: Alarco, J.A., Talbot, P.C. and Mackinnon, I.D.R. (2018) Electron Density Response to Phonon Dynamics in MgB 2 : An Indicator of Superconducting Properties.Modeling and Numerical Simulation of Material Science, 8, 21-46.https://doi.org/10.4236/mnsms.2018.82002Modeling and Numerical Simulation of Material Science many-body Schrodinger wave equation for multi-atom structures to determine [17].The linear response within the Local Density Approximation (LDA) and Generalized Gradient Approximation (GGA) with norm-conserving pseudo-potentials, a plane-wave basis set and a dense k-grid provide the most consistent model outcomes for the MgB 2 structure [4] [5] [6].Calculations are undertaken with a High-Performance Computing (HPC) facility, using numbers of cores which are multiples of the k grids defined in the reciprocal a* and b* directions of the MgB 2 Brillouin zone.The number of HPC cores is typically between 60 and 80. Modeling and Numerical Simulation of Material Science demonstrates that electron density migrates from regions with longer B-B distances to regions with shorter B−B distances.This outcome is not unexpected because an increase in B−B repulsion requires a corresponding response in electronic bond charge population to reduce the overall energy.However, the sustained localization of the electron density differences along the highly covalent B−B bonds is informative.Similar visualizations of electron density differences can be obtained for smaller atom displacements with appropriate adjustment of the iso-value.

Figure 1 .
Figure 1.Calculated representations of electron density differences (in blue) for MgB 2 viewed along a direction equivalent to the c-axis: (a) the symmetric structure without atom displacements; (b) and (c), with two perpendicular in-plane displacements of boron atoms (red dashed arrows show direction) corresponding to the E 2g modes.Mg atoms are in green, boron atoms are pink.Note the substantial modulation of electron density aligned with planar boron atoms for displacements along E 2g directions.

Figure 2 (
b) shows that the boron planes, with the superimposed electron density difference, are closer to Mg atoms in the plane immediately above the boron plane than to the Mg plane immediately below.This outcome is also expected for the A 2u mode of the MgB 2 structure.The view of electron density difference along the c-axis (Figure 2(a))

Figure 2 .
Figure 2. Electron density difference distributions (in blue) in MgB 2 for atom displacements corresponding to the A 2u mode; (a) viewed along the c-axis direction and (b) viewed along the a-axis direction perpendicular to (a).Red dashed arrows indicate direction of atom displacement(s).

Figure 3 .
Figure 3.View along the c-axis direction of MgB 2 showing the electron density difference (in blue) for displacements corresponding to the E 1u mode (red dashed arrows denote direction of displacement).Note the minor shift of Mg atoms in the direction of the displacement; light blue grid lines are guides only.

Figure 4 .
Figure 4. View along the nominal c-axis direction of MgB 2 showing the electron density difference (light blue) with superimposed electron localization functions (dark blue) for a series of incremental atom displacements corresponding to one of the E 2g modes.The relative displacement values (D x = ±0.00n) are enumerated for each unit cell configuration.Mg atoms are green and boron atoms are pink; red arrows are guides to regions of significant change.

Figure 5 .
Figure 5. Electronic band structures for MgB 2 calculated using the LDA functional with B-B atom displacements along an E 2g direction corresponding to (a) D x = 0.0; (b) D x = −0.002;(c) D x = −0.003;(d) D x = −0.004;(e) D x = −0.006and (f) D x = −0.008.Note that for a displacement around −0.006, the lower parabolic branch of the electronic band at G becomes tangential to the Fermi level, and for −0.008 this branch is already below the Fermi level.

[ 4 ]
Figure 6.(a) Fermi energies using the LDA and GGA functionals with Δk = 0.01 A −1 (and lattice expansion factor 1.000x) for MgB 2 with various displaced B-B bond distances; (b) Enlarged view of the Fermi energies in the region within the red dotted rectangle.These data are obtained using the LDA functional calculated at a k-grid, Δk = 0.01 A −1 .

[ 2 ]
[28].With an increase in the relative displacement of atom(s), the internal warped tubular cylinder shows a reduced radius which is exaggerated at the reciprocal lattice origin, G.With further increase in atom displacement, this feature manifests as "necking" in the centre of the cylinder at G as shown in Figures7(b)-(c).Ulti-Modeling and Numerical Simulation of Material Science mately, the Fermi surface collapses into a thin line at the mid-point of the reciprocal cell when the relative atom displacement, D x = −0.006as shown in Figure7(c).This relative atom displacement is the point at which the corresponding σ band becomes almost tangential to the Fermi level as identified in Figure5.With further atom displacement, the σ Fermi surface breaks apart into two approximate paraboloid sections that face in opposite directions as shown in Figure 7(d).

Figure 8 .
Figure 8. Phonon dispersion plots for MgB 2 calculated with the LDA functional for Δk = 0.01 A −1 and primitive symmetry for various displaced B−B bond distances: (a) D x = 0.000; (b) D x = 0.004; (c) D x = 0.006 and (d) D x = 0.008.E 2g vibration modes are in red; B 2g modes are in blue.

Figure 9 .
Figure 9. Idealised schematic of the σ Fermi surfaces for MgB 2 represented in an extended Brillouin zone schema viewed (a) down the c-axis and (b) perspective view at an angle to the c-axis; (c) and (d) show equivalent resonant states for (a) and (b) relative to the reciprocal space vector G. Warping of the σ Fermi surfaces is neglected in this schematic.

Figures 4 -
Figures 4-8 also highlight the duality of real and reciprocal space interactions between Fermi energy, phonon frequencies and electron density.This duality is because of a direct connection to the directional bond charge distribution that modulates the repulsive B−B forces determining the key E 2g phonon frequencies and the apparent strong localization of the electron density at the Fermi energy along the highly covalent B−B bonds.Similar undulating shapes of the PD curves [4], and the dependence of the Fermi energy on the relative atom displacement shown in Figure 6 reflects this duality.The E 2g frequency of the appropriate branch increases as the B-B distance is displaced from the equilibrium position.The simultaneous drop in the Fermi energy seems to account for the energy that is passed onto the generated phonons, as a result of an electron jumping to alternative bonds responsible for the modified electron density differences (see Figure4and Figure6).In the PD, this manifests as a higher E 2g frequency with an associated reduced PD anomaly.Once the displacement exceeds the "gap" limit, it will generate such energetic phonons that other modes become involved and the system cannot maintain coherency.A mechanical analogy to this condition is as if a natural spring is overstretched and then released, causing a disturbed movement back to the equilibrium position.A schematic of the bond charge modulation pattern in the boron plane of MgB 2 is shown in Figure10.The charge modulation has a complex pattern that only appears as the atoms deviate from their equilibrium positions.As the restoring forces bring atoms back to their equilibrium positions, the modulated charge pattern disappears.The dynamic origin of the bond charge modulation observed in MgB 2 may be different to that of a standard charge density wave (CDW) because of the specific crystallography of the structure.The restoring forces from phonons or lattice vibrations are poised to re-equilibrate geometrical and electronic structures to no modulation.For motion in the reverse direction, a new electron density difference is generated and then disappears in a resonant fashion.In a CDW, charge modulation is expected to be more persistent or

Figure 10 .
Figure 10.Schematic of the bond charge modulation pattern in the boron plane of MgB 2 based on DFT calculations for displaced atom configurations.Blue spheres represent boron atoms; dark blue bonds represent regions of higher electron density and light blue bonds represent lower electron density.Distortion along the x direction (equivalent to an E 2g mode direction) for the average hexagonal unit cell is accentuated in this schematic.

Table 1 .
Parameters for MgB 2 calculated using the LDA functional.