The effect of strain and pressure on the electron-phonon coupling and superconductivity in MgB₂—Benchmark of theoretical methodologies and outlook for nanostructure design

Beginning with the discovery by Nagamatsu et al. in 2001 of $MgB2$ being a high-temperature superconductor with $Tc=39$ K,^1,2 a flurry of both experimental and theoretical studies on different properties of this compound was initiated. $MgB2$ has an $AlB2$-type hexagonal crystal structure with the space group of P6/mmm with alternating layers of pure Mg and graphene-like B, as seen in Fig. 1, with experimental lattice parameters $a=3.086$ Å and $c=3.524$ Å.¹ Each unit cell consists of one Mg atom and two B atoms. Before its discovery, the record-holding binary system was $Nb3$Ge with $Tc=23.2$ K, making the sudden jump in $Tc$ impressive.

Both the electronic and phonon band structures hold crucial information for understanding superconductivity. More exactly, the electronic states near the Fermi energy $ϵF$ are vital, as they are the ones that couple to the underlying vibrating lattice. Due to the difference in electronegativity, Mg will donate electrons to the B-layer, effectively making it isoelectronic to graphene. Consequently, the states at $ϵF$ are predominantly of B-character. The electronic band structure reveals two $σ$-bands from the valence band with $px$ and $py$ character and one $π$-band from the conduction band with $pz$ character. These make up the Fermi surface that consists of two narrow coaxial cylindrical sheets in the $Γ$–$A$ out-of-plane direction, from the $px,y$ $σ$-bands, and two distinct 3D tubular networks along in-plane directions, arising from bonding and antibonding $pz$ $π$-bands, see, e.g., the work of Kortus et al.³ for visualization.

Measurement of phonon band structure by inelastic x-ray scattering⁴ reveals Kohn anomalies⁵ in the optical $E2g$ branch coming from in-plane boron modes. The location of the Kohn anomalies is determined by the in-plane and out-of-plane nesting that occurs on the Fermi surface.^5,6 Theoretical studies have shown that the Kohn anomalies can directly be tied, by their depths, to the superconducting transition temperature not only in $MgB2$ but also in other $AlB2$-type diborides.^7–10 However, more studies are needed that compare methods for calculating $Tc$ based on the strength of the Kohn anomaly to calculations that explicitly treat the electron–phonon coupling (EPC), in particular, for more than just one given case, i.e., in a material where the different methodological approaches can be studied as a function of an external condition. We suggest that pressure in $MgB2$ is such an excellent example.

The effect of positive hydrostatic pressure on superconductivity in $MgB2$ has been measured and is found to suppress $Tc$⁠.^11–15 Zhang and Zhang¹⁶ studied the effect of negative pressure on $Tc$ using semi-empirical models and found that it enhances $Tc$⁠. On the contrary, Gerbi and Singh¹⁷ found in their first-principles calculations that positive hydrostatic pressure increases $Tc$⁠. The study is the only one showing such pressure behavior. Moreover, the effect of anisotropic strain on MgB$2$ superconductivity has been explored experimentally^18–20 and theoretically,^21–23 out of which only Ref. 22 treats EPC from first principles (for a monolayer). The general finding is that tensile strain boosts $Tc$ in $MgB2$⁠.

At ambient conditions, the theory of superconductivity in the Migdal–Éliashberg^24,25 formalism is known to perform well in estimating the superconducting transition in $MgB2$⁠, in both isotropic and anisotropic approximations.²⁶ Additionally, this theory relates $Tc$ to the phonon vibration frequency, which, in turn, is inversely proportional to the square root of the atomic mass. As such, one guiding principle is that materials with light atoms should have high $Tc$’s. $MgB2$⁠, just having one Mg atom and two light B atoms per unit cell, and no magnetism involved, is a perfect compound to use for methodological benchmarking.

In this work, we benchmark different approaches that predict $Tc$ of $MgB2$ under hydrostatic pressure and anisotropic stress and strain. The first studied approach relies on phonon dispersions and the pressure-induced evolution of the Kohn anomaly on the optical $E2g$ branch. Phonons derived from both finite displacements and linear response methods are compared and differences in the identified Kohn anomalies along $Γ$–$q$ (⁠$q=K,M,H,L$⁠) between the methodologies are quantified.

Then, we explicitly calculate EPC in the formalism of Migdal–Éliashberg^24,25 as a function of pressure and use the finite phonon linewidth to identify exactly where in the phonon dispersion the coupling occurs. $Tc$ is derived in both isotropic and anisotropic approximations and is compared to experimental values from the literature.

Lastly, the recent development of thin film deposition methods, synthesis of diborides, and theoretical calculations on the mixing and clustering thermodynamics of diborides²⁷ opens up for attempts to design nanostructures with $MgB2$ where both the $a$ and $c$ lattice parameters can be strained simultaneously. As such, we calculate the EPC and $Tc$ over a $5×5$ grid of individually stressed and strained $a$ and $c$ lattice parameters, with discrete steps of 2% and 7% in terms of equilibrium $a$ and $c$ lattice parameters, respectively. In particular, we suggest that $Tc$ can be boosted by introducing tensile lattice strain along both the $c$-axis, and in the $ab$-plane, by co-deposition of $MgB2$ with a combination of $YB2$ and either $ZrB2$ or $HfB2$ as (Zr,Y)$B2$ and (Hf,Y)$B2$⁠.

Electronic and phonon properties in $MgB2$ are derived from first-principles calculations in two approaches. First, in the framework of density functional theory (DFT),^28–31 as implemented in the Vienna ab initio Simulation Package (VASP).^32,33 Second, from density functional perturbation theory (DFPT),^34–37 as implemented in the Quantum-ESPRESSO (QE)^38–40 and Electron–Phonon coupling using Wannier functions (EPW)^41,42 simulation packages.

For VASP, the Projector-Augmented-Wave (PAW) method⁴³ is used to expand the electronic wave function in plane waves. The generalized gradient approximation (GGA) functional is used for calculating the exchange-correlation energies, as proposed by Perdew, Burke, and Ernzerhof (PBE96).⁴⁴ To ensure a sufficient energy- and force-convergence, a plane-wave energy cutoff of 600 eV is used. For the unit cell, a $k$-point mesh of $17×17×17$ is used when sampling the Brillouin zone in the Monkhorst–Pack scheme,⁴⁵ and a $5×5×5$ mesh is used for $6×6×6$ supercells in phonon calculations. The finite displacement (FD) method, as implemented in PHONOPY,⁴⁶ is used to calculate phonon frequencies and band structures. Atomic displacements of 0.01 Å from their equilibrium positions are performed for a symmetry reduced set of displacements, using the Parlinski–Li–Kawazoe method.⁴⁷ 

For QE and EPW, EPC and superconducting properties are explicitly calculated from first-principles using both linear-response and maximally localized Wannier functions. The mutual formalism is that of the Migdal–Éliashberg theory.^24,25 Electrons are described with both norm-conserving (NC)^48–50 and Vanderbilt ultrasoft (US)⁵¹ pseudopotentials (PP). A cutoff of 400 Ry is used for the charge density, and 55 and 80 Ry are used for USPP and NCPP electronic wave functions, respectively. In both cases, the charge density is integrated on a $Γ$-centered $12×12×$ 12 $k$-point mesh and a Gaussian smearing of 0.02 Ry is applied. First-order EP matrix elements are calculated using DFPT on a $Γ$-centered $6×6×6$ $q$-point mesh. Electron–phonon Wannier interpolation was performed on $48×48×48$ $k$- and $q$-meshes.

A brief description of central equations in the Migdal–Éliashberg^24,25 formalism of superconductivity is given below. If we let $n$ and $ν$ denote band indices of electrons and phonons, respectively, the electronic response upon the scattering of an electron from state $nk$ to $mk+q$ while emitting or absorbing a phonon $νq$ with frequency $ω$ is described by the matrix element

where $M$ is the reduced mass and $1/2Mωνq$ is the zero-point phonon amplitude, and $∂V/∂uνq$ is the induced potential per unit displacement, with $uνq$ being the phonon normal coordinate. The matrix element $g$ is central to electron–phonon theory. For instance, it is used to calculate the EPC strength associated with a specific phonon $νq$⁠,

where $N(ϵF)$ is the electronic density of states at the Fermi energy and $ΩBZ$ is the volume of the Brillouin zone. The energy of electron band $n$ with wavevector $k$ is given by $ϵnk$ and so on. The delta functions restrict the electron scattering to the Fermi surface. From this, one can calculate the total EPC $λ=∑νqwqλνq$ (⁠$wq$ being the BZ weights that normalize to 1), a quantity discussed in the results below. Furthermore, the phonon linewidth, corresponding to the imaginary part of the phonon self-energy, can easily be obtained from $γνq=2πN(ϵF)ωνq2λνq$⁠. This quantity identifies what modes and in which directions the EP scattering is most amplified, leading to decreased phonon lifetimes.

The lattice parameters $a$ and $c$ at hydrostatic pressures, ranging from $−$20 to 25 GPa in steps of 5 GPa, has been derived using VASP-PBE96, QE-US, and QE-NC approaches and is presented in Table I. Overall, the values are in excellent agreement, with VASP-PBE96 being slightly higher at the pressures below $−$10 GPa. Experimental values are derived from quadratic fit of measurements from the work of Bordet et al.,⁵² where pressures from ambient conditions up to 40 GPa were considered. Values in parentheses are extrapolated from the quadratic fit.

There is a certain type of phonon anomalies that is named Kohn anomalies, after Walter Kohn who first brought them into light in his paper from 1959.⁵ The message is that lattice vibrations of nuclei in metallic compounds are partially screened by the conduction electrons. This screening can change rapidly on certain surfaces in phonon $q$-space defined by $|q+K|=2kF$⁠, where $K$ is some reciprocal lattice vector and $kF$ is the Fermi momentum. Therefore, on these surfaces, the frequencies $ω$ vary abruptly with $q$⁠. The location of the anomalies, i.e., the location of these surfaces, is determined by the shape of the electronic Fermi surface. Kohn anomalies in the optical phonon branches of $MgB2$ and other $AlB2$-type structures have previously been theoretically studied at ambient conditions and are associated with its superconducting transition temperature.^7–10 The Kohn anomaly appears on the boron $E2g$ phonon branch related to in-plane B–B bond stretching vibrations. As will be shown later in this work, we find Kohn anomalies along the $Γ$–$q$ (⁠$q=K,M,H,L$⁠) directions and explicit evidence of strong EPC at their location. We also find strong EPC along $Γ$–$A$⁠, but without a clear sign of a Kohn anomaly. Calandra et al. note the presence of a weak anomaly at $q≈0.8ΓA$ that other works have failed to capture.⁶ This highlights intricacy of these phonon dispersion features. The Kohn anomalies along $Γ$–$q$ (⁠$q=K,M,H,L$⁠) are attributed to in-plane nesting on the Fermi surface, while along $Γ$–$A$ it is due to out-of-plane nesting between two coaxial holelike $σ$-band cylinders.⁶ See Ref. 3 for visualization of the Fermi surface of $MgB2$⁠.

In the works of Alarco et al.⁹ and Mackinnon et al.,¹⁰ the authors use the following equation that connects the anomaly depth $δ$ to the thermal energy $Tδ$⁠,

where $n$ is the degrees of freedom per atom, $N$ is the number of atoms per unit cell, $Z$ is the number of formula units per unit cell, $kB$ is Boltzmann’s constant, and $kBTδ2$ is the average thermal energy contribution per degree of freedom, as dictated by the equipartition theorem. For $MgB2$ with three atoms per unit cell, each with three degrees of freedom, Eq. (3) can be simply reduced to the following expression for the thermal energy:

This quantity has been shown to correlate with experimentally measured superconducting transition temperature $Tc$⁠, not only for $MgB2$⁠, but also for metal substituted $Mg1−xMxB2$⁠, where M = Al, Sc, and Ti.^9,10 This provides proof of concept that $Tδ$ can be used to estimate $Tc$⁠.

To further link the Kohn anomaly to superconductivity in $MgB2$⁠, we can borrow the ratio $2Δ=3.50kBTc$ from the seminal paper by Bardeen et al.⁵³ Here, $Δ$ denotes the superconducting gap size. Letting $Tc→Tδ$ from the discussion above and combining with Eq. (4) gives the estimate $Δ≈0.39δ$⁠. As soon will be shown, our calculations average to $δ=14.12$ meV at ambient conditions. Therefore, one could expect $Δ≈5.51$ meV. Experimentally, the superconducting gap size $Δ$ of the $σ$ and $π$ bands have been measured by Souma et al.⁵⁴ to $6.5±0.5$ and $1.5±0.5$ meV, respectively. The authors further note that the $σ$ band is dominant in the superconductivity of $MgB2$ and that the $π$ band is less important, having much weaker coupling to phonons. As such, the rough estimate of 5.51 meV is comparable to the $6.5±0.5$ meV of the $σ$-band.

In this work, as one of the approaches to predict $Tc$ of stressed/strained $MgB2$⁠, we now explore how the Kohn anomaly $δ$ changes with hydrostatic pressure. Using both FD and DFPT approaches to calculate phonon dispersions, the Kohn anomalies identified along the four q-directions $Γ$–$q$ (⁠$q=K,M,H,L$⁠) are studied as a function of pressure. Figure 2 summarizes the Kohn anomaly properties and their pressure dependence.

Figure 2(a) shows the differences in the phonon band structures of $MgB2$ at ambient conditions between the two calculation methods of FD and DFPT. All phonon branches are, for all practical purposes, identical except for the anomalous $E2g$ branch. This shows that numerical convergence has been achieved. The difference in the $E2g$ branch is mostly seen in the relative depth, here defined as $d$⁠, and is measured from the bottom of the Kohn anomaly and up to the value at $Γ$⁠, as marked by the arrows on the right side of $Γ$⁠. In FD, the phonons are derived from direct force calculations of all symmetrically inequivalent static displacements within a large supercell. In DFPT, one essentially calculates the linear response in the charge density in the unit cell, induced by the presence of a phonon with wave vector $q$⁠, in order to find the phonon dispersions. The discrepancy could be that the two methods are not equally sensitive to the intricate details of the Kohn anomalies.

The directional and pressure dependence of $d$ is shown in Fig. 2(b). For DFPT, the values are essentially invariant with pressure, up until 25 GPa where the phonon dispersion curves get distorted. For FD, the anomaly is more responsive to pressure, and after 10 GPa it becomes increasingly shallow relative to $Γ$⁠. For both cases, $Γ$–$L$ is found to form the deepest anomalous kink.

Figure 2(c) displays the anomaly depth $δ$⁠, which measures the vertical distance from the bottom of the Kohn anomaly up to its cusp,^9,10 along the different directions as a function of pressure. At ambient conditions, the average over all directions and between FD and DFPT turns out to be $δ=14.12$ meV. It becomes evident that increasing the pressure suppresses the Mexican hat-esque feature. In fact, at 25 GPa, both FD and DFPT predict that the Kohn anomalies are fully destroyed by the pressure. This further explains why the description of $d$ breaks down at high pressure. On the other hand, negative pressure, corresponding to tensile strain, effectively stretches the upper cusp (located at $∼$ 80 meV at zero pressure) away from the bottom of the anomaly. This is beneficial for enhancing the superconducting transition temperature according to Eq. (4).

$Tδ$ is then trivially calculated for each direction-dependent $δ$ and is presented in Fig. 3. From the average over all directions and both methods $⟨δ⟩=14.12$ meV at ambient conditions, we find $Tδ=36.4$ K. Individually, both methods are close to the experimental $Tc=39$ K at $p=0$ GPa and when averaging over all directions FD overestimates by just 2.54 K and DFPT underestimates by 7.70 K

Both methods capture the same trend with $Tδ$ decreasing with pressure. The green (⁠$Γ$–$M$⁠) and black (⁠$Γ$–$K$⁠) curves, both corresponding to in-plane q-directions, are overall on top of each other for both FD and DFPT. Meanwhile, the diagonal directions $Γ$–$H$ in red and $Γ$–$L$ in blue differ slightly. In terms of the equilibrium volume $V0$⁠, the pressure of $−$20 GPa corresponds to the volume $V=1.185V0$⁠, which increases $Tδ$ to 77.39 (79.48) K according to FD (DFPT). On the opposite end at 25 GPa (V = 0.873$V0$⁠) both methods predict that the Kohn anomalies vanish and thus also $Tδ$ is completely suppressed. Finally, we conclude that comparing to calculations that explicitly treat electron–phonon coupling, using phonon dispersion curves alone to estimate superconducting transition temperature $Tc$ is computationally very cheap.

To aid in the visualization, $γνq$ is presented as Gaussian broadening of the phonon frequency $ωνq$ in the phonon band structure for all branches of equilibrium $MgB2$ in Fig. 4(a), all calculated using the QE software.

This reveals that the anomalous band in the 60–80 meV range along $Γ$–$K$⁠, $Γ$–$M$⁠, $Γ$–$H$ and $Γ$–$L$⁠, originating from B atoms, indeed have strong EPC and contribute to superconducting properties. $Γ$–$A$ of the same band is also found to have strong coupling, despite the Kohn anomaly noted by Calandra et al.⁶ at $q≈0.8ΓA$ not being captured in our phonon calculations. Moreover, looking at the electronic band structure in Fig. 4(b) reveals three important boron bands highlighted in color; two $σ$-bands (blue and green) and one $π$-band (red). The two $σ$-bands with $px,y$ character straddle the $Γ$–$A$ line with negligible dispersion near $ϵF$⁠. Furthermore, they cross $ϵF$ parabolically along $Γ$–$q$ (⁠$q=K,M,H,L$⁠) providing the majority of the contribution to $N(ϵF)$ and the EPC. The features seen in Fig. 4 are evidencing the presence of Cooper-Pairs, from the coupling of optical $B$–$B$ bond-stretching vibrations in the 60–80 meV range (at ambient pressure) and the holes at the top of the $σ$-bands.

Another key quantity in the theory is the Migdal–Éliashberg spectral function, which in the isotropic approximation is given by

where more exactly $F(ω)$ is the phonon density of states and $α$ denotes an average over all phonons of energy $ω$ in the BZ. Partial $F(ω)$ as a function of hydrostatic pressure from $−$20 to 25 GPa, in steps of 5 GPa, is presented in Fig. 5(a). Mg is shown by dashed lines and B by solid lines. As pressure is increased from $−$20 GPa (black) to 25 GPa (red), the Mg and B peaks in the density of states become increasingly separated. At 25 GPa, the phonon bands of Mg and B are almost completely separated, with minor overlap in the 40–60 meV range. On the other hand, at $−$20 GPa with high tensile strain, the B reaches down to even the 20–30 meV region. The overall shift of each individual peak is related to the atomic mass, where the lighter B atoms evidently are more sensitive to pressure effects.

The full spectral function $α2F(ω)$ is shown in Fig. 5(b) for the same pressure set. Here, the cumulative electron–phonon mass enhancement parameter $λ(ω)$ is also shown. It is expressed as

where in the upper integration limit $ω→∞$ it equates to the total EPC $λ$ mentioned above. From our calculations, we find that the total coupling decreases with pressure. Using QE (EPW), at the lowest considered pressure of $−20$ GPa it goes from its maximum $λ=1.498$ (1.259), to 0.649 (0.621) at ambient pressure, and down to 0.476 (0.459) at 25 GPa. Following from the discussion on the effect of pressure on phonon frequencies and that the integrand in Eq. (6) is weighted by $ω′−1$⁠, it can be seen as a general principle that tensile strain (negative pressure) lowers the energy of the phonon modes, resulting in enhanced EPC.

Furthermore, as the figure shows, spectral function $α2F(ω)$ does not follow the overall shape of the phonon DOS $F(ω)$⁠, a feature not usually seen in superconductors, compare, e.g., the examples of high pressure face-centered cubic (fcc) and body-centered tetragonal (bct) B, and high pressure hexagonal close packed (hcp) Fe presented by Bose and Kortus.⁵⁵ The most striking feature of the spectral function is that the position of its peak aligns with the location of the anomalous B $E2g$ optical branch, dwarfing all contributions from other phonons—cf. the light-blue $p=0$ curve in Fig. 5(b) with Fig. 4(a).

Thus, if $α2F(ω)$ is known, the total EPC is easily calculated. It can then be used to estimate $Tc$ using the Allen–Dynes form of the McMillan Equation, replacing the prefactor $ΘD/1.45$ by $ωlog/1.2$⁠, found to be accurate for known materials with $λ<1.5$⁠,^56–58 

where $ωlog$ is the logarithmically averaged phonon frequency and $μ∗$ is the Anderson-Morel Coulomb pseudopotential,⁵⁹ which is difficult to calculate from first-principles. The latter can be empirically measured from tunneling experiments and typically take on values in the range of 0.10–0.20,^56–58 introducing a slight uncertainty in the determination of $Tc$ in this procedure, especially at moderate coupling strengths. For the purpose of this work, it is taken as a fitting parameter where 0.16 is selected as the “standard choice.” The results of this isotropic approximation are presented in Fig. 6.

For comparison, the direction-averaged $Tδ$ from both finite displacements and linear response phonons of Fig. 3 has been included, here marked by colored stars connected by dotted lines. Black symbols show experimental values taken from the literature^11–15 that if all fitted to a simple linear equation reveals that the transition temperature has the average slope of $dTc/dP=−0.98$ K/GPa from 0 to 25 GPa. Kohn anomaly estimates using $Tδ$ are shown to perform well at ambient conditions, but overestimate the slope of which experimental values drop with pressure.

The first attention is brought to the comparison between McMillan–Allen–Dynes $Tc$ [Eq. (7)] as calculated using QE and EPW both with the choice $μ∗=0.16$⁠, shown by blue and red solid squares, respectively. At ambient pressure, the QE (EPW) estimate is 9.45 (8.54) K, with the difference coming from EPW having systematically lower $λ$ than QE according to our calculations. Here, the results from NCPP are shown. A comparison between USPP and NCPP revealed virtually the same estimations, within 0.5 K from each other on average.

Since $μ∗$ can be freely varied in this framework, we used QE as a trial case and found that $μ∗=0.01$ allows for a good fit to the experimental value at ambient pressure, resulting in $Tc=38.7$ K, shown by blue open squares. Varying $μ∗$ changes the absolute values of $Tc$⁠, while the relative change with respect to the ambient pressure value, i.e., the pressure trend, is widely preserved. The fact that the standard choice of $μ∗$ severely underestimates $Tc$⁠, unless the unconventionally small value of 0.01 is selected, is here taken as a sign that $MgB2$ is not very well-described as an isotropic material. This is not surprising as it has been well established herein and in the literature that the superconductivity is mainly governed by in-plane vibrations.

Instead, we go beyond the isotropic approximation and calculate the closing of the distribution of the anisotropic superconducting gap $Δnk$ on the Fermi surface, obtained both from the real axis using Padé approximants⁶⁰ and the imaginary axis for comparison. In both cases, the standard choice of $μ∗=0.16$ was employed. Our calculations show that both the real and imaginary axes’ approaches yield the same values. In Fig. 6, the red solid circles represent the results from calculations on the real axis.

Regardless of the theoretical approach, the overall trend of a $Tc$ decreasing with pressure is mutual. There is, however, the work by Gerbi and Singh¹⁷ that also has investigated the effect of hydrostatic pressure on superconductivity of $MgB2$ using EPW. They find controversially that $Tc$ increases with pressure already in the range 0–25 GPa covered in this work and that it, for instance, is enhanced to 95 K at a pressure of 650 GPa. It is an interesting outlier in the literature that is well worth noting.

Within the first year or two of the discovery of superconductivity in $MgB2$⁠, the initial experiments that investigated the effect of hydrostatic pressure on $Tc$ used anvils to achieve isotropically compressed conditions. As seen in Fig. 6, all experiments clearly indicate that compressive strain leads to a decrease in the superconducting transition temperature. Literature search reveals several experimental and theoretical efforts that look into using strain to tune $Tc$ in $MgB2$⁠, not only by compression, but also by anisotropic tensile strain.

The experimental work of Serquis et al.,¹⁸ and later with the continuation by Liao et al.,¹⁹ found that 1.1% compressive strain due to the presence of 5% Mg vacancies reduces $Tc$ by around 2 K. Pogrebnyakov et al.²⁰ synthesized $MgB2$ with $∼$0.55% in-plane biaxial tensile strain achieved by epitaxial growth on SiC that increased $Tc$ by roughly 2 K up to 41.8 K. It was confirmed by Raman scattering measurements that the increase was due to softening of the bond-stretching $E2g$ phonon mode, following the strain-induced lowering of the corresponding phonon frequency.

On the theoretical side, Zheng and Zhu²¹ used a semi-empirical model with fitted parameters to study $Tc$ as a function of the strained lattice. Their findings indicate that separately increasing both $a$ and $c$ leads to $Tc$ enhancements, with the most significant way being increasing the volume with a fixed $c/a$ ratio. Furthermore, they suggested several substrates that can provide in-plane biaxial tensile strain via lattice mismatch. Some examples being SiC, $Si1+xC1−x$⁠, AlN, GaN, and $AlxGa1−x$N. For instance, AlN provides 1% strain, compared to the 0.5% from SiC, and enhances $Tc$ by 7%–8%. Bekaert et al.²² used DFPT to study strained $MgB2$ monolayers and found that reducing in-plane Mg–Mg distances with 4.5% biaxial compressive strain reduces $Tc$ by 11 K, from 20 K of the relaxed single monolayer, while 4.5% tensile strain boosts it up to 53 K. It is noted that the response to strain is much higher in $MgB2$ than in superconducting graphene.²² Zhai et al.²³ studied how strain alters the covalency in transition metal diborides and discussed implications for superconductivity. They note how stretching uniaxially along the $c$-axis can raise the energy of the $σ$-band and favorably reduce the frequency of the $E2g$ phonon. Here, it should be noted that biaxial strain induced by epitaxial growth has limitations in achieving increased $Tc$ in $MgB2$⁠. This is because an increased $a$- or $c$-parameter tends to lead to a decrease in the other with a compensating effect on $Tc$⁠. This is in addition to the effect of loss of epitaxial coherency by misfit dislocations in thick layers.

Using mechanical strain to change the superconducting properties of a material has also been studied in other materials systems. For instance, Cheng et al.⁶¹ performed DFPT calculations to probe EPC under strain engineering in $AlB2$⁠. Similar to what is found for $MgB2$⁠, they find uniaxial tensile strain along the $c$-axis to increase $Tc$ for $AlB2$⁠, while the opposite is indicated for compressive strain. Furthermore, they find that straining is more effective than carrier doping in tuning superconductivity.⁶¹ Bozovic et al.⁶² measured $La2CuO4$ under tensile strain on $SrTiO3$ substrates and found that $Tc$ is brought up from 25 to 40 K. On $LaSrAlO4$ substrates, providing compressive strain, it is enhanced to 51.5 K. They point out that the variation in $Tc$ in $La2CuO4$ does not originate from strain alone, but is also caused by sensitivity to oxygen content. Li et al.⁶³ investigated buckled triangle borophene and $β12$ borophene from first-principles using DFPT and found that strain can be used as a tool to tune $Tc$ for both boron polymorphs. Herrera et al.⁶⁴ measured an increase in the transition temperature of SrTiO$3$ attributed to tensile strain-induced modification of phonon modes. Lastly, Ghini et al.⁶⁵ found that superconductivity in single crystal FeSe is enhanced under compressive uniaxial strain.

The literature is in this sense rich in studies on strain engineering, but to the best of our knowledge, there exists no work with explicit treatment of EPC from first-principles that treats both anisotropic compressive and tensile strain, or any combinations thereof, and their implications for superconducting properties of $MgB2$⁠. It is, therefore, our intention to investigate this matter further here.

In the work of Alling et al.,²⁷ the configurational mixing and clustering tendencies of a large set of ternary metal diborides were studied. For the purpose of this study, we select the subset of $MgxM1−xB2$ alloys with metals $M$ that do not want to form solid solutions with Mg, but cluster instead. These metals are $M=V$⁠, Ti, Hf, Zr, and Y. In contrast to the aforementioned substrates like SiC and AlN, these all form $AlB2$-type diborides that are isostructural to $MgB2$⁠. We propose that these compounds can potentially be used in cleverly designed epitaxial growth to create layered $MgxM1−xB2$ structures that are under compressive or tensile strain along preferable crystal directions. For instance, if biaxial tensile or compressive strain in the $ab$-plane is desired, epitaxial growth along a common $c$-axis could allow for sandwiching of out-of-plane ordered pure Mg and suitable $M$ (with the flat boron sublayers being spectators in the clustering). This idea is similar to what has recently been experimentally proven to be possible in the work of Dahlqvist et al.⁶⁶ We also suggest co-deposition with intermediate temperature to allow some, but limited, surface diffusion. This allows for creating 3D intertwined coherent nanostructures or columns that are strained from all directions. Furthermore, in such mixtures, co-deposition can result in a metastable solid solution that upon annealing can undergo spinodal decomposition resulting in a coherent nanostructure with strained $MgB2$ domains.

In Fig. 7, the VASP-derived lattice parameters of isostructural $AlB2$-type diborides predicted to have clustering tendencies with respect to $MgB2$ are presented. For comparison, $AlB2$ itself is also included but has ordering tendencies. To sufficiently probe the lattice parameter space of $MgB2$⁠, and to capture the values of the closest-lying clustering candidates, a $5×5$ $(a,c)$ grid was constructed. The chosen grid spacings are, in terms of the equilibrium lattice parameters, 2% in the $a$-direction and 7% in the $c$-direction. This is represented by the purple points in the figure, corresponding to different pressures. Out of the metal diborides, one can readily identify that $HfB2$ and $ZrB2$ are candidates for achieving biaxial in-plane tensile strain, with an $∼2%$–$4%$ $a$-parameter mismatch, while $c$ remains matched within 1%. Using $VB2$ or $TiB2$ results in compression of both $a$ and $c$⁠. $YB2$ could potentially be used to achieve high tensile strain, but dynamical stability of $MgB2$ at such cell shape was not tested in this study. One can also consider co-deposition of $MgB2$ with a combination of $ZrB2$ and $YB2$ or $HfB2$ and $YB2$⁠, as (Zr,Y)$B2$ and (Hf,Y)$B2$ tend to form solid solutions while they all tend to avoid mixing with $MgB2$⁠.²⁷ Considering the complexity of controlling the B to metal content in thin film deposition, such experiments would benefit from recent advances in combinatorial growth, allowing mapping of the compositional space in one single experiment.

Despite not being explicitly shown here, the phonon band structure was derived to check for dynamical stability for all 25 grid points of pure $MgB2$ using both finite displacement and linear response methods. Finite displacements predict that all strained structures are dynamically stable, while the linear response is indicative of stretching the $c$-axis by $+14%$ results in imaginary frequencies at the $Γ$-point at 0 K. It is not investigated further in this study if temperature-driven anharmonic effects can stabilize such features. Regardless, synthesis under the right conditions is capable of promoting the growth of metastable thin films when atomic diffusion is hindered, as, for instance, in physical vapor deposition of (Ti,Al)$B2$⁠.^67–69 Available grid points are used to show qualitative predictions of what is possible with strain engineering.

In Fig. 8, the total EPC coupling $λ$ and the absolute change in $Tc$ relative to the equilibrium value, calculated using McMillan–Allen–Dynes Eq. (7) with $μ∗=0.16$⁠, are presented for the anisotropically strained $MgB2$⁠. A green point and transparent plane is used to guide the eye to the value corresponding to unstrained $MgB2$⁠. In both subfigures, it is clear that compressing either $a$ or $c$ alone, or in a combination as in the case of hydrostatic pressure, the EPC strength and consequently $Tc$ are suppressed. The most minute change to the lattice, captured by the grid, is achieved by stepping either $±$2% along the $a$-direction. A 2% compression of $a$ alone is here predicted to reduce $Tc$ by 2.6 K, while on the contrary a 2% expansion increases it by 5.5 K. Compressing $c$ alone by 7% is predicted to reduce $Tc$ by 8.3 K, and a 7% expansion increases it by 14.0 K. Stretching $a$ while compressing $c$ is found to suppress superconductivity. For instance, stretching $a$ by 2% while compressing $c$ by 7% lowers $Tc$ by 9.3 K. This is in contrast to the reverse case, whereby compressing $a$ by 2% and stretching $c$ by 7% increases $Tc$ by 11.3 K. In line with the indication from the semi-empirical model of Zheng and Zhu,²¹ the grid with explicit treatment of EPC shows how individually stretching either $a$ or $c$ alone or in unison enhances the superconductivity in $MgB2$⁠. Simultaneously stretching $a$ by 2% and $c$ by 7% enhances $Tc$ by 18.2 K. If $MgB2$ can remain stable at conditions corresponding to an increase in $a$ by 4% and $c$ by 14%, $Tc$ is here predicted to be increased by 38.4 K, effectively doubling the value of unstrained $MgB2$⁠.

The effects on superconducting transition temperature in $MgB2$ by hydrostatic pressure and anisotropic stress and strain conditions have been investigated. Different theoretical methodologies with varying complexity have been benchmarked.

First, two strictly phonon-based methods of finite displacements and linear response were employed to capture how the Kohn anomaly in the optical $E2g$ branch, corresponding to in-plane B–B bond stretching modes, evolves under hydrostatic pressure. The related Kohn anomaly depth $δ$ was found to decrease under compressive strain and increase under tensile strain. The thermal energy $Tδ$ tied to the anomaly has in previous studies been directly linked to superconducting transition temperature $Tc$ for $MgB2$ and other $AlB2$-type diborides. We find that enhancement of $Tc$ strongly favors negative pressures and that sufficiently large compression of the lattice will fully suppress the Kohn anomaly. For $MgB2$⁠, using phonon dispersion curves alone to estimate $Tc$ is computationally non-expensive and surprisingly accurate for its cost. Kohn anomalies are identified along $Γ$–$q$ (⁠$q=K,M,H,L$⁠) directions and arise due to Fermi surface nesting. Averaged across all directions, the Kohn anomaly approach estimates $Tc$ using finite displacements (linear response) to 41.54 (31.30) K at ambient pressure. This approach somewhat overestimates the magnitude of the slope $dTc/dP$⁠. The results presented herein motivate applying this methodology to other materials systems with anomalous phonon branches and studying their pressure behavior.

The calculations following the Migdal–Éliashberg formalism of EPC are computationally more expensive, but allow for explicit treatment of superconducting properties. Visualizing the phonon dispersion $ωνq$ broadened by the phonon linewdith $γνq$ reveals explicit evidence of EPC on the $E2g$ branch in the range 60–80 meV (at ambient pressure), along $Γ$–$q$ (⁠$q=A,K,M,H,L$⁠). These locations correspond to where both $σ$-bands graze the Fermi energy. Calculation of the Éliashberg spectral function $α2F(ω)$ as a function of pressure captures how the location of the anomalous branch moves during lattice strain and how the total EPC varies. Increasing pressure is found to stiffen the anomalous phonon frequency and decrease the total electron–phonon coupling, effectively suppressing the superconductivity in $MgB2$⁠. The superconducting transition temperature is estimated using the isotropic McMillan–Allen–Dynes formula and is found to underestimate $Tc$ by over 20 K when a typical value of the pseudopotential $μ∗=0.16$ is used. We show that $μ∗=0.01$ boosts $Tc$ to 38.7 K and that the pressure trend closely follows that of the majority of experimental reference values. Furthermore, we calculate $Tc$ by the closing of the superconducting gap $Δnk$ with $μ∗=0.16$ as a function of pressure and find the value 37 K at ambient pressure. Furthermore, the pressure trend is in excellent agreement with experiments in the literature, in particular, the work of Tissen et al.¹³ 

Literature review reveals a plethora of studies on the topic of using stress and strain engineering to enhance $Tc$ in $MgB2$ and a handful of other selected compounds, e.g., $AlB2$⁠, $La2CuO4$⁠, borophenes, $SrTiO3$⁠, and FeSe. Whether it is compressive or tensile strain that enhances $Tc$ is system dependent. The metal diborides $MB2$ (⁠$M=V,Ti,Hf,Zr,Y$⁠) are identified as potential candidates for clever nanostructure design as they are predicted to have clustering tendencies with respect to $MgB2$⁠. A $5×5$ $(a,c)$ grid is constructed to probe the lattice parameter space in the vicinity of the $MB2$ candidates. The total EPC and change in $Tc$ relative to the equilibrium value is calculated for all 25 points. The results show that biaxial in-plane tensile strain, uniaxial out-of-plane tensile strain, or a combination of both enhance $Tc$ by a considerable amount. Furthermore, both $HfB2$ and $ZrB2$ are identified as candidates for achieving tensile strain in cleverly designed 3D intertwined coherent nanostructures or columns that are strained in both $a$ and $c$⁠, something not achieved by planar epitaxy alone. $YB2$ would allow for the highest level of tensile strain, but dynamical stability of $MgB2$ was not investigated at those lattice parameters. A detailed study of the effects of quaternary diboride alloys and their nanostructures is fit for future investigations.

See the supplementary material for additional phonon dispersion relations $ωνq$⁠, with corresponding phonon linewidths $γνq$⁠. Presented is the hydrostatic pressures series ranging from $−$20 to 25 GPa, in steps of 5 GPa. Additionally, the nine grid points in the top right quadrant of Fig. 7, where anisotropic tensile strain is predicted to enhance the superconductivity the most, are also presented.

Financial support from the Knut and Alice Wallenberg (KAW) Foundation, through Project Grant No. KAW 2015.0043 is greatly acknowledged. B.A. acknowledges financial support from the Swedish Research Council (VR) through International Career Grant Nos. 2014-6336 and 2019-05403, from Marie Sklodowska Curie Actions, Cofund, Project INCA 600398, and from the Knut and Alice Wallenberg Foundation (Wallenberg Scholar Grant No. KAW-2018.0194), as well as support from the Swedish Foundation for Strategic Research through the Future Research Leaders 6 program, No. FFL 15-0290. J.R. and B.A. acknowledge the support from the Swedish Government Strategic Research Area in Materials Science on Functional Materials at Linköping University (Faculty Grant SFOMatLiU No. 2009 00971). A.E. gratefully acknowledges the support from Office of the Permanent Secretary, Ministry of Higher Education, Science, Research and Innovation: Research Grant for New Scholar (RGNS) through Grant No. RGNS 63-013. All calculations were performed using supercomputer resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre (NSC).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.