Electric field manipulation of surface states in topological semimetals

TL;DR

This paper analyzes how perpendicular electric fields influence surface states in topological semimetal thin films, revealing shifts in dispersion, Fermi velocity renormalization, and altered surface state decay, affecting quantum transport.

Contribution

It introduces an analytical perturbation theory approach to study electric field effects on surface states in Weyl and Dirac semimetals, validated by numerical calculations.

Findings

Electric fields shift dispersion relations.

Fermi velocity is renormalized by electric fields.

Surface state decay properties are modified.

Abstract

We investigate the consequences of applying electric fields perpendicularly to thin films of topological semimetals. In particular, we consider Weyl and Dirac semimetals in a configuration such that their surface Fermi arcs lie on opposite edges of the films. We develop an analytical approach based on perturbation theory and a single-surface approximation and we compare our analytical results with numerical calculations. The effect of the electric field on the dispersion is twofold: it shifts the dispersion relation and renormalizes the Fermi velocity, which would, in turn, have direct effects on quantum transport measurements. Additionally, it modifies the spatial decay properties of surface states which will impact the connection of the Fermi arcs in opposite sides of a narrow thin film.

Tables

Table 1. Table 1: Parameters for the Hamiltonian of Na 3 ​ Bi subscript Na 3 Bi \mathrm{Na_{3}Bi} extracted from Ref.  Wang  et al. , 2012 .

$c_0=-\mathrm{0.063\hspace{0.17em}82}{\mathrm{eV}}^2$	$m_0=-\mathrm{0.086\hspace{0.17em}86}{\mathrm{eV}}^2$
$c_1=8.7536\mathrm{eV}{\AA }^2$	$m_1=-10.6424\mathrm{eV}{\AA }^2$	$v=2.4598\mathrm{eV}\AA$
$c_2=-8.4008\mathrm{eV}\AA$	$m_2=-10.3610\mathrm{eV}\AA$

Equations

$${\mathcal{H}}_{\zeta }=(m_0-m_1{\mathbf{k}}^2){\sigma }_z+v{k}_z{\sigma }_x+\zeta v{k}_y{\sigma }_y,$$

$$\mathcal{H}_{D}({\mathbf{k}})=\left(\begin{array}[]{cc}\mathcal{H}_{\zeta=+1}({\mathbf{k}})&0\\ 0&\mathcal{H}_{\zeta=-1}({\mathbf{k}})\end{array}\right)~{}.$$

$$E_b=\pm \sqrt{{\left(m_0-m_1|\mathbf{k}{|}^2\right)}^2+v^2(k_y^2+k_z^2)}.$$

$${\psi }_s=\frac{A_s}{\sqrt{2}}\left(e^{-{\lambda }_1(\eta z+w)}-e^{-{\lambda }_2(\eta z+w)}\right)\left(\begin{array}{c}1\\ \eta i\end{array}\right),$$

$${\lambda }_1=\Delta +\sqrt{F},{\lambda }_2=\Delta -\sqrt{F}.$$

$$\Delta \equiv v/(2m_1),F\equiv k_x^2+k_y^2-R^2+{\Delta }^2,$$

$$k_x^2+k_y^2<R^2.$$

$$E_s=\eta \zeta v{k}_y.$$

$$r_{\mathrm{trans}}^2\equiv R^2-{\Delta }^2.$$

$${\mathcal{H}}_{\zeta }={ϵ}_0(\mathbf{k})1_2+M(\mathbf{k}){\sigma }_z+v(\zeta k_x{\sigma }_x-k_y{\sigma }_y),$$

$${\psi }_s$$

$${\Phi }_{\eta }$$

$$F$$

$$\Delta$$

$$E_s=\epsilon (k_z)+\eta \zeta v{C}_3{k}_x,$$

$$\epsilon (k_z)=C_1+C_2{k}_z^2.$$

$${\mathcal{H}}_f=efz{1}_2,$$

$$\delta E_s^1=⟨{\psi }_s^0|{\mathcal{H}}_f|{\psi }_s^0⟩,$$

$$\delta E_s^1=\eta ef\left[{\Gamma }_0+{\Gamma }_1(k_x,k_y)\right],$$

$${\Gamma }_0$$

$${\Gamma }_1(k_x,k_y)$$

$$E_s(k_x=0)=\eta {\gamma }_0+\eta \zeta v_f^{\text{1PT}}{\tilde{k}}_y+\mathcal{O}(f^2,k_y^2),$$

$$v_f^{\text{1PT}}\equiv v+ef{\left[\zeta \frac{\partial {\Gamma }_1(0,k_y)}{\partial k_y}\right]}_{k_{\mathrm{shift}}}.$$

$${\psi }_b^0=\sqrt{\frac{1}{L_x{L}_{y}w}}{e}^{i{q}_{x}x}{e}^{i{q}_{y}y}\sin \left[q_z(z+w)\right]{\Phi }_b,$$

$$\delta E_s^2=\sum _b\frac{|⟨{\psi }_s^0|H_f|{\psi }_b^0⟩{|}^2}{E_s^0-E_b^0},$$

$$\delta E_s^2=-e^2{f}^2\eta \zeta {\Gamma }_2(k_x,k_y),$$

$${\Gamma }_2(k_x,k_y)\equiv \sum _b{}^{{}^{\prime }}\frac{|I{|}^2}{w}\frac{2v{k}_y}{(E_b^0{)}^2-(E_s^0{)}^2},$$

$$I=A_s{\int }_{-w}^{w}z\sin \left[q_z(z+w)\right]\left(e^{-{\lambda }_1(\eta z+w)}-e^{-{\lambda }_2(\eta z+w)}\right)dz,$$

$$E_s(k_x=0)$$

$$+\eta ef\left[{\Gamma }_0+{\Gamma }_1(0,k_y)\right]+\mathcal{O}(f^3).$$

$$v_f^{\text{2PT}}\equiv v+{\left[\zeta ef\frac{\partial {\Gamma }_1(0,k_y)}{\partial k_y}-e^2{f}^2\frac{\partial {\Gamma }_2(0,k_y)}{\partial k_y}\right]}_{k_{\mathrm{shift}}}.$$

$$\delta E_s^1=ef\eta \left[{\Gamma }_0+{\Gamma }_1(k_x)\right],$$

Full text

I Introduction

Topological materials have attracted great interest in the last decade since they exhibit new fundamental phenomena and hold great promise for far-reaching technological applications. A hallmark of topological materials is quantized response functions and the existence of protected gapless surface states, which arise due to the nontrivial topology of the bulk states by virtue of the bulk-boundary correspondence Shen (2017). Nontrivial topology can be characterised by topological invariants according to the symmetry class of the system enabling a complete classification in both gapped and gapless systems Chiu et al. (2016). The gapped case was the first under study, starting the fruitful field of topological insulators (TIs) Bernevig (2013); Shen (2017). On the other hand, gapless systems assemble the family of topological semimetals (TSMs), where the valence and conduction bands only touch at a zero-measure set of points in the Brillouin zone (BZ). In particular, topological Weyl and Dirac semimetals (WSMs and DSMs) are three-dimensional phases of matter in which these isolated touching points, dubbed Weyl and Dirac points respectively, are protected by topology and symmetry Burkov (2018); Armitage et al. (2018).

In WSMs, near the Weyl node, the dispersion relation appears as a three-dimensional analogous to graphene and can be described by an anisotropic version of the Weyl equation. The low-energy quasiparticles behave, then, as relativistic Weyl fermions. The conduction and valence bands are individually non-degenerate and hence either time reversal symmetry or inversion symmetry need to be broken. The Weyl nodes are monopoles of Berry curvature and the charge associated with them is called chirality Burkov (2018); Armitage et al. (2018). Due to the Nielsen-Ninomiya no-go theorem Nielsen and Ninomiya (1983), Weyl nodes always come in pairs and can annihilate only in pairs. Hence, the robustness of the Weyl nodes is quantified by the separation of the nodes in reciprocal space. For a given momentum between the Weyl nodes, the 3D WSM can be mapped onto a 2D TI Shen (2017); Turner and Vishwanath (2013). This correspondence leads to non trivial Chern numbers and to protected states localized at the surface. The protected states, called Fermi arcs, lie in open contours at the Fermi energy. Several materials have been predicted and confirmed to be topological WSMs with Fermi arcs by means of ARPES experiments. The most representative is the $\mathrm{TX}$ family where $\mathrm{T=Ta/Nb}$ and $\mathrm{X=As/P}$, $\mathrm{TaAs}$ being a particular case Huang et al. (2015); Weng et al. (2015); Xu et al. (2015); Lv et al. (2015). These materials belong to the so-called Type I WSMs and are characterised by a discrete point-like Fermi surface. However, it is possible to have anisotropies in the dispersion so that the Fermi surface is not a point but it becomes an open surface, thereby leading to electron and hole pockets. Materials within this category are referred to as Type II Weyl semimetals Huang et al. (2016).

The DSMs are obtained when both time reversal and inversion symmetry are present. Consequently, the Dirac points have a four-fold degeneracy and the net Chern number of the nodes is zero. Each Dirac point can be constructed by imposing two Weyl nodes with opposite chirality and, in order to be topologically protected, it must be stabilised by additional symmetries as the up-down parity symmetry Gorbar et al. (2015a, b) or by space-group symmetries Young et al. (2012); Gibson et al. (2015); Steinberg et al. (2014). This is the case of compounds like $\mathrm{A_{3}Bi}$ where $\mathrm{A=Na,K,Rb}$ and $\mathrm{Cd_{3}As_{2}}$, in which ARPES experiments ratify the existence of the Fermi arcs Neupane et al. (2014); Liu et al. (2014); Zhang et al. (2014); Xu et al. (2015). The manufacturing of high quality thin films Hellerstedt et al. (2016); Zhang et al. (2014) and ultrathin films of $\mathrm{Na_{3}Bi}$ Collins et al. (2018) has made this material one of the most promising candidates for technological applications of topological properties. For example, monolayer and bilayer films of $\mathrm{Na_{3}Bi}$ have bulk bandgaps greater than $300\text{\,}\mathrm{meV}$, suggesting that topological properties in these thin films will survive at room temperature Niu et al. (2017); Collins et al. (2018).

The topological response of WSMs and DSMs comprises the manifestation of the chiral anomaly in a large negative magnetoresistance, in the presence of both electric and magnetic fields, and in an anomalous Hall effect (in WSMs), due to the transport of the surface states.  Shen (2017); Armitage et al. (2018) Therefore, the renormalization of the Hamiltonian parameters due to time-dependent external fields has been the subject of intense research. Narayan (2015); Chan et al. (2016) In this case, new surface states with interesting properties may appear. González and Molina (2016) Overall, the understanding of the effect of external fields on the topological phases and transport phenomena is a field of great interest both from a first principles standpoint and from the perspective of possible applications. Hence, a special interest resides in understanding the effect of the external fields in the most direct manifestation of topology, i.e. the surface states.

The present manuscript studies the effects of an electric field applied perpendicularly to the surface and to the direction joining the line of nodes. In Sec. II, we introduce the theoretical framework: a minimal model for a generic WSM and a low-energy effective model for the DSM $\mathrm{Na_{3}Bi}$. Section III describes the effect of an external electric field applied perpendicularly to the direction along which the two nodes are aligned. We provide analytic results obtained by means of perturbation-theory techniques. In order to extend our results to non-perturbative regimes, we provide numerical calculations that match the analytic results. Finally, in Sec. V we finish with some conclusions and a brief analysis of the experimental feasibility for measuring the predicted behavior in thin films of the topological semimetal $\mathrm{Na_{3}Bi}$.

II Model Hamiltonians

II.1 Minimal model for a WSM and DSM

The minimal setting for describing a Weyl semimetal consists of two Weyl points that are realised in a time-reversal-symmetry-breaking scenario while preserving inversion symmetry Armitage et al. (2018); Shen (2017). A generic low-energy Hamiltonian that meets these requirements can be written as González and Molina (2017)

[TABLE]

where $\zeta$ is $\pm 1$ depending on the chirality, $\sigma_{i}$ with $i=x,y$ and $z$ are the Pauli matrices and ${\mathbf{k}}=(k_{x},k_{y},k_{z})$ are momentum operators. In addition, $m_{i}$ with $i=\{0,1\}$ account for mass parameters and $v$ is the Fermi velocity. From the previous Hamiltonian, a DSM can be build from two copies with opposite chirality, which are time-reversal partners. Thus, if no chirality-mixing term is considered, the Hamiltonian can be written in the following block diagonal form

[TABLE]

In bulk, the dispersion relation of Hamiltonian (1) is given by

[TABLE]

The valence band and the conduction band touch at the aforementioned Weyl points, located at ${\mathbf{k_{W}}}_{\pm}=(\pm\sqrt{m_{0}/m_{1}},0,0)$. The Weyl points are monopoles of Berry curvature and have Chern number equivalent to their chirality. The nonzero Chern number leads, according to the bulk-boundary correspondence, to surface states named Fermi arcs Hatsugai (1993).

For the analytic approach, we consider semi-infinite geometries in the perpendicular direction to the node separation by introducing a single surface in the $Z$-direction. In order to explore how the location of the boundary affects the dispersion, we shall let the boundary sit at $z=-\eta w$, where $\eta=\pm 1$ indicates the position of the surface with respect to the plane $z=0$ and $w>0$. Bear in mind that we are considering a single surface but are allowing for both signs of $\eta$ to see its manifestation in the dispersion. Surface states are obtained from Hamiltonian (1) by using the ansatz $\psi_{s}\sim e^{i{\mathbf{k}}_{\bot}\cdot{\mathbf{r}}}e^{-\lambda(\eta z+w)}\Phi$, where $\Phi$ is a constant spinor, ${\mathbf{k}}_{\bot}=(k_{x},k_{y},0)$ and $\lambda$ is a complex number with a nonzero real part. For simplicity, from now on we make implicit the plane-wave dependence $\exp(i{\mathbf{k}}_{\bot}\cdot{\mathbf{r}})$. A solution satisfying Dirichlet boundary conditions is given by González and Molina (2017)

[TABLE]

where $A_{s}$ is a normalization factor and

[TABLE]

Here we have defined

[TABLE]

with $R\equiv\sqrt{m_{0}/m_{1}}$. Surface states occur whenever $\Re[\lambda_{1,2}]>0$, which implies $F<\Delta^{2}$. Equivalently, surface states are restricted to a circle of radius $R$ in momentum space

[TABLE]

At a given $k_{x}$ that supports surface states, the dispersion relation is linear in $k_{y}$ and depends on $\zeta$ and on the position of the surface

[TABLE]

This dispersion is represented in Fig. 1 together with the bulk bands. From the dispersion relation we can give another interpretation to the condition of existence of surface states. Indeed, the region (7) defines a circle in the plane $(k_{x},k_{y})$ outside which the surface states bands become tangent to the bulk states dispersion. This leads to a hybridisation between surface states and bulk states which prevents us to use the ansatz and to have Fermi arcs.

It is interesting to point out that there are two types of surface states, depending on whether $\lambda_{1,2}$ have or not an imaginary part. If they do, the exponential decay will be accompanied by oscillations. From equation (6) we can define a transition radius as

[TABLE]

Then, we can see that oscillatory states live in a circle of radius $r_{\mathrm{trans}}$, whereas purely exponential states are found in a planar ring of inner radius $r_{\mathrm{trans}}$ and outer radius $R$. Hereafter we shall denote the oscillatory states as type A and purely exponential states as type B. Since the type B states are closer to the bulk states in energy, these have longer decay lengths than type A Benito-Matías and Molina (2019); González and Molina (2017). Notice that $r_{\mathrm{trans}}^{2}$ can either be positive or negative, depending on the parameters of the model. If positive, both type A and B states exist. However, if negative, only type B arise. In the following, we will loosely dub Weyl semimetals hosting type A states as type A Weyl semimetals, carefully remembering that these semimetals also host type B states.

II.2 Model for the Dirac semimetal $\mathrm{Na_{3}Bi}$

In order to elucidate the generality of our results and to have a more realistic approach to the problem, we will also study the case of a low-energy effective Hamiltonian that describes $\mathrm{A_{3}Bi}$ ($\mathrm{A=Na,K,Rb}$) Wang et al. (2012) and $\mathrm{Cd_{3}As_{2}}$ Wang et al. (2013) around the $\Gamma$ point. These compounds have a single band inversion occurring near the $\Gamma$ point that has been observed by ARPES measurements Neupane et al. (2014); Liu et al. (2014); Zhang et al. (2014); Xu et al. (2015). For concreteness, we restrict our analysis to the case of $\mathrm{Na_{3}Bi}$. After a density functional theory analysis, the Hamiltonian can be cast in the form of a DSM equivalent to Eq. (2) in which the corresponding WSM Hamiltonian is replaced by Wang et al. (2012)

[TABLE]

where $\epsilon_{0}({\mathbf{k}})=c_{0}+c_{1}k_{z}^{2}+c_{2}(k_{x}^{2}+k_{y}^{2})$ and $M({\mathbf{k}})=m_{0}-m_{1}k_{z}^{2}-m_{2}(k_{x}^{2}+k_{y}^{2})$, $c_{i}$ and $m_{i}$ with $i=\{0,1,2\}$ being constant parameters and mass terms respectively. Notice that, except for the diagonal contribution, $\mathcal{H}_{\zeta}$ is a rotated version of (1) that allows for anisotropy along the $Z$-direction. Threefold rotational symmetry implies that chirality-mixing terms are of order $\mathcal{O}(k^{3})$ and we will neglect them in our analysis, thereby effectively decoupling the two Dirac nodes Wang et al. (2012). This type of Dirac semimetal is also referred to as a $\mathbb{Z}_{2}$ WSM Gorbar et al. (2015a, b).

The surface states of the model have been worked out in detail in Ref. Benito-Matías and Molina, 2019. By placing a surface termination at $y=-\eta w$ and imposing Dirichlet boundary conditions as before, the surface states take the form

[TABLE]

where $A_{s}$ is a normalisation factor. Here, $\lambda_{1,2}$ are defined in equation (5), where $\Delta$ and $F$ are now given by

[TABLE]

being $k_{x,0}\equiv{c_{2}}\Delta/{m_{2}}$ and $R^{2}\equiv{m_{0}}/{m_{2}}+\Delta^{2}\left({c_{2}}/{m_{2}}\right)^{2}$. In contrast to the previous case, the diagonal term $\epsilon_{0}({\mathbf{k}})\mathbb{1}_{2}$ leads to a dispersion that is no longer flat along the $Z$-direction (recall that this model is rotated with respect to the minimal model presented above). Instead, surface states now have the following dispersion

[TABLE]

where

[TABLE]

In these two equations, $C_{1},C_{2}$ and $C_{3}$ are a combination of the Hamiltonian parameters and are given by $C_{1}=c_{0}+c_{2}{m_{0}}/{m_{2}}$, $C_{2}=c_{1}-c_{2}{m_{1}}/{m_{2}}$ and $C_{3}=\sqrt{1-c_{2}^{2}/m_{2}^{2}}$. Notice that if we set all $c_{i}=0$ and $m_{1}=m_{2}$, all expressions reduce to those obtained in the previous section and we recover a flat band behaviour along the $Z$-direction. The classification in type A and type B states remains the same, that is, type B states are purely exponential whereas type A have an oscillatory component to the exponential decay. In $\mathrm{Na_{3}Bi}$ only type B states arise as can be demonstrated by introducing into the above definitions the parameters listed in table 1.

III Electric field

In this section we study the effect of an external electric field applied perpendicularly to the surface. The surface configurations correspond to those studied in the previous section. Earlier works in Dirac materials have proven that a renormalization of the Fermi velocity of the surface states occurs in the presence of a perpendicular electric field Díaz-Fernández A. et al. (2017); Díaz-Fernández and Domínguez-Adame (2017). Therefore, we expect similar effects in WSMs and DSMs. We approach the problem in two ways: i) by way of perturbation theory (PT) to obtain analytic results and ii) by performing numerical calculations based on the Python package Kwant Groth et al. (2014). In this way, we can compare the validity of the analytic results and study the non-perturbative regime.

III.1 Minimal model for a WSM and DSM

We begin by considering the minimal model under the influence of an electric field applied along the $Z$-direction, perpendicular to the surface. Thus, its perturbation reads

[TABLE]

where $e$ is the elementary electric charge and $f$ is the external electric field. The first order correction in perturbation theory is given by

[TABLE]

where $\psi_{s}^{0}$ are the surface states in the absence of an electric field as defined in Eq. (4). From Eq. (16), the correction to the energy is given by

[TABLE]

where we have defined

[TABLE]

Notice that these corrections are independent of chirality and of the surface state type. Equation (17) reveals that there is a constant shift in energy, together with a momentum dependent term due to the perturbation. The reasoning behind this is that the decay lengths are momentum dependent as well. In fact, we can try to provide with a simple argument. Let us set zero potential at $z=0$ such that it is $-ef\eta w$ at the surface termination. Then, one may consider a surface far from $z=0$, so that $\Gamma_{1}$ becomes negligible with respect to $\Gamma_{0}$, so long as we consider low momenta. Moreover, all terms in $\Gamma_{0}$ are negligible except for $-w$ if the surface is sufficiently far from $z=0$. In that case, the correction is simply $-ef\eta w$. This means that the potential is locally acting on the surface states by simply lowering their energy by an amount equal to the value of the potential at the surface. More generally, in first order perturbation theory we are calculating $ef\langle z\rangle$, which essentially amounts to calculating the expectation value of the position in the unperturbed surface state. Hence, both $\Gamma_{1}$ and the terms in $\Gamma_{0}$ that accompany $w$ represent a correction to $-\eta w$ due to the fact that the surface state has some extension and penetrates slightly into the bulk. Notice that $\Gamma_{1}$ presents a first order pole at the momenta located in the circle of radius $R$. This is consistent with the fact that upon approaching the edge of the circle, surface states become less and less localized, till they merge with the extended states of the bulk and the uncertainty in position becomes infinite.

Having said that, we can proceed to study the velocity renormalization, which occurs due to the squared momenta in $\Gamma_{1}$. That is, if we consider low momenta, then there is no velocity renormalization to first order, similarly to what has been observed in Refs. Díaz-Fernández A. et al., 2017; Díaz-Fernández and Domínguez-Adame, 2017. Let us focus on the dispersion relation in the $k_{x}=0$ plane, where the dispersion is linear in the absence of electric field. The first order PT gives a term that typically increases the velocity if (7) is fulfilled (being the only range in which the expression is valid as we will discuss in the following subsections). The bands of the surface states are displaced, together with the bulk bands, due to the electric field. We shall denote the new position of the zero energy surface states as $k_{\mathrm{shift}}$ along the $k_{y}$ direction. With this definition, the dispersion relation within the first order PT is now

[TABLE]

where $\widetilde{k}_{y}\equiv k_{y}-k_{\mathrm{shift}}$. In addition $\gamma_{0}$ is a constant factor and the renormalized velocity reads

[TABLE]

Notice that the velocity renormalization, up to first order PT, does not depend neither on $\zeta$ nor on $\eta$ because $k_{\mathrm{shift}}$ depends explicitly on chirality as $k_{\mathrm{shift}}=\zeta{\lvert}k_{\mathrm{shift}}{\rvert}$.

Next, we obtain the second order PT correction in a slab extended from $-w\leq z\leq w$ but infinite in the $X$ and $Y$ directions. Furthermore, we assume that the width $w$ is large enough to use the surface states obtained in the semi-infinite approximation. That is, the width is larger than the decay length of the surface states into the bulk so that states of opposite surfaces cannot hybridize. In a finite slab, the bulk states have the dispersion relation (3) with a quantized momentum in the $Z$ direction due to the finite size of the system. The general form of the bulk states of the unperturbed Hamiltonian is

[TABLE]

where $q_{z}={n\pi}/({2w})$ with $n\in\mathbb{Z}$ and $\Phi_{b}$ is a constant and normalised spinor. The second order PT is given by

[TABLE]

where the sum index runs over all bulk states and $E_{b}^{0}$ ($E_{s}^{0}$) is defined according to Eq. 3 (Eq. 8). By implementing (22) the second order correction obtained is

[TABLE]

with

[TABLE]

where the primed sum runs over positive-energy bulk states due to the symmetry of the energy spectrum and

[TABLE]

only depends on the system parameters and $q_{z}$, but not on $\eta$. Notice that ${\lvert}E_{s}^{0}{\rvert}<{\lvert}E_{b}^{0}{\rvert}$ and therefore the second order correction $\delta E_{s}^{2}$ is always negative. Hence, the corrected energy dispersion up second order in the plane $k_{x}=0$, is given by

[TABLE]

The contribution of the second order correction does not introduce new relevant effects except for changing $k_{\mathrm{shift}}$ and a introducing a reduction of the velocity as shown in the following expression

[TABLE]

In summary, we identify two main effects of the electric field in the surface states: the shifting of the momenta of the surface states and the renormalization of the velocity. Since the type of surface state (A or B) depends on these two parameters, we expect that the electric field may induce a transition between different types of states. In the following, we discuss the main features that arise in the system under these conditions by comparing analytic treatment with the numerical calculations.

The shifting of the momenta of the surface states is the most salient result of the application of an external electric field. We refer to this effect as the shifting of the cone of the dispersion relation of the surface states. Notice that the cone we are referring to is not a Dirac cone. Rather, it is the cone formed by each surface contributing states with opposite velocities in the field-free case. In a slab system, in the already mentioned limit of two decoupled surfaces, the accordance between the simulation results and the analytic ones is expected to be good as long as PT is valid. In a gapless system, an energy scale for PT is not so clear as in a gapped case. Moreover, the analytic expression of the surface states (4) is valid only for momenta that fulfil Eq. (7). Hence, we expect PT to fail for electric fields that shift the momenta to values of ${k_{x}^{2}+k_{y}^{2}}\simeq R^{2}$. In the following we refer to these momenta as the critical momenta. PT results reflect intrinsically this range of validity. In fact, $\Gamma_{1}$ and $\Gamma_{2}$ have respectively a first and a second order pole at the critical momenta. From this perspective, it is clear that second order PT is going to fail faster than first order and the shift given by the constant $\Gamma_{0}$ will have a wider convergence radius.

Despite the aforementioned limitations, an expansion parameter for the PT can be defined. If we restrict the momenta of the surface states to $k_{x}=0$ and $k_{y}\ll R$, we can consider as the PT parameter, the ratio between the electric potential at the surfaces, $-ef\eta w$, and the gap between the surface and the bulk states at the critical momenta, $v\sqrt{m_{0}/m_{1}}$. In fact, it quantifies an effective gap compared with the displacement of the bands due to the electric potential. Hence, PT will be valid if ${\lvert}efw/(vR){\rvert}\ll 1$. For the parameters in Fig. 8, we found that $f\ll$12~{}(47)\text{,}\mathrm{m}\mathrm{V}\mathrm{n}\mathrm{m}^{-1} for $v=$1~{}(4)\text{\,}\mathrm{eV}\text{\,}\mathrm{nm}, respectively. In the previously mentioned figure, we have plotted the 1PT up to values of 0.5 of the expansion parameter, while 2PT is plotted up to values of 0.1 finding a good accordance with the simulation results.

As already mentioned, in order to quantify the shifting we define the shifted momenta $k_{\mathrm{shift}}$ as the momenta at which the branches intersect in the plane of $k_{x}=0$. It is obtained by finding the intersection of the two energy branches of the surface state dispersion (III.1) at a fixed chirality. Figure 2 shows a comparison of the PT result with the simulations: the second order starts to fail for really small electric fields and does not introduce relevant corrections, whereas the first order PT reproduces very well the simulation results. Therefore we will neglect the second order corrections in the following.

Another aspect predicted by PT is the velocity renormalization. Since $v_{f}^{\mathrm{1PT}}$ is obtained as a derivative of $\Gamma_{1}$, it has a second order pole for the critical momenta (see Eq. (20)). Thus, the first order PT for the velocity converges within the same radius as the second order PT for $k_{\mathrm{shift}}$. In the minimal model, the radius of convergence of second order pole functions is utterly restricted to small $f$ (see Fig. 2), where a significant velocity renormalization is absent. Therefore, the effect can be only numerically studied in the non-perturbative regime. Figure 3 shows the evolution of the velocity and the $k_{\mathrm{shift}}$ with the electric field. Even if the PT is not valid for the regime of electric fields studied, it is worth mentioning that it captures important features of the effect as the independence of the results from the type of surface states.

The renormalization of the velocity and the shifting of the cone vertex leads to the possibility of a phase transition between types of surface states. In fact, the change in the velocity is directly related to the type transition because it modifies $r^{2}_{\mathrm{trans}}$ [see Eqs. (6) and (9)]. Moreover, the shifting of the momenta of the surface states modifies the possible momenta that can be spanned in the left-hand side of (7), typically increasing the terms on that side of the equation. To study the type transition, an accurate comparison of the evolution of $k_{\mathrm{shift}}^{2}$ and $r^{2}_{\mathrm{trans}}$ as a function of $f$ is needed. From the results exposed previously, the transition from type A to type B is expected: the shifting of the momenta is the most appreciable effect of the electric field, more noticeable than the velocity renormalization. Figure 4 shows the transition $A\to B$ and the conversion of the oscillatory decay into a purely exponential one.

The simulations confirm that the transitions $B\to A$ do not take place. To prove it, we simulate a type B WSM with $v=v_{\mathrm{lim}}$, where $v_{\mathrm{lim}}\equiv 2\sqrt{m_{0}m_{1}}$, corresponding to $r^{2}_{\mathrm{trans}}=0$. In this regime, a decreasing of the velocity turns the $r^{2}_{\mathrm{trans}}<0$ and allows for a type transition as long as $k_{\mathrm{shift}}^{2}$ does not increase beyond $r_{\mathrm{trans}}^{2}$. In Fig. 5 the evolution of $k_{\mathrm{shift}}^{2}$ and $r_{\mathrm{trans}}^{2}$ are compared for two widths: the effect increases with $w$, as expected from the PT results, but the general behavior does not change with the thickness.

To complete our analysis of the minimal models, we study the case of the DSM with Hamiltonian (2) comprising two copies of the minimal Weyl Hamiltonian (1) with different chiralities in each copy. In the absence of electric field, each boundary hosts two surface states with linear dispersion and opposite chirality. The electric field does not mix chiralities, but it breaks the spatial inversion symmetry leading to a splitting of the two cones of the opposite surfaces. In fact, the sign of the momentum shift induced by the electric field on the dispersion of the Fermi arcs depends on the chirality and $\eta$ (see Fig. 6). The cited figure shows also the bands as a function of $k_{x}$ in the plane $k_{y}=0$: the effect of the electric field in the dispersion is the coalescence and distortion of the bands. This result is expected as the crossing of the cones is shifted to non-zero $k_{y}$. It is important to notice that these cones are not the Dirac cones, rather they correspond to the intersection of the two surface states of the same chirality and opposite surfaces near the zero energy, as mentioned earlier in the text. In fact, the actual Dirac cones actually shift up or down in energy, depending on the surface they are located at.

Before concluding this section, it is interesting to notice that, as it can be seen in Fig. 6, the surface states linear branches directly enter into the bulk bands and remain unaltered despite the presence of the latter. This means that surface states and bulk states do not hybridize. We have observed this fact also by looking at the evolution of the surface states as the field is increased. To this respect, more details can be found in the Supplemental Material Sup . It is important to stress, however, that the bulk states are modified by the electric field as well. States of the conduction (valence) band transition from being completely extended through the system to localize towards regions of lower (higher) electric potential. This is in fact the reason why bulk and surface states do not hybridize, because the former localize at regions where the latter have negligible probability density.

III.2 Model for the Dirac semimetal $\mathrm{Na_{3}Bi}$

Throughout this section, we extend the previous results of the minimal model to the Hamiltonian for $\mathrm{Na_{3}Bi}$. Since the electric field does not mix the chiralities, we study the $2\times 2$ model with nondegenerate Weyl nodes given by Eq. (10). We follow the same approach as with the minimal model up to first order PT applying Eq. (16) to the surface states defined by (11a). In order to follow the shifting of the cones we set $k_{z}=0$ and to simplify the notation, we define all the quantities at zero $k_{z}$. The first order PT correction now is

[TABLE]

with $\Gamma_{0}$ a constant obtained from the model parameters and $\Gamma_{1}(k_{x})$ a function of $k_{x}$

[TABLE]

Before proceeding further, a couple of words must be said about the convergence of PT in this case. As for the minimal model, the surface states (11a) exist if $F<\Delta^{2}$ Benito-Matías and Molina (2019). The former condition implies now that

[TABLE]

These domains represent ellipses in the $k_{x}-k_{z}$ plane and imposing $k_{z}=0$ they are a constraint for the values of $k_{x}$ that must fulfil $k_{x,c}^{-}<{k_{x}}<k_{x,c}^{+}$, where $k_{x,c}^{\pm}$ is a function of the model parameters and is obtained from (29). We expect PT to converge in the range $k_{x,c}^{-}\ll{k_{x}}\ll k_{x,c}^{+}$.

In order to further simplify the problem and have an intuitive idea of the phenomena, it is interesting to study the behavior of the energy dispersion for small momenta by expanding the perturbed energy in $k_{x}\ll 1$ as follows

[TABLE]

where $\gamma_{i}~{}(i=0,1,2)$ are positive constants. Moreover the renormalized velocity is defined as previously by way of an expansion near the shifted cones as follows:

[TABLE]

From inspection of Eq. (30) we observe some features revealed by the minimal model: the bands are displaced by a constant term and they are transformed by the terms $\gamma_{1}$ and $\gamma_{2}$. However, this scenario presents important modifications with respect to the minimal model. First of all, $\Gamma_{1}$ depends on the chirality $\zeta$ and $\eta$, therefore the renormalized velocity is asymmetric. Using the numerical values of table 1 and the PT expression for $v^{\text{1PT}}_{f}$ (31), we find that the velocity increases for surface states with $\eta\zeta=+1$ and decreases for $\eta\zeta=-1$. This phenomenon is perfectly observed also in the simulations, as seen in Fig. 7.

Moreover, the bands of the surface states are not flat and they do not cross at zero energy for $k_{z}=0$. Therefore, the shifting produced by the electric field modifies the crossing of the two branches in momenta and in energy what defines a new $k_{\mathrm{shift}}$ and $E_{\mathrm{shift}}$. The evolution of these two quantities is plotted in Fig. 8 along with the velocity renormalization. The previously mentioned figure compares PT with the numerical simulations finding a good accordance as long as the shifted momenta fulfill the already discussed restrictions. We find that the PT converges if $f\lesssim 10^{-2}~{}$\mathrm{m}\mathrm{V}\mathrm{\SIUnitSymbolAngstrom}^{-1} for systems with $w\sim 10^{2}~{}$\mathrm{\SIUnitSymbolAngstrom}. Remarkably, despite considering weak electric fields, the effect on the velocity renormalization is not negligible at all.

Finally, let us pay attention to the dispersion relation for the complete Hamiltonian, i.e. the DSM Hamiltonian, as plotted in Fig. 9. The effect is analogous to the one found in the minimal model, in that the branches in black move away from each other in momentum as the bulk states approach in energy.

IV Experimental proposal

A very recent paper has reported experimental measurements in the presence of an electric field in ultra thin $\mathrm{Na_{3}Bi}$ films Collins et al. (2018). The cited article, based on a theoretical proposal Pan et al. (2015), studies a different regime from the present work: the system comprises few-layer $\mathrm{Na_{3}Bi(001)}$ films grown in the $Z$ direction, along which the electric field is applied. Its main result is the closing and opening of a gap such that topological and trivial phases are induced. The importance of this experimental realisation is relevant from the standpoint of applications of the present work. In fact, it assures the possibility of realization of the proposed scenario and gives a magnitude of the electric fields that can be implemented experimentally. In the ultrathin film set up, the electric fields were implemented using two methods: doping the surface with potassium and with scanning tunnelling spectroscopy varying the tip–sample separation. The obtained electric fields are of the order of a few $\mathrm{V}\mathrm{n}\mathrm{m}^{-1}$ in a sample of atomic width. In the system proposed in this work, the thickness of the sample must be of the order of $10^{2}$\mathrm{\SIUnitSymbolAngstrom}$$ in order to have two decoupled surfaces. In the literature these samples are labelled as thin films and have been grown by Molecular Beam Epitaxy for $\mathrm{Na_{3}Bi}$ Hellerstedt et al. (2016); Zhang et al. (2014) and $\mathrm{Cd_{3}As_{2}}$ Schumann et al. (2018). Due to this estimated thickness, we expect that lower fields can be achieved. Indeed, even for electric fields one order of magnitude lower than the ones obtained by  Collins et al. Collins et al. (2018), we find that the velocity renormalization is a non-negligible effect, as seen in Fig. 10.

V Conclusions

In summary, we have obtained a suitable method to displace the Dirac nodes of the surface states in the BZ with a tunable external electric field. Not only the position of the cones is modified but also the Fermi velocity can be altered. Moreover, both effects depend on the chirality of the node at a given surface. The renormalization of the Fermi velocity would have a direct impact on the electronic transport properties through the surface of the semimetal thin films. In the case of $\mathrm{Na_{3}Bi}$ we show that these effects would be quite significant and due to their chiral dependence, we envision the possibilities for applications to chiral electronic devices Kharzeev and Yee (2013).

For some range of parameters, the renormalization of the Fermi velocity induced by the external electric field implies a transition from type A, with oscillatory decay and very short decay lengths, to type B surface states, with longer decay lengths and pure exponential decay into the bulk. This may lead to a very large increase in the hybridization of the two surfaces in a thin film. The coupling of the opposite surfaces is a necessary ingredient of the recently observed 3D Quantum Hall effect based on Weyl orbits Zhang et al. (2019); Nishihaya et al. (2019) and we foresee important observable consequences of this transition in quantum Hall transport experiments with topological semimetals.

Acknowledgements.

We acknowledge financial support through Spanish grants PGC2018-094180-B-I00 (MCIU/AEI/FEDER, EU), FIS2015-63770-P and MAT2016-75955 (MINECO/FEDER, EU), CAM/FEDER Project No.S2018/TCS-4342 (QUITEMAD-CM) and CSIC Research Platform PTI-001.