RKKY coupling in Weyl semimetal thin films

TL;DR

This paper investigates the RKKY interaction between impurity spins on the surfaces of Weyl semimetal thin films, revealing anisotropic surface coupling and thickness-dependent bulk-mediated interactions.

Contribution

It provides analytical and numerical analysis of surface and bulk-mediated RKKY coupling in Weyl semimetal thin films, highlighting the strong, thickness-dependent interactions.

Findings

Surface spin coupling reflects Fermi arc anisotropy and spin-momentum locking.

Opposite surface coupling is mediated by both Fermi arcs and bulk states.

Coupling strength varies with film thickness, peaking at an optimal thickness.

Abstract

We consider the effective coupling between impurity spins on surfaces of a thin-film Weyl semimetal within Ruderman-Kittel-Kasuya-Yoshida (RKKY) theory. If the spins are on the same surface, their coupling reflects the anisotropy and the spin-momentum locking of the Fermi arcs. By contrast when the spins are on opposite surfaces, their coupling is mediated by the Fermi arcs as well as by bulk states. In this case the coupling is both surprisingly strong and strongly thickness dependent, with a maximum at an optimum thickness. We demonstrate our results using analytical solutions of states in the thin-film geometry, as well using a two-surface recursive Green's function analysis of the tight-binding model.

Tables

Table 1. Table 1: Assuming the couplings J i ​ i subscript 𝐽 𝑖 𝑖 J_{ii} between the impurity spins on the same surface go as ∼ R − α similar-to absent superscript 𝑅 𝛼 \sim R^{-\alpha} , the best fitted value of α 𝛼 \alpha is quoted for the result of Fig.  3 (b) for R / a 𝑅 𝑎 R/a between 50 and 80.

	$N_z$	$\alpha$
$J_{xx}$	9	$\approx$ 2.53
$J_{xx}$	17	$\approx$ 2.29
$J_{xx}$	33	$\approx$ 2.26

Table 2. Table 2: Assuming the couplings J i ​ i subscript 𝐽 𝑖 𝑖 J_{ii} between the impurity spins on two opposite surfaces go as ∼ ( R 2 + L Z 2 ) − α / 2 similar-to absent superscript superscript 𝑅 2 superscript subscript 𝐿 𝑍 2 𝛼 2 \sim(R^{2}+L_{Z}^{2})^{-\alpha/2} , the best fitted value of α 𝛼 \alpha is quoted in the left-most column for the result of Fig.  4 (b) for R / a 𝑅 𝑎 R/a between 50 and 80. For fixed L z subscript 𝐿 𝑧 L_{z} , the same can be fitted with R − β superscript 𝑅 𝛽 R^{-\beta} , which is shown in the middle column. For fixed R 𝑅 R , the results from Fig.  5 (b) can be fitted (for N z subscript 𝑁 𝑧 N_{z} between 13 and 21) with L z − γ superscript subscript 𝐿 𝑧 𝛾 L_{z}^{-\gamma} , which is shown in the right-most column.

$N_z$	$\alpha$	$\beta$
9	$\approx$ 3.28	$\approx$ 2.36
17	$\approx$ 3.33	$\approx$ 2.39
33	$\approx$ 3.24	$\approx$ 2.32

Equations

$$H(k)=\lambda \sum _{\alpha =x,y,z}{\sigma }^{\alpha }\sin k_{\alpha }-\mu +{\tau }^y{\sigma }^y{M}_k.$$

$$\displaystyle U=\frac{1}{2}\left(\begin{array}[]{cccc}-1&-i&-i&1\\ 1&-i&i&1\\ 1&i&-i&1\\ -1&i&i&1\end{array}\right),$$

$$H_{\eta }^{\prime }$$

$$+\eta {\sigma }_z(2+m-\cos k_x-\cos k_z)-\lambda {\sigma }_z\sin k_y.$$

$$H_{\eta }^{\prime }$$

$$H_{\eta \xi }^{\mathrm{low}}$$

$$\displaystyle\psi_{\rm vac}\propto\left(\begin{array}[]{c}\lambda(q_{z}+iq_{x})\\ M_{\eta}^{{\rm vac}}-E\end{array}\right)e^{i(q_{z}z+q_{x}x)},$$

$$q_z=i\kappa ,\mathrm{for}z\ge L_z,$$

$$q_z=-i\kappa ,\mathrm{for}z\le 0,$$

$$\displaystyle\psi_{>}\propto\left(\begin{array}[]{c}im_{0}+iq_{x}\\ \eta m_{0}-E\end{array}\right)e^{-m_{0}z}\approx\left(\begin{array}[]{c}i\\ \eta\end{array}\right)e^{-m_{0}z}.$$

$$\displaystyle\psi_{<}\propto\left(\begin{array}[]{c}-im_{0}+iq_{x}\\ \eta m_{0}-E\end{array}\right)e^{m_{0}z}\approx\left(\begin{array}[]{c}i\\ -\eta\end{array}\right)e^{m_{0}z}.$$

$$\psi (z=0)\propto {\psi }_{<}(z=0)\mathrm{and}\psi (z=L)\propto {\psi }_{>}(z=L_z),$$

$$\psi (z)=$$

$$\displaystyle+b\left(\begin{array}[]{c}\lambda(-q_{z}+iq_{x})\\ M_{\eta}(k_{y})-E\end{array}\right)e^{-iq_{z}z},$$

$$\displaystyle{\rm Det}\left(\begin{array}[]{cccc}i&\lambda(q_{z}+iq_{x})&\lambda(-q_{z}+iq_{x})&0\\ -\eta&M_{\eta}(k_{y})-E&M_{\eta}(k_{y})-E&0\\ 0&\lambda(q_{z}+iq_{x})e^{iq_{z}L_{z}}&\lambda(-q_{z}+iq_{x})e^{-iq_{z}L_{z}}&i\\ 0&(M_{\eta}(k_{y})-E)e^{iq_{z}L_{z}}&(M_{\eta}(k_{y})-E)e^{-iq_{z}L_{z}}&\eta\end{array}\right)=0.$$

$$\frac{\tanh \left(L_z\sqrt{(M_{\eta }/\lambda {)}^2-\chi }\right)}{L_z\sqrt{(M_{\eta }/\lambda {)}^2-\chi }}=-\frac{\lambda }{L_z\eta M_{\eta }},$$

$$\frac{\tanh (\sqrt{M_{\eta }^{\prime 2}-{\chi }^{\prime }})}{\sqrt{M_{\eta }^{\prime 2}-{\chi }^{\prime }}}=-\frac{1}{M_{\eta }^{\prime }}.$$

$$\frac{a}{b}=-\frac{M_{\eta }-E+\eta \lambda (q_x+i{q}_z)}{M_{\eta }-E+\eta \lambda (q_x-i{q}_z)}.$$

$$|\psi ⟩$$

$$N=$$

$$+\mathrm{Im}\left[(K+\eta f{)}^2(K^2+g^2)\left(\frac{e^{-2iL{q}_z}-1}{q_z}\right)\right].$$

$$N=$$

$$+[(K+\eta g{)}^2(f^2+K^2)e^{-\kappa L}$$

$$+(K+\eta f{)}^2(g^2+K^2)e^{\kappa L}]\frac{\sinh (\kappa L)}{\kappa }.$$

$$H_{\text{RKKY}}$$

$$\times G(-{\mathbf{r}}_{12};\omega +i0+)]$$

$$\equiv \sum _{i,j=x,y,z}{J}_{ij}{S}_{1i}{S}_{2j},$$

$$H({\vec{k}}_{||})=\sum _j$$

$$+{\psi }_j^{†}({\vec{k}}_{||})h_{j,j}({\vec{k}}_{||}){\psi }_j({\vec{k}}_{||})),$$

$$H$$

$$\left(\omega \mathbb{I}-H({\vec{k}}_{||})\right)G({\vec{k}}_{||},\omega )=\mathbb{I}.$$

$$\left(\omega \mathbb{I}-h^{\prime }\right)G^{\prime }=\mathbb{I},$$

$$\displaystyle h^{\prime}=\left(\begin{array}[]{cc}h^{(k)}&A^{(k)}\\ A^{\dagger(k)}&h^{(k)}\\ \end{array}\right),\quad G^{\prime}=\left(\begin{array}[]{cc}G_{11}&G_{1N_{z}}\\ G_{N_{z}1}&G_{1N_{z}}\\ \end{array}\right),$$

$$A^{(i+1)}=A^{(i)}(\omega -h^{(0)}{)}^{-1}{A}^{(i)},$$

$$h^{(i+1)}=h^{(i)}+A^{†(i)}(\omega -h^{(i)})A^{(i)},$$

Full text

I Introduction

Weyl semimetals (WSMs) are three dimensional topological systems that host an even number of band-touching points (termed as Weyl nodes) in the bulk spectrum, near which the low-energy excitations follow the relativistic Weyl equation WSMreview . Such Weyl quasiparticles have definite chirality and the chirality of these quasiparticles are given by the nature of the Weyl nodes, which can act as either sources or sinks of Berry curvature in the Brillouin zone. In a finite geometry, WSMs also host unique surface states, named as Fermi arc states, whose projected Fermi surfaces are open arcs on each of the surfaces. Numerous materials have been predicted to be suitable candidates for WSMs , and a variety of experiments demonstrate their novel character.

Correlations functions impact many properties of these systems, and are of special interest because of the unique helical nature of low-energy excitations in WSMs. Within this class of properties, the effective interaction between two localized impurity spins introduced in such a system, mediated by the WSM conduction electrons, is described by the Ruderman-Kittel-Kasuya-Yoshida (RKKY) theory RKKY , and is directly related to the spin-spin correlation function of electrons within the system. Among solid-state materials, spin-orbit coupled systems  RKKYSO1 ; RKKYSO2 , particularly topological ones, are particularly interesting in the way they mediate long-range – and sometimes controllable – couplingRKKYQSH1 ; RKKYQSH2 ; RKKYElectrical1 ; RKKY2DTI ; RKKY2DTI2 ; RKKYTI ; Reja_2017 ; Reja_2019 ; SonuRKKY among spins. Work on RKKY interactions through bulk Weyl fermions RKKYWSM1 ; RKKYWSM2 ; RKKYWSM3 show the interactions can be anisotropic and are in some circumstances weak, but generally carry signatures of the chiral nodes.

Significant attention has also been given to RKKY interactions on the surface of WSMs RKKYWSMSur1 ; RKKYWSMSur2 ; RKKYWSMSur3 . Surface states, at the Fermi energy of a WSM in a slab geometry, typically reside at wavevectors which form arcs in the surface Brillouin zone. These arcs join one Weyl node to the other, and typically disperse energetically perpendicular to a given arc, with different signs of the dispersion for each of the two physical surfaces. The essentially one dimensional character of the surface states results in a strong, highly anisotropic spin-spin correlation function, with similarly anisotropic RKKY interactions among spin impurities adsorbed on the surface of a WSM. Previous studies have been largely confined to semi-infinite geometries, for which coupling between surfaces cannot be modeled. Such coupling is potentially significant, as the Fermi arc states can be relatively weakly localized at their surfaces; moreover, the penetration length of a Fermi arc state diverges as the surface wavevector approaches the projection of a Weyl node. This can induce interesting physics due to non-negligible coupling between spins on opposite surfaces of the WSM.

In this paper, we investigate the effects of such intersurface coupling by analyzing a slab geometry of finite thickness. Specifically, we examine effective spin-spin interactions due to the RKKY mechanism for two different situations. Firstly, when two spins are on the same surface of the WSM, surface electrons on the opposite surface can participate in their coupling. The resulting RKKY interaction reflects the anisotropy of the Fermi surface, and in the thick slab limit recovers previous results in which only a single surface was modeled RKKYWSMSur1 . Secondly, when the spin impurities are on opposite surfaces of the slab, the resulting coupling depends strongly on the overlap of the Fermi arc states. We find that the resulting coupling is a non-monotonic function of the slab thickness, and a thin-film limit can be defined by observing when the coupling between the spins is strongest. In order to compute the coupling, we have developed a recursive Green’s function scheme in which the elements of a Green’s function on the surfaces can be computed essentially exactly with relatively high numerical efficiency. We show that analytical solutions for the WSM wavefunctions in a slab geometry verify the numerical results, and offer some insight into their qualitative behavior.

This manuscript is organized as follows. In Section II, we introduce the simple WSM model used for our work and find analytical solutions for wavefunctions in a slab geometry with appropriate boundary conditions. In Section III, we briefly discuss the formal expression for RKKY interactions and our numerical scheme for computing them in a slab geometry of a tight-binding model. Our numerical results are presented in Section IV, along with a comparison with analytical results. Finally we conclude with a summary and discussion in Section V.

II Weyl Semimetal - Thin Film

II.1 Model Hamiltonian

A minimal model of a WSM has two Weyl nodes at the Fermi energy and breaks time-reversal (TR) symmetry. For such a model, the low-energy Hamiltonian can be written using a two-band model. If the two bands represent spin states, then for a slab geometry, with the Weyl nodes separated along the momentum of one of the translational invariant (in-plane) directions, the surface states (Fermi arcs) are spin-polarized, resulting in completely spin-polarized surfaces of the slab. As the indirect spin-exchange interaction is only interesting when the ground-state is spin-unpolarized, the minimal model we consider must have at least two Fermi arcs on each surface, with the spin-polarizations of each oriented such that the net spin density on either surface vanishes. If the two Fermi arcs on a surface have distinct locations in the surface Brillouin zone, then one has a total of four Weyl nodes in the bulk, each with a distinct location in momentum space. If the Fermi arcs join two Dirac nodes, then the Fermi arcs will overlap in the surface Brillouin zone. This latter situation was considered, for example, in Ref. RKKYWSMSur1 . In our work we confine our studies to the former case (i.e., Weyl semimetals).

Our starting point is a model Hamiltonian defined on a cubic lattice WSMModel . The Hamiltonian preserves time-reversal symmetry (defined by the time reversal operator $T=\textit{i}\sigma_{y}K$, with $K$ the complex conjugation operator and $\sigma_{y}$ a Pauli matrix acting in the spin-space), but breaks inversion symmetry, and so has four degenerate Weyl nodes. Specifically, we take

[TABLE]

Here $M_{k}=m+2-\cos{k_{x}}-\cos{k_{z}}$, and $\tau^{\alpha}$ are Pauli matrices acting in an orbital space. For $|m|\leqslant\lambda$ the four Weyl nodes are located at $\textbf{{k}}=(0,\pm\pi/2\pm k_{0},0)$, where $k_{0}=\pi/2-\sin^{-1}(m/\lambda)$. On a given surface, the two Fermi arcs join the four Weyl nodes in a pairwise fashion, as illustrated in Fig. (1). States of the two Fermi arcs are spin-polarized along $\sigma_{x}$ in opposite directions (i.e, they are eigenvectors of $\sigma_{x}$ with opposite eigenvalues for the two Fermi arcs). Furthermore, the two Fermi arcs, at low-energy, are dispersionless along the $k_{y}$ direction and have opposite velocities along the $x$ direction.

The Hamiltonian can be brought into a block diagonal form. Writing $H^{\prime}=UHU^{\dagger}$, with the unitary matrix $U$ defined by

[TABLE]

one finds $H^{\prime}$ has two $2\times 2$ blocks, where for each block (labeled by $\eta=\pm 1$), the two-band Hamiltonian is

[TABLE]

This is a particularly useful form, in which each block individually breaks time-reversal (TR) symmetry, while $T$ maps $H^{\prime}_{+}$ to $H^{\prime}_{-}$ (and vice-versa), so that the total Hamiltonian is TR symmetric. Each of the blocks has two Weyl nodes separated in momentum space, and on a given surface they are joined by one Fermi arc. In principle, a system hosting many Fermi arcs on a surface should be structured in such a way that each joins two Weyl nodes; an effective model of such a multi-Weyl node system could be written effectively as $H=H_{1}\otimes H_{2}\otimes\cdots$, where each of the blocks contains two Weyl nodes.

To focus on the physics of the Fermi arcs, we expand the Hamiltonian to lowest non-trivial order in $k_{x}$ and $k_{z}$, writing $k_{x}\rightarrow q_{x}$ and $k_{z}\rightarrow q_{z}$. Then from Eq. (II.1) we obtain

[TABLE]

with $M_{\eta}(k_{y})=m\eta-\lambda\sin k_{y}$. The four Weyl nodes are at ${\bf K}_{\eta,\xi}=(0,\eta\frac{\pi}{2}+\xi k_{0},0)$ with $\eta,\xi=\pm 1$ and $k_{0}=\cos^{-1}(m/\lambda)$. For the $\eta=+1$ block, $M_{+}<0$ between $k_{y}\in(\pi/2-k_{0},\pi/2+k_{0})$. For a surface perpendicular to the $z$ direction, along the $k_{y}$ axis these two points are connected by a Fermi arc. For the $\eta=-1$ block, $M_{-}>0$ between $k_{y}\in(-\pi/2-k_{0},-\pi/2+k_{0})$, and again there is a Fermi arc connecting these points on the $k_{y}$ axis for the same surface. This situation is illustrated in Fig. 1. Note that for $H^{\prime}$ (i.e., after the unitary transformation), states on the Fermi arcs are eigenvectors of $\sigma_{y}$ rather than $\sigma_{x}$.

Near the Weyl nodes, if we can write the low-energy Hamiltonian in the form of $H=k_{\mu}A_{\mu\nu}\sigma_{\nu}$, then the chirality of the node is given by sgn(Det[$\mathbf{A}$]). Writing, $\textbf{{k}}=(q_{x},\eta\frac{\pi}{2}+\xi k_{0}+q_{y},q_{z})$ and expanding to first order in $q_{i}$, we arrive at the low-energy Hamiltonian

[TABLE]

with $\alpha=\sqrt{1-(m/\lambda)^{2}}$. The chiralities of the four nodes may then be written as ${\rm sgn}({\rm Det}[\mathbf{A}_{\eta,\sigma}])=-\eta\xi$.

II.2 Infinite mass boundary condition

To make progress analytically, we need to construct appropriate boundary conditions of the Dirac Hamiltonian Eq. (II.1) for a slab geometry, such that the properties of the Fermi arc can be recovered. In general boundary conditions for Dirac equation can be cumbersome [REF], but our goal is to recover the properties of the surface modes (i.e, Fermi arc states). We construct boundary conditions by taking the Hamiltonian of the vacuum (outside the slab, which extends from $z=0$ to $z=L_{z}$), similar to Eq. (8), except for the mass term, whose form is taken as $M^{\rm vac}_{\eta}=\eta m_{0}$, with $m_{0}\rightarrow\infty$. This construction is required to ensure that for momentum between the Weyl nodes the effective mass term ($M_{\eta}(k_{y})$) for the Weyl semimetal and the vacuum ($M_{\eta}^{\rm vac}$) are oppositely signed.

The eigenfunctions for the Hamiltonian $H_{\rm vac}=\lambda(q_{x}\sigma_{y}-q_{z}\sigma_{x})+M_{\eta}^{{\rm vac}}\sigma_{z}$ are

[TABLE]

with eigenvalue $E=\pm\sqrt{m_{0}^{2}+\lambda^{2}(q_{z}^{2}+q_{x}^{2})}$. For $m_{0}\gg E$, the eigenfunctions are normalizeable if

[TABLE]

with $\kappa=\sqrt{m_{0}^{2}+q_{x}^{2}-E^{2}}$. Thus, in the limit $m_{0}\rightarrow\infty$, we have $\kappa\rightarrow m_{0}$. For $z>L_{z}$,

[TABLE]

For $z<0$,

[TABLE]

At $z=0,L_{z}$, these spinors become the Fermi arc wavefunctions, and are recognizable as eigenvectors of $\sigma_{y}$.

Matching the wavefunction $\psi(z)$ within the slab to these boundary forms yields the conditions

[TABLE]

where

[TABLE]

with $q_{z}=(1/\lambda)\sqrt{E^{2}-M_{\eta}^{2}-\lambda^{2}q_{x}^{2}}$. Non-trivial solutions of Eq. 23 exists if

[TABLE]

Simplifying this condition, we obtain a transcendental equation,

[TABLE]

where $\chi=(E/\lambda)^{2}-q_{x}^{2}$. For all real solutions $\chi$ of this equation, the energy has values $E=\pm\lambda\sqrt{\chi+q_{x}^{2}}$. No solutions of Eq. 35 exist with $\chi<0$.

For bound-state solutions, i.e, when $q_{z}$ is imaginary, $\chi<(M_{\eta}/\lambda)^{2}$. The left hand side of Eq. 35 is a positive function with values between 0 and 1. Thus, such bound-state solutions are only possible when $\eta M_{\eta}<0$ as well as when $|L_{z}M_{\eta}/\lambda|>1$, i.e., when $|M_{\eta}|>\lambda/L_{z}$. Defining $\chi^{\prime}=L_{z}^{2}\chi$ and $M^{\prime}_{\eta}=\eta L_{z}M_{\eta}/\lambda$, we rewrite Eq. 35 as

[TABLE]

The various solutions of $\chi^{\prime}$ from Eq. 36 can be labeled by an index $n=0,1,..$ (with increasing values of $n$ corresponding to larger values of $\chi^{\prime}$) and the corresponding energy solutions $E_{n,\pm}(q_{x},q_{y})=\pm\lambda\sqrt{\chi_{n}+q_{x}^{2}}$ gives rise to particle-hole symmetric bands. The minimum solution of $\chi^{\prime}$ is shown in Fig. 2. The bands with $n=0$ contains all the Fermi arc states (when $k_{y}$ is between the Weyl nodes, in the Fermi arc interval) and low-energy bulk states (when $k_{y}$ is outside the interval).

The coefficient $(a/b)$ for the states, Eq. (26), can be found from the boundary conditions at $z=0$ to be

[TABLE]

We can then write down the wavefunctions. Defining $K=M_{\eta}(k_{y})-E,f=\lambda(q_{x}-iq_{z}),g=\lambda(q_{x}+iq_{z})$, one finds

[TABLE]

For real $q_{z}=\sqrt{\chi-M_{\eta}(k_{y})^{2}}$ (when $\chi>m^{2}$, $f=g^{*}$) the normalization factor has the form

[TABLE]

For purely imaginary $q_{z}=i\kappa$ (when $\chi<m^{2}$),$f=q_{x}+\kappa$, $g=q_{x}-\kappa$,

[TABLE]

These are the full solutions of the low-energy states of the WSM slab in the rotated basis (Eq. (6)). Once written in the original basis, these solutions correctly reproduce the spin configuration of the Fermi arc states.

III RKKY interaction and Recursive Green’s function

Ruderman-Kittel-Kasuya-Yoshida (RKKY) theory RKKY describes the effective coupling between two impurity spins ${\bf S}_{1}$ and ${\bf S}_{2}$ in a metal mediated by the conduction electrons. The spins, located respectively at ${\bf r}_{1}$ and ${\bf r}_{2}$, are typically treated as classical magnetic moments, and are assumed to be coupled to the electrons by $sd$ Hamiltonians, $H_{sd}=J{\bf S}_{i}\cdot{\bf s}({\bf r}_{i})$ ($i=1,2$), where ${\bf s}({\bf r}_{i})$ is the conduction electron spin density at the location of impurity spin $i$. For small $J$ the resulting impurity spin interaction becomes

[TABLE]

where $\mathbf{r}_{12}$ is the separation of the two spins and $G(\mathbf{r}_{12};\omega+i0+)$ is the real space Green’s function for the unperturbed electron system. The resultant $J_{ij}$ is essentially the electronic spin-spin correlation matrix. In all of our results we show $J_{ij}$ in the unit of $J^{2}$.

Details of the particular electron system in which the impurity spins are embedded enter the calculation through $G(\mathbf{r}_{12};\omega+i0+)$. For our WSM system, we will proceed in two ways. First, we will directly compute $G$ in momentum space from the low-energy Hamiltonian wavefunctions Eq. (42), and then Fourier transform the expression to obtain the needed real-space Green’s function. Our second approach is more numerical, and involves inverting the tight-binding model, Eq. (1). In this approach the discrete translational invariance in the $x-y$ plane of the slab geometry allows, for each two dimensional wavevector, independent computation of the Green’s function. A computation of the real space Green’s function then follows from a Fourier transform. For the results we present in the next section, we also restrict ourselves to considering impurities which are exchange-coupled to the same orbital of the two-orbital model, Eq. 1, which captures the essential physics of interest. For the case of the semi-analytical model of the last section, for which the atomic-scale structure is not included, we assume the impurities to be exchange-coupled to the conduction electrons within a small region (of thickness of one lattice spacing) on each surface.

Before proceeding to our results, we use the remainder of this section to outline the recursive Green’s function method we use for our fully numerical studies. We are interested in the coupling between impurities placed on the surfaces, so that in the computation of $G(\mathbf{r}_{12};\omega+i0+)$ one only actually needs the Green’s function for sites ${\bf r}_{1}$ and ${\bf r}_{2}$ on the slab surfaces. Following Ref. Keshav, , we can compute the two dimensional Fourier transform of this, $G_{ij}(\omega,k_{x},k_{y})$, where $i$ and $j$ label the surfaces of the slab on which ${\bf r}_{1}$ and ${\bf r}_{2}$ reside, respectively.

For a slab geometry of $N_{z}$ number of sites in the $z$ direction, we re-write the tight-binding Hamiltonian (Eq. 1) in the form

[TABLE]

where $\vec{k}_{||}=(k_{x},k_{y})$, which are good quantum numbers. This allows us to write the Hamiltonian in the form of

[TABLE]

and the Green’s function is evaluated from the equation

[TABLE]

When $N_{z}=1+2^{k}$, the above set of equations can be recast in the form

[TABLE]

with

[TABLE]

where the $h_{t}^{(k)},h_{b}^{(k)}$ and $A^{(k)}$ are found by recursively solving

[TABLE]

with $h_{t}^{(0)}=h_{b}^{(0)}=h^{(0)}=h$ and $A^{(0)}=A$. This yields the two surface Green’s functions $G(\vec{k}_{||})_{1,1}$ and $G(\vec{k}_{||})_{N_{z},N_{z}}$ as well as their connections $G(\vec{k}_{||})_{1,N_{z}}$ and $G(\vec{k}_{||})_{N_{z},1}$ without requiring a solution for the full Green’s function.

IV Numerical results

In this section, we present our results for the specific cases when (i) the magnetic impurities are on the same surface of the WSM, and (ii) when the magnetic impurities are on opposite surfaces of the WSM. In both cases we compute the RKKY interaction using the recursive Green’s function method outlined in the last section, as well as using the analytical wavefunctions of the $n=0$ bands of Sec. II, which contain the Fermi arc states, and compare the results.

IV.1 Impurities on a single surface

As discussed above, the Fermi arc states disperse in energy along $\vec{k}$ perpendicular to the arc itself. As a result these states have a highly asymmetric velocity, with $v_{y}\ll v_{x}$. The effective interaction between spin impurities placed on the same surface reflects this strong asymmetry. For a semi-infinite system, it can be shown for large $r$ that the elements of the Green’s function $G(\vec{r})$ asymptotically vanish as $\sim 1/r^{2}$ when $\vec{r}$ is on the surface and $\vec{r}\parallel\hat{y}$. By contrast, they fall off as $~{}\sim 1/r$ when $\vec{r}$ points along the $\hat{x}$ direction RKKYWSMSur1 ; the difference is a consequence of the (nearly) unidirectional dispersion of the Fermi arc energies. This results in the strongest RKKY coupling for impurities separated along the $\hat{x}$ direction, and in what follows we focus on separations along this direction. Moreover, for a thick enough sample (when the presence of the other surface may be neglected), states in each of the Fermi arcs are spin-polarized (along the direction of $\sigma_{x}$ in our model) and are chiral in their dispersion (i.e, the energy is proportional to $\pm k_{x}$ for $\eta=\pm 1$). One expects in this case that the RKKY interaction $J_{xx}$ will be vanishingly small RKKYQSH1 ; SonuRKKY ; RKKYWSMSur1 .

These expectations may be understood as follows. The Fermi arcs, for a thick sample, are exponentially confined to a surface at $z=0$ with an approximate wave function (see Eqs  (17) and (22)):

[TABLE]

where $M_{\eta}(k_{y})=m\eta-\lambda\sin(k_{y})$. These wavefunctions allow us to write an effective Green’s function from the Fermi arc on this two-dimensional surface in the form

[TABLE]

where the Fermi arcs exists between $k_{\eta,1}=\eta\frac{\pi}{2}-k_{0}$ and $k_{\eta,2}=\eta\frac{\pi}{2}+k_{0}$. The $k_{y}$ momentum can be integrated between $k_{\eta,1}$ and $k_{\eta,2}$ and performing the $k_{x}$ integral one obtains,

[TABLE]

where, the Fermi velocities of the Weyl nodes $v_{F}=\lambda$. In the limit $\delta\rightarrow 0+$,

[TABLE]

and r = $(x,y)$. This approximate form for the Green’s function is useful in determining the asymptotic behavior of the RKKY interaction, as we show briefly in the next section (see also Ref. RKKYWSMSur1, ). Note that the step function in Eq. (69) implements the chiralities of the Fermi arcs. As $\eta=\pm$, in the full 4$\times$4 orbital/spin space, the Green’s function is then

[TABLE]

With this expression it is straightforward to work out the RKKY integral (Eq. 45) and show that, in the original basis of the Hamiltonian, Eq. (1), the elements of the correlation matrix, Eq. (46), are given by

[TABLE]

Other off-diagonal components vanish. In this limit, the separation of the pairs of the Weyl nodes (taken as $\pi/a$ in our model) does not alter the result.

The symmetries apparent in Eqs. 76 can be understood from the original tight-binding system as we discuss in Appendix A, and are captured by our results, as illustrated in Fig. 3. In particular the value of $J_{xx}$ decreases with increasing thickness of the sample whereas the other diagonal components become constant with increasing thickness. In contrast to the analytical result, $J_{zz}$ and $J_{yy}$ fall off asymptotically roughly as $1/x^{2}$, as listed in Table 1. This difference from Eqs. 76 is expected due to the curvature of the Fermi arcs RKKYWSMSur1 . Off-diagonal components other than $J_{zy}$ and $J_{yz}$ (not shown) are several orders of magnitude smaller than these quantities, as expected from the above analysis. Note also the qualitative agreement between results from our semi-analytical model and the tight-binding computation.

The results for $J_{xx}$ obtained from the tight-binding computation are particularly interesting. In contrast to the result for straight Fermi arcs in a single surface system derived above, in the slab geometry $J_{xx}$ remains non-zero and falls off rather slowly (see Table 1) with the distance between the impurity spins. (Similar behavior is found in our analytical model at small thickness.) The discrepancy can be attributed to two possible effects: (i) interactions mediated by the bulk states which were not included in the simple Fermi arc analysis, and (ii) the presence of the second surface. Interestingly, Fig. 3(b) shows that $J_{xx}$ vanishes rapidly with increasing thickness, which clearly favors mechanism (ii). Fig. 3(d) also manifests a critical thickness $L_{c}$, defined by the width $L_{z}$ at which $J_{xx}$ attains its maximum value before decaying sharply with further increase. This critical thickness can be used to define a “thin-film limit” of the system, for which the effects of having two surfaces are maximal. Noting that $k_{0}$ is the only relevant momentum scale, we expect the thin film limit to scale as $L_{c}\sim 1/k_{0}$. A numerical verification of this hypothesis is presented in Appendix B.

IV.2 Impurities on opposite surfaces

When the impurity spins are put on opposite surfaces of the slab, they may communicate via electron states that are present in the bulk of the Weyl semimetal. To examine this effect numerically, we place the two spins on different surfaces of a WSM slab with various thicknesses and vary their separation along the $x$ axis (i.e., the direction in which the Fermi arcs states disperse). Results from these are illustrated in Figs. 4 and  5.

For the range of parameters we examined, the symmetry properties of the spin coupling matrix turn out to be the same as when the spins are situated in the bulk and are separated along the $z$ direction (see Appendix A). Numerically, when the two spins are located at sites $(x,y,z)=(0,0,0)$ and $(R,0,L_{z}=N_{z}a)$, the coupling between them is surprisingly strong despite the fact that they reside on different surfaces, decaying as $\sim(R^{2}+L_{z}^{2})^{-\alpha/2}$, where $\alpha$ is between 3 and 4, for fixed $L_{z}$ and increasing $R$, showing a rather slow decay of the RKKY coupling and possibility of mean-field magnetic ordering. Fitting the decay with $\sim R^{-\beta}$, where the thickness $L_{z}$ is fixed, the coupling decay in even slower manner, with $\beta$ between 2 and 3. Results from the tight-binding simulation and from the low-energy wavefunctions both support these results.

As a function of the thickness $L_{z}$, for fixed $R$ the couplings initially increase and after attaining maximum values decrease rapidly. Results for both the analytical and tight-binding approaches for varying $N_{z}$ are illustrated in Fig. 5. For fixed $R$, as a function of the thickness the coupling decays as $\sim L_{z}^{-\gamma}$ with $\gamma\approx 5$. This is the same falloff as for RKKY coupling in the bulk of the WSM RKKYWSMSur1 , suggesting that for large enough $L_{z}$ the coupling between spins on opposite surfaces is eventually dominated by the bulk states. The clear differences among the parameters $\alpha,\beta,\gamma$ capture the essential physics of the WSM system with Fermi arcs, and have been listed in the Table II.

The maximum coupling as a function of the thickness can be qualitatively understood as follows. As the thickness increases, the Fermi arcs localize increasingly firmly on the surfaces, increasing the surface density of states near the Fermi energy. This leads to an increase in the coupling between the impurity spins on the surfaces and the conduction electrons, which can mediate intersurface interactions effectively when $L_{z}$ is not too large. On the other hand, as $L_{z}$ increases, the number of conduction electron states which are sensitive to both surfaces decreases, resulting in weaker coupling between spins on opposite surfaces. With increasing thickness, the competition between these two mechanisms gives rise to a critical thickness for which the coupling between spins placed on opposite surfaces maximizes. As in our earlier argument for impurities on the same surface, with $k_{0}$ the only relevant momentum scale we expect this thickness to scale as $\sim 1/k_{0}$. We explore this in Appendix B. Again this critical thickness also defines a thin-film limit; the values of this critical thickness obtained from inter-surface coupling are of the same scale as those obtained from the intra-surface coupling.

V SUMMARY AND DISCUSSION

In this work we have examined RKKY interactions among impurity spins on the surfaces of Weyl semimetal (WSM) slabs, using both an approach in which the wavefunctions of the WSM electrons are found in an analytical form, and a more fully numerical recursive Green’s function technique. We find that Fermi arc surface states play an important role in the RKKY coupling, creating couplings that are stronger and more long-range than is found for impurities well-inside the bulk of the system. Surprisingly, even the coupling between spins on opposite surfaces can be relatively strong. As a function of film thickness, we find that the RKKY couplings are non-monotonic, with maxima that can define a “thin-film” limit, in which the effects of both surfaces are in some sense maximal. The relative strengths and signs of different components of the RKKY couplings $J_{ij}$ can be understood using a simple model in which only surface states associated with the Fermi arcs are retained, and in which the Fermi arcs are perfectly straight.

The importance of electron states with strong support on the surfaces can be examined by comparing results for geometries with Fermi arcs to ones without them. Fig.  6 illustrates RKKY coupling for spins on the same and opposite surfaces which are perpendicular to the direction of separation between the Weyl nodes in the bulk, for which surface states are not present. The generally smaller scale of the resulting couplings supports the idea that the Fermi arc states play a large quantitative role in setting the coupling scale.

The results presented to this point have been for vanishing chemical potential $\mu$, where the only extended Fermi surfaces are due to the Fermi arcs, and the Fermi energy passes directly through the Weyl nodes in the bulk. In general, when $\mu\neq 0$ and the Fermi wavevector $k_{F}\neq 0$ in the bulk, one expects $2k_{F}$ oscillations in the RKKY coupling. Results for $\mu\neq 0$ are presented in Fig.  7, for which the oscillations are apparent. The envelopes within which these oscillations occur behave rather similarly to the results for $\mu=0$.

When a system is of order or thinner than a critical thickness $\sim 1/k_{0}$, our results show that a proper treatment of RKKY interactions requires one to retain states from the Fermi arcs of both surfaces, even if the two spins reside on the same surface. For real systems, such as TaAs taas , the typical separation of Weyl nodes is rather small (of the order of $k_{0}\approx 0.1\pi/a$) and thus we expect the critical thickness to be of the order of several tens to a hundred lattice spacings. Such thicknesses are quite reasonable for thin-film semiconductor systems.

We conclude with some speculations about the kind of magnetic order these RKKY interactions might induce in the low temperature state of spin impurities on the surfaces of a WSM thin film. At large distances, the strongest couplings we find are for $J_{yy}=J_{zz}<0$ within a single surface, suggesting the system will form a planar ferromagnet in its ground state. The non-vanishing $J_{yz}$ and $J_{zy}$ couplings if large enough could induce spiral order; while at short distances these can be larger than the diagonal elements, at long distances the latter are significantly larger. Given the relatively slow spatial decay of the RKKY interaction, it seems likely that the system will favor ferromagnetism. Furthermore, the sign of coupling for impurities on different surfaces suggests that the magnetization of the two surfaces will be parallel to one another in the groundstate. In principle at low temperature such magnetic order should be detectable. Moreover, with this type of order one expects a magnetic disordering transition at finite temperature in the Kosterlitz-Thouless universality class, which might be detected in thermal measurements or via spin transport in the system. Finally, the importance of the Fermi arc states in supporting such magnetic order could be tested by comparing the behavior of slabs in which the surfaces support them to ones in which they do not. We leave the investigation of these questions to future work.

Acknowledgements – HAF thanks the NSF for support through Grant Nos. DMR-1914451, DMR-1506263 and DMR-1506460, by the US-Israel Binational Science Foundation. AK thanks support from SERB (Gov. of India) through grant ECR/2018/001443 and BRNS (Gov. of India) support through grant 58/20/15/2019-BRNS.

Appendix

A. Symmetries among the susceptibility matrix elements

In this appendix we briefly discuss the symmetries among the various coupling elements $J_{ij}$ for spins that are on the same or opposite surfaces, as well as in the bulk, based on the symmetry of the underlying Hamiltonian of the WSM.

First we consider the case when both spins, coupled to same orbital, are on the same surface of a WSM slab of thickness $N_{z}$ with spatial separation $(x=R,0,0)$. The two Green’s functions $G(\textbf{r}_{12},\omega)$ and $G(-\textbf{r}_{12},\omega)$ required to calculate the RKKY coupling between the spins for some arbitrary values of $x$ and $\omega$ have the structures

[TABLE]

where $s_{i}$ (with $i=0,1,2,3,4$) are complex numbers depending on $R$ and $\omega$. Using these in Eq. 46, we obtain all possible nonzero components of the spin-spin correlation matrix to be

[TABLE]

Next, we consider the case when both spins, coupled to the same orbital, are on the opposite surfaces of WSM slab of thickness $L_{Z}$ with lateral spatial separation $R$, the positions of the two spins are (0,0,0) and $(x=R,0,L_{z})$. The two Green’s functions $G(\textbf{r}_{12},\omega)$ and $G(-\textbf{r}_{12},\omega)$ required to calculate the RKKY coupling between the spins for arbitrary values of $R$ and $\omega$ now have the structure

[TABLE]

with $o_{i}$ (with $i=0,1,2,3,4$) complex numbers depending on $R$, $L_{z}$, and $\omega$ ($o_{2}$ is two orders magnitude smaller than other elements). Using these in Eq. 46, we find all possible nonzero components of spin-spin correlation matrix are related by

[TABLE]

Finally, we consider the case when the two spins, coupled to the same orbital, are deep in the bulk of a WSM and have spatial separation $(0,0,z=R)$. The two Green’s functions $G(\textbf{r}_{12},\omega)$ and $G(-\textbf{r}_{12},\omega)$ required to calculate the RKKY coupling between the spins for some arbitrary values of $R$ and $\omega$ have the structure

[TABLE]

where $b_{i}$, (with $i=0,1,2,3,4$) are complex numbers depending on $R$ and $\omega$. Using these in Eq. 46, we obtain the components of the spin-spin correlation matrix as similar as those in Eqs. (A. Symmetries among the susceptibility matrix elements):

[TABLE]

B. Critical thickness and the separation of the Weyl nodes

Finally, we examine in more detail the $k_{0}$ dependence of the critical thicknesses at which the couplings are maximized, as discussed in the main text. Fig. 8 illustrates numerical results for the thickness dependence of RKKY couplings for various values of the Weyl node separation $k_{0}$, for both the cases when the spins are on the same surface as well as when the spins are on the opposite surfaces. If one estimates the critical thickness $L_{c}$ at which the coupling attains its maximum, one finds that $L_{c}k_{0}\approx$ constant.