Layertronic control of topological states in multilayer metal-organic frameworks

TL;DR

This paper explores how external electric fields can control the localization of topological states in multilayer metal-organic frameworks, revealing new ways to tune electronic properties for materials design.

Contribution

It introduces a method to control layer localization of topological states in multilayer MOFs using electric fields, combining DFT and tight-binding models.

Findings

Electric fields can tune electronic state localization in multilayer MOFs.

Discovery of a V-gutter Dirac state in infinite stacked nanoribbons.

Layer localization adds a new degree of freedom in materials design.

Abstract

We investigate the layer localization control of two-dimensional states in multilayer metal-organic frameworks (MOFs). For finite stackings of (NiC4S4)3 MOFs, the weak van der Waals coupling between adjacent layers leads to a Fermi level dependent distribution of the electronic states in the monolayers. Such distribution is reflected in the topological edge states of multilayer nanoribbons. Moreover, by applying an external electric field, parallel to the stacking direction, the spacial localization of the electronic states can be controlled for a chosen Fermi energy. This localization behavior is studied comparing density functional theory calculations with a kagome lattice tight-binding model. Furthermore, for infinite stacked nanoribbons, a new V-gutter Dirac state is found in the side surfaces, which allows anisotropic current control by tuning the Fermi energy. Our results can be immediately extended to other kagome MOFs with eclipsed stackings, introducing a new degree of freedom (layer localization) to materials design.

Tables

Table 1. Table 1: Lattice parameter a 𝑎 a (Å), averaged inter-layer distance d 𝑑 d (Å), corrugation δ ​ z = ⟨ z 2 ⟩ − ⟨ z ⟩ 2 𝛿 𝑧 delimited-⟨⟩ superscript 𝑧 2 superscript delimited-⟨⟩ 𝑧 2 {\delta z}=\sqrt{\langle z^{2}\rangle-\langle z\rangle^{2}} (Å) and binding energy ( E b superscript 𝐸 𝑏 E^{b} in eV/layer) of the bulk and layered systems.

PBE
System	$a$	$d$	$\delta z$	$E^b$
1ML	14.62	-	0.00	0.000
2ML-GH	14.63	4.23	0.01	0.037
2ML-AB	14.62	3.99	0.05	0.033
2ML-AA	14.62	4.69	0.01	0.022
3ML-AAA	14.64	3.95	0.00	0.004
bulk-GH	14.63	3.79	0.00	0.055
bulk-AB	14.63	3.73	0.00	0.002
bulk-AA	14.63	3.93	0.00	0.009
LDA
System	$a$	$d$	$\delta z$	$E^b$
1ML	14.37	-	0.00	0.000
2ML-GH	14.42	2.96	0.11	0.468
2ML-AB	14.36	3.16	0.06	0.429
2ML-AA	14.44	2.96	0.08	1.051
3ML-AAA	14.43	3.11	0.04	1.138
bulk-GH	14.39	3.13	0.00	0.985
bulk-AB	14.40	3.29	0.00	0.888
bulk-AA	14.45	3.12	0.00	1.800
vdW-D3
System	$a$	$d$	$\delta z$	$E^b$
1ML	14.62	-	0.00	0.000
2ML-GH	14.61	3.26	0.11	0.679
2ML-AB	14.60	3.34	0.15	0.661
2ML-AA	14.64	3.54	0.01	0.860
3ML-AAA	14.64	3.57	0.00	1.206
bulk-GH	14.62	3.31	0.00	1.501
bulk-AB	14.60	3.46	0.00	1.404
bulk-AA	14.65	3.35	0.00	1.968
vdW-DF2
System	$a$	$d$	$\delta z$	$E^b$
1ML	14.83	-	0.00	0.000
2ML-GH	14.83	3.33	0.11	0.645
2ML-AB	14.82	3.34	0.21	0.602
2ML-AA	14.86	3.64	0.01	0.803
3ML-AAA	14.85	3.65	0.00	1.107
bulk-GH	14.84	3.46	0.00	1.360
bulk-AB	14.82	3.56	0.00	1.267
bulk-AA	14.87	3.56	0.00	1.744

Table 2. Table 2: TB parameters used to fit the DFT band structure of Fig.  9 .

System	${ϵ}_0$	$t$	$\alpha$=${\alpha }_z$	$\lambda$	$\Delta$
1L	0.50	-0.08	3.0	0.0	-
2L	0.33	-0.08	3.0	0.0	-0.29
3L	0.42	-0.08	3.0	0.0	-0.29
bulk	0.44	-0.08	3.0	0.0	-0.32

Equations

$$H_N=h_3(\mathbf{k})\otimes {\mathbb{I}}_N+\Delta {\mathbb{I}}_3\otimes \sum _{j=1}^{N-1}|j+1⟩⟨j|+c.c.,$$

$${\lambda }_{\nu }^n(\mathbf{k})=E^n(\mathbf{k})-2|\Delta |\cos \left(\frac{\nu \pi }{N+1}\right),$$

$${\varphi }_{j,\nu }=\sin \left(\frac{j\nu \pi }{N+1}\right)\equiv \sin \left(k_{\nu }{z}_j\right)$$

$$E_{XM}(k_x,k_z)=\hslash v_F{\sigma }_z{k}_x-2|\Delta |{\sigma }_0\cos (k_{z}d),$$

$$H=H_0+H_{SO}+H_L,$$

$$H_0=\sum _L\left[\sum _{i\in L}{ϵ}_0{c}_i^{†}{c}_i+\sum _{(i\ne j)\in L}{t}_{ij}{c}_i^{†}{c}_j\right],$$

$$H_{SO}=\sum _L\left[i\sum _{(i\ne j)\in L}{\lambda }_{ij}{c}_i^{†}\mathbf{\sigma }\cdot \frac{\left({\mathbf{d}}_{kj}\times {\mathbf{d}}_{ik}\right)}{|d_{kj}||d_{ik}|}{c}_j\right].$$

$$H_L=\sum _{L\ne L^{\prime }}\left[\sum _{k\in L}\sum _{l\in L^{\prime }}{t}_{kl}^z{c}_k^{†}{c}_l\right],$$

$$t_{ij}=Nt\exp (-\alpha d_{ij}),$$

$$t_{ij}^z=N_z\Delta \exp (-{\alpha }_z{d}_{ij}),$$

$${\lambda }_{ij}=N\lambda \exp (-\alpha d_{ij}).$$

Full text

I Introduction

Two-dimensional topological insulators (2DTI) have attracted broad interest for novel applications due to their characteristic dissipationless transport properties [1; 2]. While graphene was the first predicted 2DTI, its topological gap is too small to be observed experimentally [3; 4]. Recently, a new class of large gap 2D topological materials, the metal-organic frameworks (MOFs), emerged as an ideal playground to achieve different properties. That is mainly due to a rich tunability of MOFs via different combinations of metals and organic ligands to form designed materials with selected properties [5; 6; 7; 8; 9; 10]. For instance, it has been shown that Mn, Fe and Co-bis(dithiolene) are ferromagnetic half-metals, while n-doped Ru-bis(dithiolene) has a large gap Quantum Anomalous Hall phase [11; 12]. On the other hand, by keeping iron as the metal and changing the organic ligand to a coronene based molecule, a ferromagnetic semicondutor phase has been found [13]. Currently, a great number of MOFs are predicted to host topological phases [14; 15; 16; 17; 18; 19; 20; 21; 22]. Within these materials, nickel and palladium-bis(dithiolene) monolayers were experimentally synthesized [23; 24], and theoretically shown to host a 2DTI phase for n-doped systems [25]. Although the monolayer system need to be n-doped to achieve the 2DTI phase, in a recent work we have shown that for bilayers of nickel and platinum-bis(dithiolene), the topological gap falls onto the Fermi level [26].

The development of novel nano-devices based on 2DTIs is linked to the experimental control of its properties, and manipulation of its charge, spin and valley degrees of freedom [27; 28]. Recently, it has been proposed a bending strain perturbation approach to control spin currents in 2DTIs [29]. Meanwhile, ensuing research has shown layer localization control of topological states by an applied external electric field in MOF bilayer nanoribbons [26] and bilayer graphene grain boundaries [30]. Furthermore, the exploration of van de Waals heterostructures allows for designed properties by combining the stacking of two-dimensional materials [31; 32]. Within the heterostructures, the number of stacked layers and the stacking order has great impact in the electronic properties [33; 34; 35]. Therefore, in order to control the electronic properties of MOFs multilayer, a clear understanding of its stacking is required.

In this paper, we study the stability and control of the electronic properties of the MOF (NiC4S${}_{4})_{3}$ stacking, from a monolayer to the bulk limit. Our first-principles results show that the mirror symmetric AA stacking is the most stable configuration. For a finite AA stacking, each layer corresponds to a 2DTI system with helical edge states weakly coupled to its counterpart in adjacent layers. Here, the edge states are obtained from a tight-binding model of the stacking, fitted to match the first-principles data. These features allow us to propose a novel electric field control of the layer localization of the edge states. For an infinite AA stacking, in the bulk limit, the edge states from each layer combine to form a topological V-gutter-like surface states dispersion along the stacking direction. These surface states are highly anisotropic and tunable, introducing a new feature to the design of future multifunctional spintronic devices [36; 37].

II Results and Discussions

II.1 Stacking geometry

The Nickel bis(dithiolene) monolayer presents a planar geometry, composed by Ni atoms connected by carbon-sulfur molecules (C6S6), which combine to form a kagome lattice [dashed lines in Fig. 1(a)]. In order to analyze the energetic stability of multilayered (NiC4S${}_{4})_{3}$ systems, we have considered various stacking geometries. For instance, in (NiC4S${}_{4})_{3}$ bilayers (2ML), the stacking geometry is characterized by the alignment of sites X and Y of the different layers (2ML-XY), with X and Y = A, B, G, and H, as indicated in Fig. 1(a). Within this nomenclature, 2ML-AA [Fig. 1(b)] represents the eclipsed structure, 2ML-AB [Fig. 1(c)] dislocates one layer along the zigzag direction, and 2ML-GH [Fig. 1(d)] mimics the Bernal stacking of graphene.

In a previous study [26], we have found the 2ML-AA configuration to be the most stable. There, the calculations were performed using the non-local vdW-DF2 [38] correction to describe the vdW interactions. Indeed, the energetic preference of 2ML-AA was also verified in Refs. [39; 40] using different approaches to describe the long range vdW interactions [41; 42; 43]. In contrast, while Ref. [35] shows an energetic preference for the 2ML-GH configuration using the vdW-D3 [44], our results with this functional corroborates the previous works, showing a preference for the 2ML-AA case.

Notice that, despite the weak chemical interactions (essentially vdW) between the stacked layers, the electronic band structure of multilayered (NiC4S${}_{4})_{3}$ systems strongly depends on the stacking geometry. Indeed, the 2ML-AA and -GH configurations present significantly different band structures near the Fermi level, see Appendix A. Moreover, here we are interested in the topological insulator properties of (NiC4S${}_{4})_{3}$ MOFs [25], and such a dependence will reflect on the characteristics of the edge states in bilayer [26], and multilayer (NiC4S${}_{4})_{3}$.

A detailed analysis of the binding energy $E^{b}$ 111The binding energy ($E^{b}$) of (NiC4S${}_{4})_{3}$ bilayer systems is defined as, $E^{b}=E^{\rm ML}-\frac{1}{2}E^{\text{\rm 2ML-XY}},$ where $E^{\text{\rm 2ML-XY}}$ and $E^{\rm ML}$ are the total energies of the (NiC4S${}_{4})_{3}$ 2ML-XY, and (free standing) monolayer. Similar definition was applied for the 3ML-XYZ and bulk-XY systems. and geometric properties of different stackings is shown in Table 1. Here, the total energy calculations were performed considering full relaxation of the atomic positions and lattice vectors, until the force in each atom was less than 5 meV/Å. Moreover, we compare results obtained using different functionals, viz.: LDA [46], GGA-PBE [47], local vdW-D3 correction [44], and the non-local vdW-DF2 correction [38]. The wave functions were expanded in a plane-wave basis with cutoff energy of 450 eV. For BZ integration a $4\times 4\times 6$ ($4\times 4\times 1$) k-mesh for bulk (finite stacking) systems was considered. We have also considered the GGA+U method to correct the on-site repulsion of Ni d-orbitals. As shown in Appendix B, these corrections do not affect significantly our results.

Starting with the 2ML case, the PBE results show larger interlayer distances $d$, and lower values of binding energies when compared with the other functionals, thus indicating that the vdW interactions play an important role on the energetic stability and the equilibrium geometry of the multilayered (NiC4S${}_{4})_{3}$ systems. When vdW corrections are included, we find that 2ML-AA stackings are more stable than 2ML-GH 222Although our calculations have more base functions and denser k-grid than that of Ref. [35], we have also run calculations with exact the same parameters listed in this reference, but even in this case we find that the 2ML-AA stacking is more stable., with a $E^{b}$ difference of $\sim$0.15 eV/layer. The 2ML-AA stacking shows $E^{b}=$ 0.80 (0.86) eV/layer, and interlayer distance $d=$ 3.64 (3.54) Å, within the vdW-DF2 (vdW-D3) functional. In general, the bond distances obtained from the vdW-DF2 approach are larger than those obtained by the vdW-D3 approximation. Nonetheless, both vdW corrections capture the enhancement of the vertical corrugation ($\delta z$) due to the reduction of the interlayer distance in 2ML-AB, $\delta z$=0.1$\rightarrow$0.21 and 0.1$\rightarrow$0.15 Å, for vdW-DF2 and -D3, respectively. The energetic preference for the 2ML-AA configurations has been also verified within the LDA approximation. However, a comparison between the LDA and vdW results show that, within the LDA approach, the $E^{b}$ difference between 2ML-AA and -AB stackings is larger (0.62 eV/layer), and the 2ML-AA configuration shows a smaller interlayer distance (2.96 Å), which leads to a larger corrugation $\delta z$=0.08 Å.

For trilayers and the bulk systems, further binding energy results, shown in Table 1, also indicate the energetic preference for the AA stacking (3ML-AAA and bulk-AA). In comparison with the 2ML-AA case, the binding energy of the bulk-AA system increases by $\sim$1 eV/layer, followed by a reduction of the interlayer distance by 0.08 Å, $d$=3.64$\rightarrow$3.56 Å, (0.19 Å, 3.54$\rightarrow$3.35 Å) using the vdW-DF2 (vdW-D3) approximation. Within the LDA approach, we have also found the bulk-AA stacking to be the most stable. However, in contrast with the (vdW) results above, at the equilibrium geometry the interlayer distance increases by 0.16 Å when compared with the bilayer 2ML-AA, 2.96$\rightarrow$3.12 Å. These findings suggest that the vdW and LDA approximations present the same qualitative picture regarding the energetic stability of multilayered (NiC4S${}_{4})_{3}$ systems. However, the LDA approach fails to describe their equilibrium geometries. On the other hand, the PBE functional fail to describe the interlayer interaction, with the interlayer distance overestimated and smaller binding energies. Hereafter we will consider the most stable eclipsed (AA) geometry, obtained using the vdW-DF2 approach.

II.2 Electronic properties

To analyze the electronic structure of the AA stacking, we first consider the DFT results without SOC, as shown in Fig. 2 for 1ML, 2ML, 3ML, and bulk (NiC4S${}_{4})_{3}$ systems. Later on, in Section II.3, we will introduce the SOC to discuss the topological edge states.

For the single layer (NiC4S${}_{4})_{3}$, the typical kagome bands (KBs), indicated by red lines in Fig. 2(a1), are characterized by two energy bands [$E^{1}(\bm{k})$ and $E^{2}(\bm{k})$] with linear dispersion at the K-point, yielding a Dirac cone, connected to a nearly flat band [$E^{3}(\bm{k})$] at the $\Gamma$-point. The orbital projection, Fig. 2(b1), shows that the KBs are mostly composed by C-pz and S-pz orbitals, while near the center of the Brilloin zone, and along the M-K direction, the Ni-dxz and -dyz orbitals also contribute to the formation of those energy bands. For the 2ML-AA case, one finds two sets of KBs, Fig. 2(a2). There, one set of KBs is fully unoccupied, while the other is partially occupied, with the Dirac point lying slightly below the Fermi level. Such a metallic character is also seen in 3ML, Fig. 2(a3), where three sets of KBs are seen. The orbital projection on these energy bands, Figs. 2(b2) and 2(b3), reveal that, despite showing similar energy dispersion, the KBs present distinct orbital distribution. Namely, the contribution of the Ni-dxz/dyz [blue circles in Fig. 2(b1)-(b3)] orbitals are larger on high energy KBs, whereas the C-pz and S-pz orbitals [green circles in Fig. 2(b1)-(b3)] are dominant on the low energy KBs. Such an orbital distribution has been strengthened in the (NiC4S${}_{4})_{3}$ bulk system. As shown in Figs. 2(a4) and 2(b4), the contribution of the Ni-d (C-pz and S-pz) orbitals becomes dominant on the fully occupied (empty) KBs at high (low) energy, which are localized along the A-L-H-A ($\Gamma$-M-K-$\Gamma$) directions, see inset of Fig. 2(a4).

II.2.1 Effective tight-binding model

Within the eclipsed AA stacking geometry, the interlayer coupling is mostly given by weak vdW interactions between kagome sites aligned along the stacking direction $z$. Given this scenario, we propose a tight-binding (TB) model, for the stacking of $N$ layers, considering first neighbors coupling between monolayers. The resulting Hamiltonian reads

[TABLE]

where $\ket{j}$ labels the subspace of each layer $j$, $\mathbb{I}_{N}$ is a $N\times N$ identity matrix within this subspace, and ${\bm{k}}=(k_{x},\,k_{y})$ is the in-plane momentum. The $N$-layer Hamiltonian $H_{N}$ above is written in terms of the monolayer’s $3\times 3$ Hamiltonian $h_{3}({\bm{k}})$, which gives the characteristic KBs dispersion. The vdW interlayer coupling is set by $\Delta$, for which the diagonal form ($\mathbb{I}_{3}$) labels the direct coupling between equivalent atoms in the AA stacking geometry. The second term of $H_{N}$ is a tridiagonal Toeplitz matrix [49], yielding a closed form for the eigenvalues as

[TABLE]

with $E^{n}(\bm{k})$ being the kagome monolayer dispersion [Fig.2(a1)], where $n$ labels each of the three KBs for each set $\nu=1,\dots,N$. The eigenvectors are $\ket{n,\,{\bm{k}}}\otimes\ket{\nu}$, where the first term is an eigenstate of $h_{3}(\bm{k})$. To understand the second quantum number $\nu$ and the eigenstate component $\ket{\nu}$, it is interesting to notice that the tridiagonal Toeplitz form of the $\Delta$ coupling is equivalent to the Hamiltonian of a scalar $N$-sites chain, i.e., a quantum well along the stacking direction. Namely, $\ket{\nu}=\Omega\,[\phi_{1,\nu},\cdots,\phi_{j,\nu},\cdots,\phi_{N,\nu}]$, where $\Omega=\sqrt{2/(N+1)}$ is a normalization factor, and

[TABLE]

is the contribution of the state $\nu$ to the $j$-th layer. In the second form of $\phi_{j,\nu}$ above, we use the quantum well analogy to introduce $k_{\nu}=\nu\pi/L_{N}$ as the quantized $k_{z}$ momentum along the stacking direction $z_{j}=jd$, where $d$ is the interlayer distance, and $L_{N}=(N+1)d$ is the effective stacking length 333In this analogy, the quantum well chain has $N+2$ sites distant by $d$, yielding a length $L_{N}=(N+1)d$. A hard-wall boundary condition makes $\phi_{0,\nu}=\phi_{N+1,\nu}=0$, such that the real layers lie within $1\leq j\leq N$.. This analogy is useful to interpret the layer localization and the electric field effects below.

To obtain an estimate for $\Delta$, we compare the first-principles band structures with the eigenvalues $\lambda^{n}_{\nu}({\bm{k}})$ of our TB model. The first-principles energy dispersion for 2ML and 3ML, Fig. 2, show that each set of KBs are nearly rigidly shifted from one another. This matches the results of the TB model, where $\lambda^{n}_{\nu+1}({\bm{k}})-\lambda^{n}_{\nu}({\bm{k}})$ is independent of $\bm{k}$. For 2ML [Fig. 2(a2)], the two sets ($\nu=1$ and 2) of KBs are distant by $\lambda^{n}_{2}({\bm{k}})-\lambda^{n}_{1}({\bm{k}})=2|\Delta|\approx 0.56$ eV, thus $|\Delta|\approx 0.28$ eV. For the 3ML case we also find $|\Delta|\approx 0.28$ eV.

Taking the bulk limit of $N\rightarrow\infty$ layers, we can explore the quantum well analogy to express $\lambda^{n}_{\nu}({\bm{k}})\rightarrow\lambda^{n}_{\nu}({\bm{k}},k_{z})=E^{n}({\bm{k}})-2\,|\Delta|\,\cos(k_{z}d)$, assuming that now $k_{z}$ is a continuous variable, since $L_{N}\rightarrow\infty$. The $\cos(k_{z}d)$ dispersion is clearly seen in Fig. 2(a4) connecting the topmost and bottommost KBs along the $\Gamma$-A direction ($\parallel k_{z}$). From this dispersion, we obtain $|\Delta|=0.31$ eV. Such an increase of $\Delta$ (0.28 $\rightarrow$ 0.31 eV), although small, is a result of the reduction of the interlayer distance in the bulk system with respect to the few layer systems, 2ML or 3ML in Table 1.

To further explore the band width changes, and support our TB model, in Fig. 3(a) we show $|\Delta|$ as a function of the interlayer spacing $d$ for the (NiC4S${}_{4})_{3}$ bulk system. This distance can be controlled by an external pressure (uniaxial strain), thus allowing us to tune $\Delta$. Since the $\Delta$ coupling reflects the hybridization of neighboring layers, it decreases as $d$ increases. As discussed above, $\Delta$ is extracted from the bands connecting KBs along the $\Gamma$-A direction. At the equilibrium interlayer distance, $d$=3.56 Å, the energy difference between the upper- and bottommost KBs is about 1.24 eV, Fig. 3(b). Meanwhile, for $d$=3.63 (3.49) Å, this energy split reduces (increases) to 1.12 (1.44) eV [Figs. 3(c) and 3(d)]. That is, a $\sim 2$% uniaxial strain induces a $\sim 10$% change in $\Delta$.

II.2.2 Layer localization

The quantum well analogy, introduced above, enricher our discussion in the layer localization of the kagome states for finite AA stackings. First, in order to verify the validity of this model, we performed first-principles DFT calculations of the KBs for a 5ML slab. Panels L1 to L5 in Fig. 4(a) show the band structures projected on each layer $j$, and indicates the sets $\nu$ of KBs for clarity. Even though the band structure is dense at low energies, it is possible to see that: (i) the highest and the lowest KBs (${\nu\text{=1, and 5}}$) are mostly localized in the inner layers (L2, L3, L4); (ii) L3 does not contribute to ${\nu\text{=2, and 4}}$, which localizes mostly on the outer layers; (iii) the $\nu=3$ KBs contributes mostly to L1, L3 and L5. These features are more clearly seen in the circles of the left panel in Fig. 4(b), which compares the layer projections $|\phi_{j,\nu}|^{2}$ extracted from the DFT data [Fig. 4(a)], with the TB solutions. Indeed, these projections match remarkably well, thus validating the TB model.

In the presence of an external perpendicular electric field $\bm{E}_{ext}=E_{ext}\hat{z}$, the projections $|\phi_{j,\nu}|^{2}$ tend to shift to the outer layers. This is shown in Fig. 4(b) and (c), for a 5ML and 10ML slab respectively. There, we compare the eigenstate $\phi_{j,\nu}$ distributions among the layers in the absence (left panels) and presence (right panels) of $E_{ext}$. For a 5ML slab, the DFT data on the right panel of Fig. 4(b) is obtained for a nominal $E_{ext}=0.5$ eV/Å. However, due to charge rearrangement between the layers, the actual electric field is inhomogeneous across the layers. Nonetheless, the DFT data matches qualitatively well with the TB model for an effective uniform field $E_{ext}=0.2\Delta/d\approx 0.02$ eV/Å. Thus, the localization of each state follows the on-site energy shift due to $E_{ext}$, e.g., the lower energy KBs shifts towards L1.

The same behavior is extrapolated to thicker slabs, for instance a 10ML slab as shown in Fig. 4(c). Therefore $E_{ext}$ allows one to control the layer-localization of the KBs projection. In turn, combining this feature with the tunability of the Fermi energy (e.g., doping), one obtains an interesting engineering of the layer localization of topological edge states, which we discuss in the next section.

II.3 Edge states

The SOC opens an energy gap at the degenerated points ($\Gamma$ and K) of each set of KBs, yielding a QSH phase [25; 26]. Interestingly, the weak vdW coupling between the layers of our (NiC4S${}_{4})_{3}$ AA stacking leads to a series of nearly independent edge states in each layer. To model this layered QSH phase, we add the SOC to $h_{3}(\bm{k})$ of our TB model, since the large number of atoms involved in ribbon geometries makes first-principles calculations unattainable. The on-site, hopping, and SOC parameters within the single layer, and the $\Delta$ coupling between layers were obtained by fitting the first-principles DFT band structure of the bulk system. Further details can be found in the Appendix C and in Refs.  [26; 12].

II.3.1 Edge states on a finite stacking

In this section we consider a 10ML AA stacking in a ribbon geometry. Nevertheless, the following discussion applies to any finite stacking. First, in the absence of electric fields, we calculate the spin-polarized edge states and present its projections into the topmost (L10), middle (L5), and bottommost (L1) layers in Figs. 5(b1)–5(b3). These projections agree with the layer resolved distribution of the KBs discussed previously. For instance, the $\nu=1$ state (low energy) is mostly located in the middle layer. The same is true for the $\nu=10$ (high energy) state, since their KBs projections $\phi_{j,\nu}$ are similar, see left panel of Fig. 4(c). On the other hand, states within $2\leq\nu\leq 9$ are spread along many layers, as illustrated in Fig. 5(a). These results tell us that the vdW coupling between layers reflects into weakly coupled QSH edge states. Such edge states distribution allows for an unexploited control of the edge currents. For instance, tuning the Fermi level to the $\nu=1$ state, the helical currents will flow mostly through middle of the slab, while for $\nu=2$ state is distributed along its extremities, see top panel of Fig. 5(a).

On the other hand, in the presence of $E_{ext}$, the layer localization of the edge states follows that of the KBs, see right panel of Fig. 4(c). This allows for a fine tuning of the edge states contributions near the Fermi level. To illustrate, Figs. 5(c1)-(c3) show the edge state dispersions for an external field $E_{ext}=0.5\Delta/d_{z}\approx 0.05$ eV/Å. Here, we can see the lower (higher) energy edge states beginning to localize in the bottom (top) layer. Increasing the external field, $E_{ext}=0.5\rightarrow 1.0\Delta/d_{z}$, the localization effect is enhanced, compare Figs, 5(c3) and (d3). That is, the field changes the on-site energies of each layer, pushing the $\nu=1$ state towards the bottom layer, as depicted in Fig. 5(a). Consequently, the states near the Fermi level [shaded area in Figs. 5(b)-(d).] are mostly located in the bottom layer. If the field is reversed, $E_{ext}\rightarrow-E_{ext}$, the states near the Fermi level will switch the localization to the top layer. This tuning of the edge states, by $E_{ext}$ and/or changing the Fermi level, allow for an additional degree of freedom that can be useful to design electronic (layertronic) devices based on topologically protected edge currents in metal organic frameworks.

II.3.2 Edge states on the bulk stacking

The electronic confinement effects, present in the few layers systems, no longer takes place in the bulk phase. Here, the weakly coupled edge states from each layer are evenly spread along the side surfaces of the bulk staking, as shown schematically in Fig. 6(a1). Nonetheless, these side surface states carry the characteristics of the individual QSH edge states of the monolayers.

As discussed previously, the infinite superposition of the KBs gives rise to an continuum energy dispersion $\propto-2|\Delta|\cos(k_{z}d)$, where $k_{z}$ is the wave vector parallel to the stacking direction $\Gamma$-A. This bulk dispersion is projected in the 2D BZ of the side surface [Fig. 6(a2) and 6(b)], yielding a continuum of linear crossings at the edge of the Brillouin zone (X–M $\parallel k_{z}$), as shown in Fig. 6(c).

Interestingly, this peculiar V-gutter-like edge state dispersion [Fig. 6(c)] is highly anisotropic. An effective model for these states can be cast as a individual 2DTI edge states coupled along the stacking direction, wich reads

[TABLE]

where $\sigma_{0}$ is the $2\times 2$ identity matrix, $\sigma_{z}$ is third Pauli matrix in spin space, $v_{F}$ is the velocity of the edge states along $k_{x}$, and $\Delta$ is the same interlayer coupling from the bulk model. Along $k_{x}$, the states are helical with constant velocities $v_{x}=\pm v_{F}$ for spin up and down, while along $k_{z}$ the velocity is spin independent, but changes with $k_{z}$ or Fermi energy, $v_{z}=(|\Delta|d/\hbar)\sin(k_{z}d)$. This anisotropy is a consequence of the weak vdW coupling between the stacked layers. Thus, as shown in Fig. 3, the dynamics along $z$ can be further controlled by changing $\Delta$ with the uniaxial strain.

III Conclusion

In this paper, we show that electronic layer localization can be engineered in multilayer MOFs by controlling the Fermi level (e.g., doping) and an external electric field $E_{ext}$. This is possible due to the van der Waals (vdW) coupling between adjacent layers on eclipsed AA stackings. Here, we consider the (NiC4S${}_{4})_{3}$ MOFs, showing that AA stacking is the most stable for this material. Nonetheless, our results remain valid for other MOFs that also allow for eclipsed stackings. The interlayer vdW coupling reflects into weakly coupled topological edge states in each monolayer of the stackings. Consequently, we have shown that the layer localization of edge states can be controlled by $E_{ext}$, applied along the stacking direction, thus providing a novel degree of freedom to the design of TI based electronics.

Furthermore, extrapolating the system to an infinite stacking of nanoribbons, the weak vdW coupling between layers leads to a new V-gutter-like Dirac states in the side surfaces. This new surface states are highly anisotropic, showing a trivial velocity along the stacking direction, while along the monolayer edge a helical topological current is always present for a wide range of energy.

IV ACKNOWLEDGMENTS

The authors acknowledge financial support from the Brazilian agencies CNPq, and FAPEMIG, and the CENAPAD-SP and Laboratório Nacional de Computação Científica (LNCC-SCAFMat) for computer time.

Appendix A Band structure for bilayer stacking

For the (NiC4S${}_{4})_{3}$ bilayer, the stacking order has great impact in the electronic properties. For instance, in the mirror symmetric AA stacking, Fig. 2(a2), we find sets of the characteristc kagome band dispersion. On the other hand, for the non-mirror symmetric stacking, Fig. 7(a) and (b), the coupling between the layers leads to a different hybridization between the states, where the characteristic kagome bands are no longer present.

Appendix B DFT+U

We have performed GGA+U calculations [51] in order to verify the effects of the Hubbard on-site repulsion term for the d-orbitals of Ni in (NiC4S${}_{4})_{3}$. We find that our main results are not affected by this correction. The comparison between the $U=0$ eV and $U=3$ eV cases in Fig. 8(a) and (b) shows that the Ni d orbitals goes slightly down in energy for a finite $U$, which gives a negligible change in the charge state of Ni atoms, going from $+0.66\,e$ for $U=0$ to $+0.61\,e$ for $U=3$ eV. Nevertheless, Fig. 8 shows that the Kagome band dispersions are not significantly affected by $U$. Since this is the main characteristic that leads to the tunability of the layer contributions (the sets of kagome bands), we can conclude that our results are robust, and the Hubbard repulsion term can be neglected here.

Appendix C TB model parameters from DFT

The considered tight-binding (TB) Hamiltonian is given by

[TABLE]

where $H_{0}$ contains the on-site energy $\epsilon_{0}$ and hoppings within each layer $L$,

[TABLE]

with $i$ and $j$ indicating sites in the same layer, and sum over $L$ runs over all layers. $H_{SO}$ is the intrinsic spin-orbit coupling,

[TABLE]

The last term, $H_{L}$ contains the coupling between layers,

[TABLE]

where $k$ and $l$ indicates the same site, but in different layers. Here, $c_{i}^{\dagger}=(c_{i\uparrow}^{\dagger},c_{i\downarrow}^{\dagger})^{T}$, $c_{i}=(c_{i\uparrow},ci_{i\downarrow})^{T}$, $c_{i\chi}^{\dagger}$ ($c_{i\chi}$) is the creation (annihilation) operators of an electron in the $i$-th site with spin $\chi$, and ${\bm{\sigma}}$ are the spin Pauli matrices. The ${\bm{d}}_{ij}$ vectors connect the $i$-th to the $j$-th site. For the hoppings ($t_{ij}$ and $t_{ij}^{z}$) and spin-orbit couplings ($\lambda_{ij}$), we consider a distance dependent parametrization, which allow us to control the hopping range. This parametrization reads

[TABLE]

[TABLE]

and

[TABLE]

All hoppings are normalized in relation to the first neighbor hopping, by the factor $N=\exp(\alpha\,d_{nn})$ and $N_{z}=\exp(\alpha_{z}\,d_{z})$, where $d_{z}$ is the interlayer distance. In Fig. 9 we show a comparison between the DFT and TB model, using the parameters shown in Tab. 2. Note that the TB parameter $\Delta$ is consistent with the one obtained through the analytic model.