Tight-binding model with sublattice-asymmetric spin-orbit coupling
for square-net nodal line Dirac semimetals

Topological semimetals (TSMs) are a class of materials whose crossings between valence and conduction bands are topologically protected [1]. There are two main subcategories of TSMs: nodal point semimetals (NPSs), and nodal line semimetals (NLSs). NPSs have a discrete set of crossing points at $E_{F}$, while NLSs display a continuous crossing line. Initial, theory-lead searches for NLSs called for materials with nonsymmorphic symmetry [2, 3] such as ZrSiS [4], ZrSiSe, and ZrSiTe [5]. Those materials display unusually-high magnetoresistance [6, 7, 8, 9, 10], high carrier mobilities [11, 7, 12], and linear band dispersions around $E_{F}$ [4, 6, 13].

ZrSiX (X=S, Se, or Te) belongs to the nonsymmorphic P4/nmm space group. As shown in Fig. 1, two atoms of each species are present in their u.c., with their Si atoms arranged in stacked square nets [14, 15, 16, 17, 18, 19, 20]. The electronic structure of two-site square lattices is described by $\ket{p_{x}}$ and $\ket{p_{y}}$ orbitals in various models. Indeed, such a basis set was proposed by Luo and Xiang to study a Bi-based topological insulator [21]. Later, Klemenz and coworkers discovered nodal line semimetals in square nets [22, 23, 24], which are 2D symmorphic lattices, and presented a symmetry analysis of such nets relying on $\ket{p_{x}}$ and $\ket{p_{y}}$ orbitals in Ref. [24] (Fig. 2(a)). Since then, other two additional teams have implemented two-site $p_{x}p_{y}$ models to study the effect of SOC in similar nets [25, 26].

Recent quantum oscillation experiments have shown that the Berry phase generated for ZrSiS depends on the orientation of the magnetic field [8, 9] and Yang et al. suggested that a key agent in such dependence is SOC [27]. It is natural to include SOC in $p_{x}p_{y}$ models to understand the observations made in Ref. [27]. It is not possible, nevertheless, to gap the crossings at the Fermi level in these models upon the customary implementation of SOC [28]. Luo and Xiang [21], Deng et al. [25], and Aryal et al. [26] independently realized that $p_{x}p_{y}$ square net models can be gapped at zero energy through a combination of (i) an on-site SOC (this is, one in which neighboring atoms do not contribute to SOC), and (ii) a lowering of symmetry such that the two atoms in the u.c. turn inequivalent. Such asymmetry was created by a relative displacement of one atom along the 2D plane (Fig 2(b)) [25], or by an out-of-plane relative displacement (buckling, Fig 2(c)) [21]. On the other hand, Aryal et al. incorporated an on-site asymmetry between sublattices directly onto the SOC interaction ad hoc (meaning that no microscopic justification was offered to motivate their choice for SOC interaction).

Given that the Si square nets on ZrSiS, ZrSiSe, and ZrSiTe are in fact not distorted in the way done in Refs. [21] and [25], and that Aryal’s work did not provide justification for their SOC, we seek the source of an asymmetric SOC that could be added to the $p_{x}p_{y}$ models to produce a band gap at zero energy, and thus make it consistent with ab initio results containing SOC. The search for that mechanism from ab initio information is our guiding principle, and the motivation for this work.

We have previously worked with 2D models for topological insulators relying on $\ket{p_{x}}$ and $\ket{p_{y}}$ orbitals also [29], and we possess a deep knowledge of ab initio-informed TB formulations [30] as provided in the SIESTA code [31, 32]. We have written extensions to this code for quantum transport [33] and ab initio molecular dynamics [34] calculations. We are also closely familiar with SOC on this code, down to implementing it [35, 36, 37, 38]. Our knowledge of TB models [39, 40, 41] makes the extraction of TB parameters [42] from SIESTA possible, to the point of turning on-site SOC interactions for analysis at will. This experience underpins the description of a SOC-induced mechanism to break sublattice symmetry.

The article is organized as follows: Sec. II contains a brief overview of the $p_{x}p_{y}$ models, and three previous implementations of SOC on them. We demonstrate that the leading contribution of SOC on ZrSiS is due to on-site interactions in Sec. III. To simplify the model, an on-site SOC will be used from that point onwards. A DFT-based TB model is developed in Sec IV by a gradual removal of numerical atomic orbitals from the basis set; we determine the orbitals that contribute to the electronic dispersion of the material around $E_{F}$ the most, while preserving the SOC band gaps. We also determine Zr to contribute to the on-site SOC the most. Applying the Löwdin partitioning scheme to this auxiliary TB model, we effectively induce SOC-gaps into the $p_{x}p_{y}$ model provided in Ref. [24] that break sublattice symmetry in Sec. V. Conclusions are presented in Sec. VI.

Although SOC is not constrained to be an on-site effect [36, 44], Cuadrado et al. have shown that this approximation is usually sufficient [37]. We shall carry a detailed study of SOC in ZrSiS relying on the DFT SIESTA code [31]. The code produces TB representations of the electronic structure with Hamiltonian matrix elements calculated at the ab initio level onto a basis of numeric atomic-like orbitals [32]. This TB representation is generated without the need for post-processing (read Wannierization).

A so-called double-$\zeta$-plus polarization basis set (DZP, or standard) is used in SIESTA calculations [32]. It contains two radial functions per angular momentum in the valence, and one extra (“polarizing”) channel for an unoccupied (excited) orbital momentum channel. For Si and S, the DZP basis set includes two orbitals for the $s-$channel, six for the $p-$channel, and five for the (polarizing) $d-$channel, leading to thirteen atomic-like orbitals for each of those two species [32]. For the Zr atom, the basis contains two radial functions for the $s-$ and $d-$channels, and one radial function for the $p-$channel, leading to 15 orbitals per atom. This way, ZrSiS is described with a basis set containing 82 localized orbitals.

SIESTA uses relativistic pseudopotentials [45] and SOC can be introduced in two gradual ways: either (i) through the on-site approximation, which only couples orbitals within the same atom [36], or (ii) by including matrix elements that couple orbitals belonging to the same atom but also to other atoms, i.e. off-site terms. We dubbed this second approach a full SOC.

The electronic dispersion of ZrSiS with and without SOC corrections and a standard (DZP) basis set is depicted in Fig. 4. The first Brillouin zone is displayed in Fig. 4(a), Fig. 4(b) contains an electronic dispersion without SOC, while Fig. 4(c) and Fig. 4(d) are band structures with on-site and full SOC, respectively. Visual comparison of Figs. 4(c) and 4(d) indicates that the predominant contribution of SOC comes from on-site interactions (in particular, the gaps at the $\Gamma-X$ segment for the on-site and full SOC are both approximately $44$ meV).

Fig. 4 shows that contributions beyond the on-site interaction are small already, but we shall be more rigorous and obtain a band structure associated to a system only perturbed by the off-site SOC next. Even though SIESTA does not explicitly provide this option, it stores Hamiltonian and overlap matrices which can thus be examined and manipulated at will. We obtain the off-site SOC dispersion by defining its associated Hamiltonian and overlap matrices as:

where $H_{0}$ ($S_{0}$), $H^{onsite}$ ($S^{onsite}$), $H^{offsite}$ ($S^{offsite}$) and $H^{full}$ ($S^{full}$) are the Hamiltonian (overlap) matrices without SOC, with on-site SOC, under off-site, and under full SOC, respectively.

Fig. 5(a) depicts the band structure when only the off-site SOC is turned on. Fig. 5(b) is a zoom-in around the vicinity of the crossing along the $\Gamma-X$ segment, where the bands with on-site and full SOC have been included for comparison. The energy gap with just off-site SOC is $\sim 0.3$ meV, a value two orders of magnitude smaller than the gap of roughly 20 meV reported in Refs. [4, 46] and the one seen in Fig. 4. Consequently, the gap is due to the on-site SOC predominantly.

The fact that the band gaps around $E_{F}$ in ZrSiS NLS are predominantly due to an on-site SOC invites to extend the basis of the $p_{x}p_{y}$ model beyond Si $\ket{p_{x}}$ and $\ket{p_{y}}$ orbitals, which will permit observing SOC-gaps without lowering the symmetry of the Si square net [21, 25], and without the need to postulate a SOC interaction that breaks sublattice symmetry ad-hoc [26]. We will develop a TB Hamiltonian with less than 82 orbitals (but more than the four ones used in the $p_{x}p_{y}$ models) in Sec. IV. That TB model will reproduce SOC-induced energy gaps due to an on-site SOC. We will embed this SOC into the $p_{x}p_{y}$ model through Löwdin partitioning method subsequently (Sec. V).

Löwdin partitioning will now be used [49, 50] to project the on-site SOC due to Zr within the 16-orbital model onto the (4-orbital) Si $p_{x}p_{y}$ Hamiltonian at the vicinities of their crossing points ($\mathbf{k}_{16}^{cross}$ and $\mathbf{k}_{p_{x}p_{y}}^{cross}$, respectively). The first steps are a diagonalization without SOC of both the 16-orbital Hamiltonian ($H_{16}\otimes\sigma_{0}$) and the $p_{x}p_{y}$ one ($H_{pxpy}^{4\times 4}\otimes\sigma_{0}$), and a rearrangement of the diagonalized Hamiltonians, such that the two double degenerate bands around $E_{F}$ (that we call “low-energy” bands) appear at the four uppermost entries:

(the curly $\mathcal{H}$ emphasizes that those Hamiltonians without SOC are all diagonal, and $\Delta\mathbf{k}$ is a small excursion away from the $k-$points where the crossing takes place). The point is for the “low-energy blocks”, $\mathcal{H}_{p_{x}p_{y}}^{low}$ and $\mathcal{H}_{16}^{low}$, to describe the same two double degenerate bands despite them belonging to different models:

up to a scaling factor.

To be specific, let $U$ and $V$ be the unitary matrices that diagonalize the unperturbed $p_{x}p_{y}$ and 16-orbital Hamiltonians, respectively, such that the eigenvalues are ordered at $E_{F}$, with eigenvectors $\ket{\psi_{1},\pm}$ and $\ket{\psi_{2},\pm}$ (plus and minus signs label the two orthogonal spin components) (Fig. 9(a)) taking the uppermost four entries:

To add SOC to $\mathcal{H}_{16}$ we must express it in the same basis:

where the superscript in $H^{Zr}_{SOC}$ emphasizes that it is only Zr on-site SOC what is being taken into account. As the low-energy bands are gapped under Zr on-site SOC in the 16-orbital model, its low-energy block will be perturbed, however, $\mathcal{P}$ in general will not be block diagonal, which implies that $\mathcal{H}=\mathcal{H}_{16}+\mathcal{P}$ will not be block diagonal either.

Löwdin partitioning method transforms arbitrary perturbations into a block-diagonal form, nevertheless, so that the vector spaces associated to the low and high energy bands are uncoupled again. Furthermore, the zeroth-order Löwdin partition scheme is to neglect non-block-diagonal terms, giving rise to the gapped electronic dispersion of the low-energy block shown in Fig. 9(b), which is a promising step towards a 4-orbital $p_{x}p_{y}$ model with SOC. Our numerical results lead to a SOC perturbation of the low-energy matrix block of the 16-orbital Hamiltonian with SOC [41]:

at the vicinity of the band crossing, and in the basis of eigenvectors $\{\ket{\psi_{1},+},\ket{\psi_{1},-},\ket{\psi_{2},+},\ket{\psi_{2},-}\}$, which are explicitly provided in Table 3. Here, $\eta$ is a complex parameter that depends on $\mathbf{k}$; $\eta(\mathbf{k}_{16}^{cross})=(2+30i)$ meV. We add this perturbation onto the low-energy block $\mathcal{H}_{p_{x}p_{y}}^{low}$ of the $p_{x}p_{y}$ model in Fig. 9(c), thus finally embedding it with a DFT-based, effective SOC coupling.

Now, given that the SOC perturbation only affected the bands that produced the crossing, the remaining two bands in the $p_{x}p_{y}$ dispersion are still unmodified. In Fig. 10(a), we show the bands of the $p_{x}p_{y}$ model once $\mathcal{P}^{low}$ has been added to its low-energy block. Löwdin procedure only addresses how the dispersion around $E_{F}$ is modified, however, a characteristic of ZrSiS is the symmetry-protected crossing at $X$, which is not observed in Fig. 10(a). To recover such degeneracy, the two high-energy bands in the $p_{x}p_{y}$ model (Fig. 3) must also be modified. Analogously to the low-energy block, we add an high-energy SOC Hamiltonian:

which is now written in the basis $\{\ket{\psi_{3},+}$, $\ket{\psi_{3},-}$, $\ket{\psi_{4},+}$, $\ket{\psi_{4},-}\}$. Thus, the effective on-site SOC inherited from Zr in the $p_{x}p_{y}$ model can be expressed as:

and the Hamiltonian under the effective on-site SOC expressed in the eigenbasis of the unperturbed $p_{x}p_{y}$ Hamiltonian happens to be:

with an associated dispersion degenerate at $X$ depicted in Fig 10(b) (see a similar degeneracy in Figs. 4(c), 4(d), 6(b), and 6(d)).

The last question is whether SOC makes the Si $A$ and $B$ sublattices inequivalent. To answer this, we transformed the $p_{x}p_{y}$ SOC Hamiltonian $\mathcal{P}_{p_{x}p_{y}}$ back into the original basis ($\{\ket{p_{x}^{A},+},\ket{p_{y}^{A},+},\ket{p_{x}^{B},+},\ket{p_{y}^{B},+}$, $\ket{p_{x}^{A},-},\ket{p_{y}^{A},-},\ket{p_{x}^{B},-},\ket{p_{y}^{B},-}\}$), to get:

where $Q$ has the form:

along the $\Gamma-X$ segment: Zr’s on-site SOC did produce an inequivalence among A and B sublattices on the $p_{x}p_{y}$ model, as the hopping from $\ket{p_{x}^{A},\pm}$ to $\ket{p_{y}^{A},\pm}$ has opposite sign with respect to the hopping from $\ket{p_{x}^{B},\pm}$ to $\ket{p_{y}^{B},\pm}$ at the crossing at $\Gamma-X$. Thus, Zr is responsible for the lowering of symmetry that Refs. [21] and [25] assigned to the atoms within their square nets, and which Ref. [26] postulated ad hoc. Fig. 11 displays the band structures of the $p_{x}p_{y}$ models with SOC plus the sublattice asymmetric perturbations discussed in Sec. II for comparison with our model.

The correct, actual SOC Hamiltonian for ZrSiS–obtained with DFT guidance here–can also be expressed in terms of the Pauli matrices introduced in Sec. II, and it is our main contribution:

The picture that arises as a result is as follows: Fig 12 displays the trajectories that an electron originally with spin up follows to connect with atoms on the same sublattice. The electron will first hop to a Zr atom, which will be in charge of switching its spin through on-site SOC. Then, it will hope again to a Si atom of the same sublattice to produce an effective hopping term that acquires the value of $-\operatorname{Re}(\eta)$ in the vicinity of $\mathbf{k}_{cross}$. Observe that, depending on the direction of hopping (horizontal versus vertical), the electron will pass either under or over the Si layer. Since the relative positions of neighboring Zr atoms depends on the sublattice, the presence of Zr facilitates the asymmetry observed in Eq. (21).

Informed by DFT, we have systematically showed how to embed Zr on-site SOC into the $p_{x}p_{y}$ model for ZrSiS. In this manner, we not only connect the observations made in previous works, but also, we show that on-site SOC does not require the addition of extra terms to gap the dispersion around the Fermi level.

In this work, we followed a systematic procedure to induce a gap in the $p_{x}p_{y}$ model. Relying on DFT all along, we concluded that (i) on-site SOC is the main responsible to the SOC-gaps reported on the band structure of ZrSiS around $E_{F}$. Furthermore, (ii) we determined that Zr’s orbitals play a key role on the on-site SOC, and that coupling among opposite spins is crucial to split degenerate bands around $E_{F}$. Then, (iii) we developed a DFT-based TB model larger than $p_{x}p_{y}$, but smaller than the DZP one, which allowed to treat SOC carefully. At that point, (iv) we applied a Lödwin decomposition to extract an effective SOC, and ended up with a $p_{x}p_{y}$ dispersion gapped at $E_{F}$ and consistent with DFT, without arbitrary lowering the symmetry of the square net, nor postulating the entries of the SOC interaction in an ad hoc manner. This model extends the usefulness of the minimal square-net nodal line semimetal $p_{x}p_{y}$ TB models.