Electronic Transport in Weyl Semimetals with a Uniform Concentration of Torsional Dislocations

Abstract

In this article, we consider a theoretical model for a type I Weyl semimetal, under the presence of a diluted uniform concentration of torsional dislocations. By means of a mathematical analysis for partial wave scattering (phase-shift) for the T-matrix, we obtain the corresponding retarded and advanced Green’s functions that include the effects of multiple scattering events with the ensemble of randomly distributed dislocations. Combining this analysis with the Kubo formalism, and including vertex corrections, we calculate the electronic conductivity as a function of temperature and concentration of dislocations. We further evaluate our analytical formulas to predict the electrical conductivity of several transition metal monopnictides, i.e., TaAs, TaP, NbAs, and NbP.

1. Introduction

Weyl semimetals (WSMs) constitute a remarkable example of three-dimensional, gapless materials with nontrivial topological properties—first proposed theoretically [1,2,3,4,5,6,7] and, more recently, discovered experimentally on TaAs crystals [8].

In a WSM, the band structure possesses an even number of Weyl nodes with linear dispersion, where the conduction and valence bands touch. These nodes are monopolar sources of Berry curvature, and hence are protected from being gapped since their charge (chirality) is a topological invariant [7]. In the vicinity of these nodes, low energy conducting states behave as Weyl fermions, i.e., massless quasi-particles with pseudo-relativistic Dirac linear dispersion [4,5,6,7]. In Weyl fermions, conserved chirality determines the projection of spin over their momentum direction, a condition referred to as “spin-momentum locking”. While Type I WSMs are Lorentz covariant, this symmetry is violated in Type II WSMs, where the Dirac cones are strongly tilted [9].

The presence of Weyl nodes in the bulk spectrum determines the emergence of Fermi arcs [8], the chiral anomaly, and the chiral magnetic effect, among other remarkable properties [9]. Therefore, considerable attention has been paid to understand the electronic transport properties of WSMs [10,11,12]. For instance, there are recent works on charge transport [13] in the presence of spin–orbit coupled impurities [14], electrochemical [15] and nonlinear transport induced by Berry curvature dipoles [16]. Somewhat less explored are the effects of mechanical strain and deformations in WSMs. From a theoretical perspective, it has been proposed that different types of elastic strains can be modeled as gauge fields in WSMs [17,18,19]. In previous works, we have studied the combined effects of a single torsional dislocation and an external magnetic field on the electronic [20,21] and thermoelectric [20,22] transport properties of WSMs, using the Landauer ballistic formalism in combination with a mathematical analysis for the quantum mechanical scattering cross-sections [23].

In this work, we extend our previous analysis to study the case of a diluted, uniform concentration of torsional dislocations and its effects on the electrical conductivity of type I WSMs. In contrast to our former studies [20,21,22], here we employ the Kubo linear-response formalism at finite temperatures. This requires explicitly calculating the retarded and advanced Green´s functions for the system, including the multiple scattering events due to the random distribution of dislocation defects in the form of a disorder-averaged self-energy term. For this purpose, we first analyze the scattering phase shift arising from a single torsional dislocation, and then we obtain the corresponding (retarded and advanced) Green’s function in terms of the T-matrix elements by solving analytically the Lippmann–Schwinger equation. We further extend this analysis, by incorporating the effect of a random distribution of such dislocations, with a concentration $n_d$, in the form of a disorder-averaged self-energy into the corresponding Dyson’s equation. Finally, we analyze the correction due to the scattering vertex, and by including this additional contribution, we calculate the electrical conductivity from the Kubo formula, as a function of temperature and concentration of dislocations. We present explicit evaluations of our analytical expressions for the electrical conductivity as a function of temperature and concentration of dislocations $n_d$, for several materials in the family of transition metals’ monopnictides, i.e., TaAs, TaP, NbAs and NbP, where the corresponding microscopic parameters, estimated by ab initio methods, were reported in the literature [24,25,26].

2. Scattering by a Single Dislocation

As a continuum model for a type I WSM under the presence of a single dislocation defect, as depicted in Figure 1, we consider the Hamiltonian [22]

$${\widehat{H}}^{\xi }=\xi \hslash v_F\sigma ·(\mathbf{p}+e{\mathbf{A}}^{\xi })+{\sigma }_0{V}_0\delta (r-a)\equiv {\widehat{H}}_0^{\xi }+{\widehat{H}}_1^{\xi },$$

(1)

where

$$\begin{array}{cc}\hfill {\widehat{H}}_0^{\xi }& =\xi v_F\sigma ·\mathbf{p},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\widehat{H}}_1^{\xi }& =\xi e{v}_F\left(\sigma ·\widehat{\varphi }\right)\frac{1}{2}{B}^{\xi }r\Theta (a-r)+V_0\delta (r-a){\sigma }_0.\hfill \end{array}$$

(3)

Here, $\xi =\pm$ labels each of the Weyl nodes located at ${\mathbf{K}}_{\pm }=\pm \mathbf{b}/2$. The expression in Equation (2) is the free-particle Hamiltonian, whereas the expression in Equation (3) represents the interaction with the dislocation, where torsional strain is described as a pseudo-magnetic field inside the cylinder [21,22,23], as well as the lattice mismatch effect at the boundary of the dislocation, modeled as a repulsive delta barrier on its surface [22].

The “free” spinor eigenfunctions for the defect-free reference system satisfy

$${\widehat{H}}_0^{\xi }\left|{\Phi }_{\mathbf{k},\lambda }\right.〉={\mathcal{E}}_{\lambda ,\mathbf{k}}^{(0,\xi )}\left|{\Phi }_{\mathbf{k},\lambda }\right.〉,$$

(4)

where the energy spectrum is given by

$${\mathcal{E}}_{\lambda ,\mathbf{k}}^{(0,\xi )}=\lambda \xi \hslash v_F\left|\mathbf{k}\right|,$$

(5)

and $\lambda = \pm 1$ is the band (helicity) index. When projected onto coordinate space, these spinor eigenfunctions have the explicit form

$${\Phi }_{\lambda ,\mathbf{k}}\left(\mathbf{r}\right)=\sqrt{\frac{\left|\mathbf{k}\right|+\lambda k_z}{2\left|\mathbf{k}\right|}}\left(\begin{array}{c}1\\ \frac{\lambda \left|\mathbf{k}\right|-k_z}{k_x-i{k}_y}\end{array}\right)e^{i\mathbf{k}·\mathbf{r}},$$

(6)

and constitute an orthonormal basis for the Hilbert space.

If we now consider the (elastic) scattering effects induced by the torsional dislocation modeled by Equation (3), we need to look for the eigenvectors $\left|{\Psi }_{\lambda ,\mathbf{k}}\right.〉$ of the total Hamiltonian in Equation (1) with the same energy as in Equation (5). The answer is provided by the solution to the well-known Lippmann–Schwinger equation

$$\left|{\Psi }_{\mathbf{k},\lambda }\right.〉=\left|{\Phi }_{\mathbf{k},\lambda }\right.〉+{\widehat{G}}_{R,0}^{\xi }\left(E\right){\widehat{H}}_1^{\xi }\left|{\Psi }_{\mathbf{k},\lambda }\right.〉,$$

(7)

where the free Green’s function can be expressed in a coordinate-independent representation form via the resolvent,

$${\widehat{G}}_{R/A,0}^{\xi }\left(E\right)={\left[E-{\widehat{H}}_0^{\xi }\pm i{\eta }^{+}\right]}^{-1}.$$

(8)

Here, the index $R/A$ stands for retarded and advanced, respectively. As shown in detail in Appendix A, in the coordinate representation, the corresponding free Green’s function is given by the explicit matrix form ${\mathbf{G}}_{R,0}^{\xi }\left(\mathbf{r},{\mathbf{r}}^{\prime };k\right)=\delta (z-z^{\prime }){\mathbf{G}}_{R,0}^{\xi }\left(\mathbf{x},{\mathbf{x}}^{\prime };k\right)$, where $\mathbf{r}=(\mathbf{x},z)$ and

$$\begin{array}{c}\hfill {\mathbf{G}}_{R,0}^{\xi }\left(\mathbf{x},{\mathbf{x}}^{\prime };k\right)=-\frac{\lambda \xi ik}{4\hslash v_F}\left[\begin{array}{cc}{H}_0^{\left(1\right)}\left(k|\mathbf{x}-{\mathbf{x}}^{\prime }|\right)& i\lambda e^{-i\phi }{H}_1^{\left(1\right)}\left(k|\mathbf{x}-{\mathbf{x}}^{\prime }|\right)\\ i\lambda e^{i\phi }{H}_1^{\left(1\right)}\left(k|\mathbf{x}-{\mathbf{x}}^{\prime }|\right)& H_0^{\left(1\right)}\left(k|\mathbf{x}-{\mathbf{x}}^{\prime }|\right)\end{array}\right].\end{array}$$

(9)

Here, $H_0^{\left(1\right)}\left(z\right)$ and $H_1^{\left(1\right)}\left(z\right)$ are the Hankel functions, and $\mathbf{x}=(x,y)$ is the position vector on any plane perpendicular to the cylinder’s axis.

For the scattering analysis, we need the retarded resolvent for the full Hamiltonian, which is defined as the solution to the equation

$$\left(E+i{\eta }^{+}-{\widehat{H}}^{\xi }\right){\widehat{G}}_R^{\xi }\left(E\right)=\widehat{I}.$$

(10)

Combining Equation (10) with Equation (8), we readily obtain

$$\begin{array}{cc}\hfill {\widehat{G}}_R^{\xi }\left(E\right)& ={\widehat{G}}_{R,0}^{\xi }\left(E\right)+{\widehat{G}}_{R,0}^{\xi }\left(E\right){\widehat{H}}_1^{\xi }{\widehat{G}}_R^{\xi }\left(E\right)\hfill \\ \hfill & ={\widehat{G}}_{R,0}^{\xi }\left(E\right)+{\widehat{G}}_{R,0}^{\xi }\left(E\right){\widehat{T}}^{\xi }\left(E\right){\widehat{G}}_{R,0}^{\xi }\left(E\right),\hfill \end{array}$$

(11)

where we introduced the standard definition of the T-matrix operator ${\widehat{T}}^{\xi }\left(E\right)$ that can be formally expressed in closed form by

$$\begin{array}{cc}\hfill {\widehat{T}}^{\xi }\left(E\right)& ={\widehat{H}}_1^{\xi }+{\widehat{H}}_1^{\xi }{\widehat{G}}_{R,0}^{\xi }\left(E\right){\widehat{T}}^{\xi }\left(E\right)\hfill \\ \hfill & ={\widehat{H}}_1^{\xi }{\left(\widehat{I}-{\widehat{G}}_{R,0}^{\xi }\left(E\right){\widehat{H}}_1^{\xi }\right)}^{-1}.\hfill \end{array}$$

(12)

Using this definition, along with the property ${\widehat{H}}_1^{\xi }|{\Psi }_{\mathbf{k},\lambda }\rangle ={\widehat{T}}^{\xi }|{\Phi }_{\mathbf{k},\lambda }\rangle$, we obtain the Lippmann–Schwinger Equation (7) in the coordinate representation

$${\Psi }_{\mathbf{k},\lambda }\left(\mathbf{r}\right)={\Phi }_{\mathbf{k},\lambda }\left(\mathbf{r}\right)+\int d^3{r}^{\prime }\int d^3{r}^{″}〈\mathbf{r}|{\widehat{G}}_{R,0}^{\xi }\left(E\right)|{\mathbf{r}}^{\prime }〉〈{\mathbf{r}}^{\prime }\left|{\widehat{T}}^{\xi }\left(E\right)\right|{\mathbf{r}}^{″}〉{\Phi }_{\mathbf{k},\lambda }\left({\mathbf{r}}^{″}\right).$$

(13)

As shown in detail in Appendix B, by considering the asymptotic behavior of the Hankel functions, $H_{\nu }^{\left(1\right)}\left(x\right)\sim \sqrt{\frac{2}{\pi x}}{e}^{i\left(x-\frac{\nu \pi }{2}-\frac{\pi }{4}\right)}$ (for $x\to \infty$), Equation (13) can be reduced to the x–y plane and takes the explicit asymptotic expression

$${\Psi }_{{\mathbf{k}}_{\Vert },\lambda }\left(\mathbf{x}\right)\sim \frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\ \lambda \end{array}\right)e^{ikx}-\frac{\lambda \xi }{2\hslash v_F}\sqrt{\frac{ik}{\pi }}{T}_{{\mathbf{k}}_{\Vert }^{\prime }{\mathbf{k}}_{\Vert }}^{(\lambda ,\xi )}\left(\begin{array}{c}1\\ \lambda e^{i\varphi }\end{array}\right)\frac{e^{ikr}}{\sqrt{r}},$$

(14)

where, as we explain in the Appendix, the particles have only momenta perpendicular to the defect’s axis, i.e., ${\mathbf{k}}_{\Vert }=(k_x,k_y)$. Comparing this last result with our previous reported expression for the scattering amplitude [27],

$$\left[\begin{array}{c}{f}_1\left(\varphi \right)\\ f_2\left(\varphi \right)\end{array}\right]=\frac{e^{-\frac{i\pi }{4}}}{\sqrt{4\pi k}}\sum _{m=-\infty }^{\infty }\left[\begin{array}{c}{e}^{im\varphi }\\ \lambda e^{i(m+1)\varphi }\end{array}\right]\left(e^{2i{\delta }_m}-1\right),$$

(15)

we identify $T_{{\mathbf{k}}^{\prime }\mathbf{k}}^{(\lambda ,\xi )}=-2\lambda \xi \hslash v_F\sqrt{\pi /ik}{f}_1\left(\varphi \right)$. Therefore, we arrived at an explicit analytical expression for the T-matrix elements in terms of the phase shift ${\delta }_m\left(k\right)$ for each angular momentum channel $m\in \mathbb{Z}$

$$T_{{\mathbf{k}}_{\Vert }^{\prime }{\mathbf{k}}_{\Vert }}^{(\lambda ,\xi )}=-\frac{2\lambda \xi \hslash v_F}{k}\sum _{m=-\infty }^{\infty }{e}^{i{\delta }_m\left(k\right)}sin{\delta }_m\left(k\right)e^{im\varphi },$$

(16)

where $\varphi$ is the angle between ${\mathbf{k}}_{\Vert }$ and ${\mathbf{k}}_{\Vert }^{\prime }$, and the analytical expression for the phase shift is given in Appendix B by Equation (A27).

Figure 2. Pictorial description of the scattering event on a plane perpendicular to the cylindrical defect axis.

Figure 2. Pictorial description of the scattering event on a plane perpendicular to the cylindrical defect axis.

3. Scattering by a Uniform Concentration of Dislocations

Let us now consider a uniform concentration $n_d=N_d/A$ (per unit transverse surface) of identical cylindrical dislocations, as depicted in Figure 3, represented by the density function

$$\rho \left(\mathbf{x}\right)=\sum _{j=1}^{N_d}\delta (\mathbf{x}-{\mathbf{X}}_j),$$

(17)

where ${\mathbf{X}}_j$ is the position of the $j^{th}$-dislocation’s axis. The Fourier transform of this density function is thus given by the expression

$$\tilde{\rho }\left({\mathbf{k}}_{\Vert }\right)=\int d^{2}x{e}^{-i{\mathbf{k}}_{\Vert }·\mathbf{x}}\rho \left(\mathbf{x}\right)=\sum _{j=1}^{N_d}{e}^{-i{\mathbf{k}}_{\Vert }·{\mathbf{X}}_j}.$$

(18)

The operator that plays the role of a scattering potential for this distribution of dislocation defects is

$$V\left(\mathbf{x}\right)=\int d^2{x}^{\prime }\rho \left({\mathbf{x}}^{\prime }\right)H_1^{\xi }(\mathbf{x}-{\mathbf{x}}^{\prime })=\sum _{j=1}^{N_d}{H}_1^{\xi }(\mathbf{x}-{\mathbf{X}}_j),$$

(19)

where $H_1^{\xi }$ was defined in Equation (3) as the contribution from a single dislocation. The matrix elements of the scattering operator Equation (19) in the free spinor basis defined by Equation (4) are

$$\begin{array}{c}\hfill 〈{\Phi }_{{\mathbf{k}}_{\Vert },\lambda }\left|V\left(\mathbf{x}\right)\right|{\Phi }_{{\mathbf{k}}_{\Vert }^{\prime },{\lambda }^{\prime }}〉={\left[\tilde{V}({\mathbf{k}}_{\Vert }-{\mathbf{k}}_{\Vert }^{\prime })\right]}_{\lambda {\lambda }^{\prime }},\end{array}$$

(20)

where $\tilde{V}\left({\mathbf{k}}_{\Vert }\right)$ is the Fourier transform of $V\left(\mathbf{x}\right)$:

$$\begin{array}{cc}\hfill \tilde{V}\left({\mathbf{k}}_{\Vert }\right)& ={\int }_{{\mathbb{R}}^2}{d}^{2}x{e}^{-i{\mathbf{k}}_{\Vert }·\mathbf{x}}V\left(\mathbf{x}\right)=\sum _{j=1}^{N_d}{\int }_{{\mathbb{R}}^2}{d}^{2}x{e}^{-i{\mathbf{k}}_{\Vert }·\mathbf{x}}{H}_1^{\xi }(\mathbf{x}-{\mathbf{X}}_j)\hfill \\ \hfill & ={\tilde{H}}_1^{\xi }\left({\mathbf{k}}_{\Vert }\right)\tilde{\rho }\left({\mathbf{k}}_{\Vert }\right).\hfill \end{array}$$

(21)

Then, the matrix elements of the potential in Equation (20) become

$${\left[\tilde{V}\left({\mathbf{k}}_{\Vert }\right)\right]}_{\lambda {\lambda }^{\prime }}={\left[{\tilde{H}}_1^{\xi }\left({\mathbf{k}}_{\Vert }\right)\right]}_{\lambda {\lambda }^{\prime }}\tilde{\rho }\left({\mathbf{k}}_{\Vert }\right).$$

(22)

Let us also introduce the configurational average of a function $f\left({\mathbf{X}}_j\right)$ over the statistical distribution of dislocations as

$$〈f〉={\int }_{{\mathbb{R}}^2}{d}^2{X}_{j}P\left({\mathbf{X}}_j\right)f\left({\mathbf{X}}_j\right),$$

(23)

where $P\left({\mathbf{X}}_j\right)$ is the normalized distribution function for the defects in the sample. In particular, for a uniform distribution, we have $P\left({\mathbf{X}}_j\right)=1/A$, where A is the area of the plane normal to each cylinder’s axis. Now, the full retarded Green’s function under the presence of several dislocations represented by the operator $\widehat{V}$ given in Equation (19) satisfies the equation

$${\widehat{G}}_R^{\xi }\left(E\right)={\widehat{G}}_{R,0}^{\xi }\left(E\right)+{\widehat{G}}_{R,0}^{\xi }\left(E\right)\widehat{V}{\widehat{G}}_R^{\xi }\left(E\right).$$

(24)

The configurational average, as defined in Equation (23), of the full Green’s function in this last equation can be written as

$$〈{\widehat{G}}_R^{\xi }\left(E\right)〉={\widehat{G}}_{R,0}^{\xi }\left(E\right)+{\widehat{G}}_{R,0}^{\xi }\left(E\right){\Sigma }_R^{\lambda ,\xi }\left(E\right)〈{\widehat{G}}_R^{\xi }\left(E\right)〉.$$

(25)

This is the Dyson’s equation with the retarded self-energy ${\Sigma }_R^{\lambda ,\xi }\left(E\right)$ that can be explicitly solved to yield

$$〈G_R^{\lambda ,\xi }\left({\mathbf{k}}_{\Vert }\right)〉=\frac{1}{E-\lambda \xi \hslash v_F\left|{\mathbf{k}}_{\Vert }\right|-{\Sigma }_R^{\lambda ,\xi }\left({\mathbf{k}}_{\Vert }\right)}.$$

(26)

The effect of the statistical distribution of dislocations’ is entirely determined by the function $\tilde{\rho }\left({\mathbf{k}}_{\Vert }\right)$. In the perturbative expansion of the full Green’s function, we encounter nth-products of the form $\tilde{\rho }\left({\mathbf{k}}_{\mathbf{1}}\right)\tilde{\rho }\left({\mathbf{k}}_{\mathbf{2}}\right)\cdots \tilde{\rho }\left({\mathbf{k}}_n\right)$. The configurational average of a single factor is given by

$$\begin{array}{cc}\hfill 〈\tilde{\rho }\left({\mathbf{k}}_{\Vert }\right)〉& =〈\sum _{j=1}^{N_d}{e}^{-i{\mathbf{k}}_{\Vert }·{\mathbf{X}}_j}〉=\sum _{j=1}^{N_d}{\int }_{{\mathbb{R}}^2}{d}^2{X}_j\frac{1}{A}{e}^{-i{\mathbf{k}}_{\Vert }·{\mathbf{X}}_j}=\frac{N_d}{A}{\left(2\pi \right)}^2{\delta }^{\left(2\right)}\left({\mathbf{k}}_{\Vert }\right),\hfill \end{array}$$

(27)

Similarly, for the product of two factors, we obtain

$$\begin{array}{cc}\hfill 〈\tilde{\rho }({\mathbf{k}}_{\Vert }^1)\tilde{\rho }({\mathbf{k}}_{\Vert }^2)〉& =〈\sum _{j=1}^{N_d}\sum _{l=1}^{N_d}{e}^{-i{\mathbf{k}}_{\Vert }^1·{\mathbf{X}}_j-i{\mathbf{k}}_{\Vert }^2·{\mathbf{X}}_l}〉=〈\sum _{j=l}{e}^{-i({\mathbf{k}}_{\Vert }^1+{\mathbf{k}}_{\Vert }^2)·{\mathbf{X}}_j}+\sum _{j\ne l}{e}^{-i{\mathbf{k}}_{\Vert }^1·{\mathbf{X}}_j-i{\mathbf{k}}_{\Vert }^2·{\mathbf{X}}_l}〉\hfill \\ \hfill & =\frac{N_d}{A}{\left(2\pi \right)}^2{\delta }^{\left(2\right)}({\mathbf{k}}_{\Vert }^1+{\mathbf{k}}_{\Vert }^2)+\frac{N_d(N_d-1)}{A^2}{\left(2\pi \right)}^4{\delta }^{\left(2\right)}({\mathbf{k}}_{\Vert }^1){\delta }^{\left(2\right)}({\mathbf{k}}_{\Vert }^2),\hfill \end{array}$$

(28)

and we have a similar behavior for higher order products. Now, notice that, for $N_d\gg 1$, we have $N_d(N_d-1)\approx N_d^2$, $N_d(N_d-1)(N_d-3)\approx N_d^3$ and so on. We define the concentration of defects, i.e., the number of dislocations per unit of area perpendicular to the cylinder’s axis as $n_d=N_d/A$. As discussed in standard references [28,29], for small concentrations $n_d\ll 1$, the scaling discussed before ensures that the total Green’s function in Equation (25) can be calculated accurately by the sequence of diagrams for the retarded self-energy in momentum space as given in Figure 4, an approach well known as the non-crossing approximation (NCA). This series of diagrams corresponds to the configurational average of the T-matrix over the random distribution of dislocations after Equation (23)

$${\Sigma }_R^{\lambda ,\xi }\left(E\right)=〈{\widehat{T}}^{\xi }\left(E\right)〉=n_d{T}_{{\mathbf{k}}_{\Vert }{\mathbf{k}}_{\Vert }}^{(\lambda ,\xi )}.$$

(29)

Using the expression in Equation (16) for the T-matrix elements, for ${\mathbf{k}}_{\Vert }={\mathbf{k}}_{\Vert }^{\prime }$ that implies $\varphi =0$, we have that the real part of the self-energy

$$\mathrm{Re}{\Sigma }_R^{\lambda ,\xi }\left({\mathbf{k}}_{\Vert }\right)=-\frac{2\lambda \xi n_d\hslash v_F}{k}\sum _{m=-\infty }^{\infty }cos{\delta }_m^{(\lambda ,\xi )}\left(k\right)sin{\delta }_m^{(\lambda ,\xi )}\left(k\right)$$

(30)

contains an infinite sum over highly oscillatory terms that converges to zero. Therefore, no contribution arises from the real part of the self-energy. The imaginary part, on the other hand, leads to the definition of the relaxation time,

$$\frac{1}{{\tau }^{(\lambda ,\xi )}\left(k\right)}=-\frac{2\lambda \xi }{\hslash }{n}_d\mathrm{Im}{T}_{{\mathbf{k}}_{\Vert }{\mathbf{k}}_{\Vert }}^{(\lambda ,\xi )}.$$

(31)

3.1. Electrical Conductivity in the Linear-Response Regime

In order to arrive at the definition of the electrical conductivity in the linear response regime, let us first consider a single Fourier mode for an external electric field $\mathbf{E}$, in the gauge

$$\mathbf{E}=-\frac{\partial }{\partial t}\mathbf{A}(\mathbf{r},t),$$

(32)

where $\mathbf{A}(\mathbf{r},t)=\mathbf{A}(\mathbf{r},\omega )e^{-i\omega t}$ is the vector potential. Then, $\mathbf{E}=i\omega \mathbf{A}$. In the linear response formalism, the current is given by the expression

$$j_{\alpha }(\mathbf{r},\omega )=\int d^3{r}^{\prime }{\sigma }_{\alpha \beta }(\mathbf{r},{\mathbf{r}}^{\prime };\omega )E_{\beta }({\mathbf{r}}^{\prime }.\omega ),$$

(33)

where the conductivity tensor is defined by

$${\sigma }_{\alpha \beta }(\mathbf{r},{\mathbf{r}}^{\prime };\omega )=\frac{1}{i\omega }{K}_{\alpha \beta }(\mathbf{r},{\mathbf{r}}^{\prime };\omega ).$$

(34)

In the Kubo formalism, the tensor $K_{\alpha \beta }$ is defined in terms of the retarded current-current correlator as follows:

$$K_{\alpha \beta }(\mathbf{r},t;{\mathbf{r}}^{\prime },t^{\prime })=i{\hslash }^{-1}\theta (t-t^{\prime })\mathrm{Tr}\left\{\widehat{\rho }\left[{\widehat{j}}_{\alpha }(\mathbf{r},t),{\widehat{j}}_{\beta }({\mathbf{r}}^{\prime },t^{\prime })\right]\right\},$$

(35)

where $\widehat{\rho }$ is the statistical density matrix operator. As shown in detail in Appendix D, the Fourier transform of this tensor to the frequency domain is given by

$$K_{\alpha \beta }^{\xi }(\mathbf{r},{\mathbf{r}}^{\prime };\omega )=e^2{v}_F^2{\int }_{-\infty }^{\infty }\frac{d{E}^{\prime }}{2\pi }{\int }_{-\infty }^{\infty }\frac{dE}{2\pi }\frac{f_0\left(E^{\prime }\right)-f_0\left(E\right)}{\hslash \omega +E-E^{\prime }+i{\eta }^{+}}\mathrm{Tr}\left[{\sigma }_{\alpha }{𝓐}^{\xi }(\mathbf{r},{\mathbf{r}}^{\prime };E^{\prime }){\sigma }_{\beta }{𝓐}^{\xi }({\mathbf{r}}^{\prime },\mathbf{r};E)\right].$$

(36)

Here, $f_0\left(E\right)={\left[e^{(E-\mu )/kT}+1\right]}^{-1}$ is the Fermi distribution, and we introduced the (disorder-averaged) spectral function

$$\begin{array}{c}\hfill {\mathcal{A}}^{\lambda ,\xi }\left(k\right)=i\left[〈G_R^{\lambda ,\xi }\left({\mathbf{k}}_{\Vert }\right)〉-〈G_A^{\lambda ,\xi }\left({\mathbf{k}}_{\Vert }\right)〉\right]=\frac{2\left(\frac{\hslash }{2{\tau }^{(\lambda ,\xi )}\left(k\right)}\right)}{{\left(E-{\mathcal{E}}_k^{\lambda ,\xi }\right)}^2+{\left(\frac{\hslash }{2{\tau }^{(\lambda ,\xi )}\left(k\right)}\right)}^2}\end{array}$$

(37)

that clearly represents a Lorentzian distribution whose spectral width is defined by the inverse of the relaxation time (see Appendix C for the details). After some algebraic manipulations, we obtain the conductivity tensor at finite frequency and temperature

$$\begin{array}{cc}\hfill \Re {\sigma }_{\alpha \beta }^{\xi }(\mathbf{r},{\mathbf{r}}^{\prime };\omega )=& -\frac{e^2\hslash v_F^2}{2\pi }{\int }_{-\infty }^{\infty }dE\left\{\frac{f_0(E+\hslash \omega )-f_0\left(E\right)}{\hslash \omega }\right\}\hfill \\ \hfill & \times \mathrm{Tr}\left[{\sigma }_{\alpha }{𝓐}^{\xi }(\mathbf{r},{\mathbf{r}}^{\prime };E+\hslash \omega ){\sigma }_{\beta }{𝓐}^{\xi }({\mathbf{r}}^{\prime },\mathbf{r};E)\right].\hfill \end{array}$$

(38)

Using the coordinates representation of the spectral function given in Appendix C, after Equation (A33), we can read off the Fourier transform to momentum space of the conductivity

$$\begin{array}{cc}\hfill {\sigma }_{\alpha \beta }^{\xi }(\mathbf{q};\omega )=& -\frac{e^2\hslash v_F^2}{2\pi }\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{\int }_{-\infty }^{\infty }dE\left\{\frac{f_0(E+\hslash \omega )-f_0\left(E\right)}{\hslash \omega }\right\}\hfill \\ \hfill & \times \sum _{\lambda ,{\lambda }^{\prime }}\mathrm{Tr}\left\{{\sigma }_{\alpha }\left({\sigma }_0+\lambda \frac{\sigma ·({\mathbf{k}}_{\Vert }+\mathbf{q})}{|{\mathbf{k}}_{\Vert }+\mathbf{q}|}\right){\sigma }_{\beta }\left({\sigma }_0+{\lambda }^{\prime }\frac{\sigma ·{\mathbf{k}}_{\Vert }}{|{\mathbf{k}}_{\Vert }|}\right)\right\}\hfill \end{array}$$

$$\begin{array}{c}\hfill \times {\mathcal{A}}^{\lambda ,\xi }\left(\right|{\mathbf{k}}_{\Vert }+\mathbf{q}|;E+\hslash \omega ){\mathcal{A}}^{{\lambda }^{\prime },\xi }\left(\right|{\mathbf{k}}_{\Vert }|;E).\end{array}$$

(39)

We are interested in the DC conductivity, so we take the limit $\mathbf{q}\to \mathbf{0}$ first and then the limit $\omega \to 0$. After a long calculation (details in the Appendix D), the result is

$$\begin{array}{cc}\hfill {\sigma }_{\alpha \beta }^{(\lambda ,\xi )}\left(T\right)& ={\delta }_{\alpha \beta }\frac{e^2\hslash v_F^2}{{\pi }^3}{\int }_0^{\infty }dk{\int }_{-\infty }^{\infty }dE\left(-\frac{\partial f_0\left(E\right)}{\partial E}\right)〈G_R^{\lambda ,\xi }\left({\mathbf{k}}_{\Vert }\right)〉〈G_A^{\lambda ,\xi }\left({\mathbf{k}}_{\Vert }\right)〉{\mathbf{k}}_{\Vert }·{\mathbf{k}}_{\Vert }.\hfill \end{array}$$

(40)

3.2. Vertex Corrections

The self-energy contribution modifies the definition of the retarded and advanced Green´s functions in Equation (40), as depicted by the double lines in Figure 5b. However, there are also scattering processes involving links between the two internal Green function lines, as depicted in Figure 5a. When considering such diagrams with cross-links, as in Figure 5a, we must include the vertex correction as depicted in Figure 5b.

Taking into account the vertex correction, the conductivity becomes

$$\begin{array}{cc}\hfill {\sigma }_{\alpha \beta }^{(\lambda ,\xi )}\left(T\right)& ={\delta }_{\alpha \beta }\frac{e^2\hslash v_F^2}{{\pi }^3}{\int }_0^{\infty }dk{\int }_{-\infty }^{\infty }dE\left(-\frac{\partial f_0\left(E\right)}{\partial E}\right)〈G_R^{\lambda ,\xi }\left({\mathbf{k}}_{\Vert }\right)〉〈G_A^{\lambda ,\xi }\left({\mathbf{k}}_{\Vert }\right)〉{\mathbf{k}}_{\Vert }·{\Gamma }_{RA}({\mathbf{k}}_{\Vert },E),\hfill \end{array}$$

(41)

where the vertex function ${\Gamma }_{RA}({\mathbf{k}}_{\Vert },E)$ is given as the solution to the Bethe-Salpeter equation as depicted in Figure 6. Then, we have

$$\begin{array}{c}\hfill {\Gamma }_{RA}({\mathbf{k}}_{\Vert },E)={\mathbf{k}}_{\Vert }+n_d\int \frac{d^2{k}^{\prime }}{{\left(2\pi \right)}^2}〈G_R^{\lambda ,\xi }({\mathbf{k}}_{\Vert }^{\prime })〉〈G_A^{\lambda ,\xi }({\mathbf{k}}_{\Vert }^{\prime })〉{\left|T_{{\mathbf{k}}_{\Vert }^{\prime }{\mathbf{k}}_{\Vert }}^{(\lambda ,\xi )}\right|}^2{\Gamma }_{RA}({\mathbf{k}}_{\Vert }^{\prime },E).\end{array}$$

(42)

The iterative solution of Equation (42) for ${\Gamma }_{RA}({\mathbf{k}}_{\Vert },E)$ shows that the vertex function must be of the form

$${\Gamma }_{RA}({\mathbf{k}}_{\Vert },E)=\gamma ({\mathbf{k}}_{\Vert },E){\mathbf{k}}_{\Vert }.$$

(43)

Then, we obtain a secular integral equation for the scalar function $\gamma ({\mathbf{k}}_{\Vert },E)$ that in the low concentration limit becomes

$$\begin{array}{c}\hfill \gamma ({\mathbf{k}}_{\Vert },E)=1+n_d\frac{2\pi }{\hslash }\int \frac{d^2{k}^{\prime }}{{\left(2\pi \right)}^2}{\tau }^{(\lambda ,\xi )}\left(k^{\prime }\right){\left|T_{{\mathbf{k}}_{\Vert }^{\prime }{\mathbf{k}}_{\Vert }}^{(\lambda ,\xi )}\right|}^2\delta (E-\lambda \xi \hslash v_F{k}^{\prime })\gamma ({\mathbf{k}}_{\Vert }^{\prime },E)\frac{{\mathbf{k}}_{\Vert }·{\mathbf{k}}_{\Vert }^{\prime }}{k^2}.\end{array}$$

(44)

In the limit of low concentrations, we use the result in Appendix D, Equation (A46), to obtain

$$\begin{array}{c}\hfill {\sigma }_{\alpha \alpha }^{(\lambda ,\xi )}\left(T\right)=\frac{2e^2{v}_F^2}{{\pi }^2}{\int }_0^{\infty }dk{k}^2{\left(-\frac{\partial f_0\left(E\right)}{\partial E}\right)}_{E=\lambda \xi \hslash v_{F}k}{\tau }^{(\lambda ,\xi )}\left(k\right)\gamma ({\mathbf{k}}_{\Vert },\lambda \xi \hslash v_{F}k).\end{array}$$

(45)

At low temperatures, an exact solution is possible since the derivative of the Fermi distribution takes a compact support at the Fermi energy. Therefore, we can evaluate $\gamma \left(k\right)$ and ${\tau }^{(\lambda ,\xi )}\left(k\right)$ at the Fermi momentum $k_F^{\xi }$, to obtain

$$\gamma (k_F^{\xi })=\frac{{\tau }_1^{(\lambda ,\xi )}(k_F^{\xi })}{{\tau }_1^{(\lambda ,\xi )}(k_F^{\xi })-{\tau }^{(\lambda ,\xi )}(k_F^{\xi })},$$

(46)

where we defined (for $cos{\varphi }^{\prime }={\mathbf{k}}_{\Vert }·{\mathbf{k}}_{\Vert }^{\prime }/k^2$)

$$\frac{1}{{\tau }_1^{(\lambda ,\xi )}(k_F^{\xi })}=n_d\frac{2\pi }{\hslash }\int \frac{d^2{k}^{\prime }}{{\left(2\pi \right)}^2}{\left|T_{{\mathbf{k}}_{\Vert }^{\prime }{\mathbf{k}}_{\Vert }}^{(\lambda ,\xi )}\right|}^{2}cos{\varphi }^{\prime }\delta (\hslash v_F{k}_F^{\xi }-\hslash v_F{k}^{\prime }).$$

(47)

After the substitution of $\gamma ({\mathbf{k}}_{\Vert },E)$ from Equation (46) into Equation (45), we finally obtain a closed analytical expression for the bulk electrical conductivity, as a function of temperature and concentration of dislocations $n_d$

$$\begin{array}{cc}\hfill {\sigma }_{\alpha \alpha }^{(\lambda ,\xi )}\left(T\right)& =\frac{2e^2{v}_F^2}{{\pi }^2{k}_{B}T}{\tau }_{\mathrm{tr}}^{(\lambda ,\xi )}(k_F^{\xi }){\int }_0^{\infty }dk{k}^2{f}_0\left({\mathcal{E}}_{{\mathbf{k}}_{\Vert }}^{\lambda ,\xi }\right)\left[1-f_0\left({\mathcal{E}}_{{\mathbf{k}}_{\Vert }}^{\lambda ,\xi }\right)\right]\hfill \\ \hfill & =-\frac{4}{{\pi }^2{v}_F}\left(\frac{e^2}{\hslash }\right){\left(\frac{k_{B}T}{\hslash }\right)}^2{\tau }_{\mathrm{tr}}^{(\lambda ,\xi )}(k_F^{\xi }){\mathrm{Li}}_2\left(-e^{\frac{\hslash v_F{k}_F}{k_{B}T}}\right),\hfill \end{array}$$

(48)

where ${\mathrm{Li}}_2\left(x\right)$ is the polylogarithm of order 2. Here, the total transport relaxation time is defined by

$$\begin{array}{cc}\hfill \frac{1}{{\tau }_{\mathrm{tr}}^{(\lambda ,\xi )}(k_F^{\xi })}& =\frac{1}{{\tau }^{(\lambda ,\xi )}(k_F^{\xi })}-\frac{1}{{\tau }_1^{(\lambda ,\xi )}(k_F^{\xi })}\hfill \\ \hfill & =\frac{2\pi n_d}{\hslash }\int \frac{d^2{k}^{\prime }}{{\left(2\pi \right)}^2}\delta (\hslash v_F{k}_F^{\xi }-\hslash v_F{k}^{\prime }){\left|T_{{\mathbf{k}}_{\Vert }^{\prime }{\mathbf{k}}_{\Vert }}^{(\lambda ,\xi )}\right|}^2(1-cos{\varphi }^{\prime }).\hfill \end{array}$$

(49)

Using the analytical expression for the T-matrix elements in Equation (16), we obtain a closed expression for the transport relaxation time in terms of the scattering phase shifts ${\delta }_m\left(k\right)$

$$\frac{1}{{\tau }_{\mathrm{tr}}^{(\lambda ,\xi )}(k_F^{\xi })}=\frac{2n_d{v}_F}{k_F^{\xi }}\sum _{m=-\infty }^{\infty }{sin}^2\left[{\delta }_m(k_F^{\xi })-{\delta }_{m-1}(k_F^{\xi })\right].$$

(50)

From Equation (48), we can investigate the zero temperature $T\to 0$ and high temperature $T\gg \hslash v_F{k}_F/k_B$ limits, respectively. In the zero temperature limit, we obtain

$${\sigma }_{\alpha \alpha }^{(\lambda ,\xi )}(T\to 0)=\frac{2}{{\pi }^2}{k}_F^{\xi 2}\left(\frac{e^2}{\hslash }\right)v_{F,\alpha }{\tau }_{\mathrm{tr}}^{(\lambda ,\xi )}(k_F^{\xi }),$$

(51)

a constant that depends on the microscopic material properties (such as $v_F$), as well as on the concentration of dislocations $n_d$ through the relaxation time.

On the other hand, in the high-temperature limit $T\gg \hslash v_F{k}_F/k_B$, we obtain a quadratic dependence on temperature

$${\sigma }_{\alpha \alpha }^{(\lambda ,\xi )}(T\gg \hslash v_F{k}_F/k_B)=\frac{1}{3v_{F,\alpha }}\left(\frac{e^2}{\hslash }\right){\left(\frac{k_{B}T}{\hslash }\right)}^2{\tau }_{\mathrm{tr}}^{(\lambda ,\xi )}(k_F^{\xi }),$$

(52)

where the overall constant depends on the microscopic parameters for each material, as well as on the concentration of dislocations through the relaxation time.

4. Results

In this section, we apply the theory and analytical expressions obtained in the previous section to calculate the electrical conductivity of several materials in the family of transition metals’ monopnictides, i.e., TaAs, TaP, NbAs, and NbP. For an estimation of the concentration of defects $n_d$ in real crystal systems, Ref. [24] reports that the native concentration of dislocations in the lattice of the materials TiO₂ and SrTiO₃ varies in the range $n_d\sim {10}^5$–${10}^7{\mathrm{cm}}^{-2}$. These concentrations can be enhanced using different treatments up to ${10}^{13}{\mathrm{cm}}^{-2}$, close to the rendering amorphous limit. The microscopic/atomistic parameters involved in our theory are obtained from ab-initio studies for WSM materials, as reported in Refs. [25,26]. In particular, the later reference identifies anisotropies in the Fermi velocities and density of charge carries at different Weyl nodes and bands. Using these results for the densities of carriers, we compute the Fermi momentum at each Weyl node, i.e., $k_F^{\xi }$, as displayed in

Table 1. Values of $k_F^{\xi }$ computed from the carrier densities reported in Ref. [26].

Material	${\mathit{k}}_{\mathit{F}}^{+}[{\mathrm{nm}}^{-1}]$	${\mathit{k}}_{\mathit{F}}^{-}[{\mathrm{nm}}^{-1}]$
TaAs	0.23	0.05
TaP	0.50	0.09
NbAs	0.46	0.03
NbP	1.04	0.15

.

In what follows, for definiteness, we shall assume that the axis of the defects is along the crystallographic z-direction and that we are measuring the conductivity along the x-direction. Then, we use the reported x-components of the Fermi velocities [25,26]. We have different Fermi velocities $v_{F,x}^{(\lambda ,\xi )}$, for the conduction band ($\lambda =+ 1$) and for the valence band ($\lambda =-1$), and for each of the Weyl nodes ($\xi =\pm$), respectively. Actually, for the valence band, Refs. [25,26] report the hole velocity. Their results are presented in

Table 2. Values of the Fermi velocity $v_{F,x}^{(\lambda ,\xi )}$ in units of ${10}^5$ m/s, as reported in Ref. [26]. Notice that, for the valence bands ($\lambda =-1$), they report the hole velocity.

Material	${\mathit{v}}_{\mathit{F},\mathit{x}}^{(+,+)}$	${\mathit{v}}_{\mathit{F},\mathit{x}}^{(-,+)}$	${\mathit{v}}_{\mathit{F},\mathit{x}}^{(+,-)}$	${\mathit{v}}_{\mathit{F},\mathit{x}}^{(-,-)}$
TaAs	3.2	−5.3	2.6	−4.3
TaP	3.7	−5.4	2.0	−3.9
NbAs	3.0	−4.8	2.5	−3.2
NbP	3.0	−5.1	1.7	−2.4

.

Now, in order to study the additional effect of the torsional dislocations, we follow our previous work [22] to estimate the geometrical parameters involved in the model. We assume that the dislocations are cylindrical regions along the z-axis with radius a. Here, we further assume that the defects possess an average radius of $a=15$ nm. The simple relation between the torsional angle $\theta$ (in degrees) and the pseudo-magnetic field representing strain is $B_S{a}^2=1.36\theta {\tilde{\varphi }}_0$ [22], where the modified flux quantum in this material is approximately ${\displaystyle {\tilde{\varphi }}_0\equiv \frac{\hslash v_F}{e}=\frac{1}{2\pi }\frac{v_F}{c}\frac{hc}{e}=\frac{1}{2\pi }\frac{1.5}{300}·4.14\times {10}^5\mathrm{T}{\mathring{\mathrm{A}}}^2\approx 330 \mathrm{T}{\mathring{\mathrm{A}}}^2}$. In this work, we have chosen a torsion angle $\theta ={15}^{\circ }$. The lattice mismatch effect at the surface of the dislocation cylinders is modeled by a repulsive delta-potential, with strength $V_0$, expressed in terms of the “spinor rotation” angle $\alpha =V_0/\hslash v_F$. According to our previous work [22], a realistic choice is $\alpha =3\pi /4$.

With all of these parameters fixed, we can compute the transport relaxation time for each material from Equation (50). Our results are presented in

Table 3. Computed values for the total relaxation time and the total transport relaxation time for each material after Equation (50). We consider a concentration of dislocations $n_d={10}^{11} {\mathrm{cm}}^{-2}$.

Material	$\mathit{\tau }$ [${10}^{-13}$ s]	${\mathit{\tau }}_{\mathbf{tr}}$ [${10}^{-13}$ s]
TaAs	2.2	2.6
TaP	2.4	3.2
NbAs	2.2	3.1
NbP	2.4	4.2

.

Now, we compute the conductivity along the x-direction ${\sigma }_{xx}$. In what follows, we simply call it $\sigma \left(T\right)$, as a function of temperature. The total conductivity is the sum over nodes and bands

$$\sigma \left(T\right)=\sum _{\xi =\pm 1}\sum _{\lambda =\pm 1}{\sigma }_{xx}^{(\lambda ,\xi )}\left(T\right),$$

(53)

where ${\sigma }_{xx}^{(\lambda ,\xi )}\left(T\right)$ is given in Equation (48), including the vertex correction. Our results for $T=0$ are presented in

Table 4. Computed values for the total conductivity ${\sigma }_0=\sigma (T=0)$ at zero temperature for each material. We consider a value of $n_d={10}^{11} {\mathrm{cm}}^{-2}$.

Material	${\mathit{\sigma }}_0$ [${10}^3$ ${\mathbf{\Omega }}^{-1} {\mathrm{cm}}^{-1}$]
TaAs	1.5
TaP	7.9
NbAs	7.5
NbP	34.6

.

The conductivity as a function of temperature, for the transition metals’ monopnictides TaAs, TaP, NbAs and NbP, is presented in Figure 7 for all of them compared, and individually in the panel Figure 8.

Now, let us study the conductivity behavior with respect to the density of dislocations $n_d$. In Figure 9, we present a plot of the natural logarithm of the conductivity versus temperature for three different concentrations of dislocations.

The total conductivity as a function of the concentration of defects and at zero temperature is presented in Figure 10.

Finally, a plot of the resistance, defined as the inverse of conductivity, as a function of the dislocations’ density is presented in Figure 11.

5. Discussion and Conclusions

In this work, we have studied the effect of a distribution of mechanical defects, i.e., torsional dislocations, over the electrical conductivity of the family of transition metals monopnictides TaAs, TaP, NbAs and NbP. Our theory is based on the mathematical analysis of the scattering phase shifts from a single defect, as stated in our previous work [20,21,22,23,27]. We extended this previous analysis to develop a Green´s function formalism, in order to represent the scattering due to a finite concentration of randomly distributed defects. Within the non-crossing approximation for the self-energy, we solved explicitly for the disorder-averaged retarded Green´s function that allows us to calculate the electrical conductivity in the Kubo linear-response formalism. We obtained general analytical expressions in terms of the parameters involved in the low-energy model representing the family of materials, and using the ab-initio estimations for such parameters, we provided a characterization of the conductivity as a function of temperature and concentration of defects for the transition metal monopnictides TaAs, TaP, NbAs and NbP. As a universal feature, we identified a $\sim T^2$ temperature dependence for $T\gg \hslash v_F{k}_F/k_B$, where the pre-factor depends on material-specific microscopic parameters as well as in the concentration of dislocations $n_d$ through the scattering relaxation time. Our results do not involve the electron–phonon scattering effects that will presumably contribute at finite temperatures. However, those can be included via Mathiessen’s rule in an overall relaxation time combining Equation (50) for ${\tau }_{tr}$ with a separate theoretical estimation for the electron–phonon relaxation time ${\tau }_{e-ph}$, as follows: ${\tau }^{-1}={\tau }_{tr}^{-1}+{\tau }_{e-ph}^{-1}$. Since the electron–phonon interaction that determines the magnitude of ${\tau }_{e-ph}$ is an entirely different physical mechanism, it deserves an analysis on its own, to be communicated in a separate article which is under current development.