Superlattices of Gadolinium and Bismuth Based Thallium Dichalcogenides as Potential Magnetic Topological Insulators

Using relativistic spin-polarized density functional theory calculations we investigate magnetism, electronic structure and topology of the ternary thallium gadolinium dichalcogenides TlGdZ2 (Z= Se and Te) as well as superlattices on their basis. We find TlGdZ2 to have an antiferromagnetic exchange coupling both within and between the Gd layers, which leads to frustration and a complex magnetic structure. The electronic structure calculations reveal both TlGdSe2 and TlGdTe2 to be topologically trivial semiconductors. However, as we show further, a three-dimensional (3D) magnetic topological insulator (TI) state can potentially be achieved by constructing superlattices of the TlGdZ2/(TlBiZ2)n type, in which structural units of TlGdZ2 are alternated with those of the isomorphic TlBiZ2 compounds, known to be non-magnetic 3D TIs. Our results suggest a new approach for achieving 3D magnetic TI phases in such superlattices which is applicable to a large family of thallium rare-earth dichalcogenides and is expected to yield a fertile and tunable playground for exotic topological physics.

Intrinsic MTIs and heterostructures on their basis have emerged recently as a promising alternative to observe the above-mentioned quantized topological effects at elevated temperatures [2,3,. Discovery of the first representative of the MTI class, the antiferromagnetic (AFM) TI [50] MnBi 2 Te 4 [2,[34][35][36][37], has led to an appearance of an entire family Gd-Te-and -Tl-Te-Bi-Te-units, creates a desired SOC-driven bulk band gap inversion, making these systems potential 3D MTIs both intrinsically and in the forced FM state, which can be realized under an external magnetic field. Taking into account a great variety of the TlGdZ 2 -like compounds, selenides, tellurides and even sulfides, containing instead of Gd other rare earth (RE) elements, such as Pr, Nd, Sm, Tb, Dy, Ho, Er, and Tm [71,72,77], our study uncovers a new large family of potential magnetic topological materials of the (TlREZ 2 ) m /(TlXZ 2 ) n type (X = Bi, Sb and Z = S, Se, Te). In these systems, by varying RE, X, Z, and m, n-parameters, it should be possible to tune not only the exchange coupling (both intra-and interlayer), but also the magnetic anisotropy and, hence, the overall spin structure. Moreover, the electronic properties such as the fundamental band gap size and character (inverted or not), 4 f levels position, as well as the SOC strength can be varied as well. Thus, the TlREZ 2 /TlXZ 2 superlattices are expected to provide a fertile playground for the realization of new and exotic physics. In a broader sense, our results suggest an alternative approach for achieving a 3D MTI state. Namely, apart from doping non-magnetic TIs with transition metal elements or looking for a material combining magnetism and topology intrinsically, these properties can be realized in a superlattice of the isostructural stoichiometric magnetic trivial and non-magnetic topological insulators.

Calculation Details
The calculations were performed within the DFT using the projector augmented-wave method [78], implemented in VASP [79,80]. The generalized gradient approximation (GGA) was used to describe the exchange-correlation energy [81]. The Hamiltonian contained scalar-relativistic corrections and the SOC was taken into account by the second variation method [82]. We used the GGA + U approach [83,84] to describe the strongly localized Gd 4 f states. The values U = 6.7 eV and J = 0.7 eV were used, typically yielding a reasonable description of the electronic structure of the Gd-containing compounds [85,86]. The electronic structures, calculated with VASP, were verified against those obtained within the full-potential linearized augmented plane wave (FLAPW) method [87] within the GGA+U approach [85,88] under the fully localized limit [89] as implemented in FLEUR [90] and a good agreement was found.
For the TlGdZ 2 systems, the full crystal structure optimization, i.e., that of the a and c lattice constants, c/a ratio, and atomic coordinates, was performed in the rhombohedral unit cells, while for the superlattices the hexagonal unit cells were used. The lattice parameters obtained for the bulk TlXZ 2 (Table 1) are in good agreement with the experimental ones [71][72][73]75], being slightly larger than the latter, within 1.5 %. The intralayer exchange coupling was investigated using the ( √ 3 × √ 3R30 • ) × 1 hexagonal cells containing 3 atoms per layer, which is needed to model the intralayer non-collinear 120 • AFM structure. To study the interlayer exchange coupling, we used oblique cells having a hexagonal ab basis with the (1 × 1) in-plane periodicity and containing 8, 16, and 24 atoms in the cases of TlGdZ 2 , TlGdTe 2 /(TlBiTe 2 ) 1 , and TlGdTe 2 /(TlBiTe 2 ) 2 , respectively. Atomic positions were relaxed using a force tolerance criterion of 10 meV/Å. Surface band structures were calculated within a tight-binding approach based on the maximally localized Wannier functions [91,92] using iterative Green function method [93] implemented in WannierTools package [94]. Z 2 topological invariant ν 0 was calculated according to [95]: where n(Γ i ) ± is a number of occupied states with expectation values of parity ±1 at inversion symmetry invariant momenta points Γ i .

Magnetic State of TlGdZ 2
TlGdZ 2 has a markedly anisotropic crystal structure, in which the distance between Gd atoms in the basal ab plane is significantly smaller than that along the c-axis. In such quasi-2D magnets, there is a hierarchy of exchange interactions: the intralayer couplings are usually significantly stronger than the interlayer ones. We scrutinize these interactions by calculating the total energies of the following three spin configurations. In the first two, the interaction between the nearest neighbors in the Gd layers is assumed to be ferromagnetic, while their alignment along the c-axis can be either FM (Figure 1c) or AFM (Figure 1d), enabling evaluation of the interlayer exchange coupling. Determination of the intralayer coupling requires consideration of an additional spin structure, that could occur if the interaction between the nearest neighbors in the Gd layers was antiferromagnetic. It is wellknown that in combination with a 2D hexagonal lattice, AFM exchange coupling leads to a frustration [97], typically resolving itself into a non-collinear pattern, in which three spin sublattices form an angle of 120 • with respect to each other. Such a spin structure is shown in Figure 1e and will be hereinafter referred to as non-collinear antiferromagnetic (NCAFM).
While the character and strength of the interlayer exchange coupling can be conveniently characterized by ∆E ⊥ = E AFM − E FM , a simple difference between E NCAFM and either of E AFM and E FM does not describe the intralayer coupling precisely. Indeed, apart from a dominant contribution from the intralayer alignment, the total energies of the cases with FM Gd planes (i.e., E AFM and E FM ) also include a smaller one, coming from the interlayer exchange coupling, which is either purely AFM or FM. However, in the NCAFM case, there can be no purely AFM or FM interlayer alignment. This is because of the way this spin structure combines with the ABCABC stacking of atomic layers in TlGdZ 2 (Figure 1b,e), leading to frustration. In detail, a Gd magnetic moment from one layer, that has three nearest magnetic neighbors in each of the two adjacent Gd layers (Figure 1b), cannot simultaneously be either AFM-or FM-coupled to all of them because we necessarily assume all of the Gd layers to be NCAFM (which is to model the intralayer AFM interaction case). Therefore, both E NCAFM − E FM and E NCAFM − E AFM inevitably contain a certain contribution from the coupling between layers, since it does not cancel out. Indeed, the energetical favorableness of a particular interlayer configuration (FM or AFM) makes the corresponding difference smaller and, vice versa, its unfavorableness makes it larger, which, however, has nothing to do with the intralayer coupling since E NCAFM stays the same. Let us therefore account for this spurious contribution as follows: where δ eliminates the contribution from the interlayer coupling. By doing simple math, we get According to our total-energy calculations, ∆E || < 0 for both TlGdSe 2 and TlGdTe 2 . In the former case, the energy gain, |∆E || |, amounts to 9.5 meV per Gd pair, while in the latter it is 11.7 meV (Table 1). We thus can conclude that there is a pronounced tendency towards the non-collinear 120 • AFM structure formation in the Gd layers of the TlGdZ 2 compounds. On the other hand, our calculations show that ∆E ⊥ is negative, too, and equal to −1.3 (−2.5) meV per Gd pair for TlGdSe 2 (TlGdTe 2 ), signaling the AFM character of the interaction between neighboring Gd planes. As expected, the interlayer coupling is significantly weaker than the intralayer one, which is due to the quasi-2D magnetic character of these compounds, determined by their layered crystal structure. Indeed, the interlayer Gd-Gd distance along c is significantly larger than the in-plane one, e.g., 8 [34,45,55], that are known to be in the magnetically three-dimensional regime. The Gd local magnetic moments in both TlGdSe 2 and TlGdTe 2 are roughly equal to 7 µ B (S = 7/2), in a reasonable agreement with the available experimental data [98] and DFT calculations [76].
As has been discussed above, due to the combination of the intralayer NCAFM state and the ABCABC-type stacking of atomic planes, the interlayer exchange coupling leads to magnetic frustration. To get a deeper insight into the influence of the interlayer coupling on the TlGdZ 2 magnetic structure, different spin alignments between the 120 • -ordered planes have been considered, as is further exemplified by TlGdTe 2 . We first note that in the NCAFM configuration considered each Gd moment has one ferromagnetically aligned nearest neighbor in each adjacent Gd layer, while the other two nearest neighbor moments are rotated by ±120 • with respect to it. The presence of the FM-coupled pairs might seem energetically unfavorable since the interlayer interaction tends to align local moments antiferromagnetically because ∆E ⊥ < 0. However, if we choose an antiparallel alignment in the same pairs, the total energy stays essentially the same (within the calculation accuracy). This is because the energy decrease due to the AFM alignment is compensated by the energy increase due to the interactions with the other two nearest neighbor moments, pointing at ±60 • . The sum of the latter two moments yields a vector of the same magnitude as the first one but points in the opposite direction.
Another possibility is that the interlayer configuration is helimagnetic, which is compatible with the interlayer magnetic frustration [99,100]. It would also be in line with the results of the electron paramagnetic resonance study of TlGdSe 2 [75], pointing towards a possible helimagnetic state. To explore this, we have considered a spin structure in which the 120 • pattern of the Gd local moments rotates by an angle φ = 30 • when going from one Gd layer to another along the c-axis. Here, each Gd moment is perpendicular to that of one of the nearest neighbors in each adjacent Gd layer, while with respect to the other two nearest neighbor moments it deviates by 30 • from the collinear FM and AFM alignments. In such a situation, no energy gain arises either, as compared to the NCAFM structure. Still, we can exclude neither a possibility of the helix angle being incommensurate with the hexagonal symmetry nor an appearance of a broken-helix magnetic order with more than one rotation angle [65]. Noteworthy, our magnetic anisotropy energy calculations show that the SOC does not favor any particular direction of the local magnetic moment, either in-plane or out-of-plane. In this situation, the role of the dipole-dipole interaction becomes crucial, which generally favors locking of the magnetic moments within the basal plane and realization of the tail-chase structure, although a canting towards the out-of-plane direction or the amplitude modulations can also occur in some cases [101][102][103]. It is thus evident that the exact 3D magnetic ground state structure of both TlGdTe 2 and TlGdSe 2 is very difficult to determine using the DFT calculations. The neutron diffraction measurements, which is a standard tool used to resolve magnetic structures in experiment [65,99], are desirable for TlGdZ 2 .

Electronic Structure of TlGdZ 2 (Z = Se, Te)
The bulk band structures of TlGdZ 2 are presented in Figure 2 (the NCAFM structure, shown in Figure 1e, is assumed). It can be seen that both TlGdSe 2 and TlGdTe 2 show insulating spectra with indirect band gaps. According to the density of states calculations, the gap sizes are 1.26 eV and 0.56 eV, respectively. In both cases, the valence band maxima are formed by p-states of a chalcogen, while the conduction band minima are composed of the Tl pand Gd d-states. No signature of the SOC-driven band gap inversion is observed from the analysis of the orbital composition of the states forming the band gap edges. We further confirm the absence of the inversion by direct calculations of the density of states performed at different values of the SOC constant. Thus, we can conclude that TlGdZ 2 (Z = Se, Te) are topologically trivial magnetic semiconductors.

TlGdTe 2 /(TlBiTe 2 ) n Superlattices
Taking advantage of a relatively small lattice mismatch (∆a ≈2%) between TlGdTe 2 and the non-magnetic TI TlBiTe 2 (Supplementary Figure S1), the similarity of their atomic composition, as well as their having the same R3m crystal structure (Figure 1a), we consider bulk superlattices of the TlGdTe 2 /(TlBiTe 2 ) n type with n =1 and 2. Their structure, shown in Figure 3a,d, comprises four-layer -Tl-Te-Gd-Te-units alternated with single (n = 1) or double (n = 2) four-layer unit of -Tl-Te-Bi-Te-. We have chosen TlGdTe 2 as opposed to TlGdSe 2 since the former has a smaller bulk band gap, which should be more easily inverted due to the strong SOC of Bi thus potentially allowing the realization of a topologically non-trivial magnetic insulating system.
After a full structural optimization (see lattice parameters in Table 2), the magnetism of TlGdTe 2 /(TlBiTe 2 ) n has been studied. On the one hand, we find that the intralayer NCAFM structure is by about 11.5 meV more favorable than the FM one for both n = 1 and n = 2, similar to pure TlGdTe 2 . On the other hand, the interlayer exchange coupling for n = 1 weakens by an order of magnitude as compared to the latter (∆E ⊥ = −0.19 meV [ Table 2)] vs. −2.5 meV [Table 1)]), which is due to the increase of the Gd-Gd interlayer distance. When going from n = 1 to n = 2 it weakens by another order of magnitude. Table 2. The same as Table 1, but for TlGdTe 2 /(TlBiTe 2 ) n . Note that for n = 1, the c parameter is approximately twice as large as the one for n = 2 because the hexagonal unit cell of the former contains 24 atoms as compared to 12 of the latter. We then study the electronic structure of the TlGdTe 2 /(TlBiTe 2 ) n superlattices, as previously assuming the NCAFM configuration. In Figure 3b,e, it is seen that, as compared to pure TlGdTe 2 , the superlattices have a direct band gap (located in the A-point) and the gap size is much smaller (cf. data in Tables 1 and 2). At that, the gap of the n = 2 system (141 meV), having a higher Bi content, is larger than that of n = 1 (86 meV). The latter fact, along with the direct character of the band gap, points towards band gap inversion due to SOC, which is enhanced due to the presence of Bi. To confirm the band gap inversion we have calculated its value as a function of the SOC strength, λ. The results, presented in Figure 3c,f, show that at λ/λ 0 ≈ 0.91 (λ/λ 0 ≈ 0.84) the gap of the n = 1 (n = 2) system is closed, while at other values of λ/λ 0 it is non-zero, meaning that at λ = λ 0 it is indeed inverted around the A-point, which indicates a potential MTI state of both TlGdTe 2 /(TlBiTe 2 ) n .

Superlattice
Since, similarly to pure TlGdZ 2 , we do not know the exact magnetic ground state of TlGdTe 2 /(TlBiTe 2 ) n because of the coexistence of frustrated exchange interactions with dipolar contributions to magnetic anisotropy, we cannot reliably establish their topological classification. It should, however, be noted, that even for a complex magnetic structure, there may be a symmetry, which could give rise to a topological classification. For example, it has recently been found by neutron diffraction that EuIn 2 As 2 exhibits a low-symmetry broken-helix antiferromagnetic order with two helix turn angles, which nevertheless supports a magnetic topological-crystalline axion insulator protected by the combination of a 180 • rotation and time-reversal symmetry [65].
Taking advantage of relatively weak exchange interactions in TlGdTe 2 /(TlBiTe 2 ) n , it should in principle be possible to impose a ferromagnetically polarized state by applying an external magnetic field, as in the compounds of the MnBi 2 Te 4 family [5,6,62]. In the artificial ferromagnetic state with a magnetization perpendicular to the Gd layers, the TlGdTe 2 /(TlBiTe 2 ) n are characterized by gapped spectra (not shown) with bulk band gaps of 69 and 131 meV for n = 1 and 2, respectively. By varying the SOC constant λ we find that the band gaps are inverted, too. Moreover, the Z 4 invariant calculations show that both systems are axion insulators (Z 4 = 2), having the quantized topological magnetoelectric response in the bulk and chiral hinge modes [7,96]. Finally, the calculated Z 2 invariant shows that in the paramagnetic state these systems are time-reversal symmetric strong 3D TIs.
To further illustrate the topologically non-trivial character of our superlattices, in Figure 4 we show the surface electronic structure of TlGdTe 2 /(TlBiTe 2 ) 1 in the paramagnetic (PM) and ferromagnetic (FM) states. In the PM state, a linearly dispersing gapless surface state is clearly seen within the fundamental bulk band gap (Figure 4a), as expected for the time-reversal symmetric strong 3D TI. In contrast, when the system is driven in the FM state (which can be achieved by the out-of-plane external magnetic field), the Dirac point splits, consistently with what has been predicted for some compounds of the (Mn(Bi,Sb) 2 Te 4 ) · n((Bi,Sb) 2 Te 3 ) family in the FM state [55,[104][105][106].  Finally, we would like to discuss the similarities and differences between the here discussed MTI superlattices TlGdTe 2 /(TlBiTe 2 ) n and intrinsic MTIs of the (MnBi 2 Te 4 )·n(Bi 2 Te 3 ) and (MnSb 2 Te 4 )·n(Sb 2 Te 3 ) families. Pretty much as the latter systems, composed of the building blocks of the magnetic (Mn(Bi,Sb) 2 Te 4 ) and nonmagnetic ((Bi,Sb) 2 Te 3 ) compounds, our superlattices are made of two stoichiometric systems, too. However, while (Mn(Bi,Sb) 2 Te 4 ) · n((Bi,Sb) 2 Te 3 ) are compounds themselves, TlGdTe 2 /(TlBiTe 2 ) n can hardly be considered as such. Indeed, (Mn(Bi,Sb) 2 Te 4 ) · n((Bi,Sb) 2 Te 3 ) can be grown in the bulk single crystal form, while applying the bulk growth strategy to the Tl(Gd,Bi)Z 2 system would, most likely, result in a disordered TlGd x Bi 1−x Z 2 solid solution (found to be paramagnetic down to 2 K at x = 0.1) [107]. In this sense, TlGdTe 2 /(TlBiTe 2 ) n are somewhat similar to delta-doped systems, such as digital magnetic alloys [108], representing transition metal monolayers embedded in a matrix of semiconductor material, studied previously in the field of diluted magnetic semiconductors. We propose that TlGdTe 2 /(TlBiTe 2 ) n can be synthesized using layer-by-layer growth techniques, such as molecular beam epitaxy or atomic layer deposition, that have proven successful for growing layered van der Waals and non-van der Waals materials, complex heterostructures and superlattices [32,56,[109][110][111].

Conclusions
In summary, using ab initio and tight-binding calculations we have studied magnetism and electronic structure of the TlGdZ 2 (Z = Se, Te) compounds as well as TlGdTe 2 /(TlBiTe 2 ) n superlattices. Our results suggest a complex magnetic ground state in all these systems, which is due to the combination of the antiferromagnetic exchange interaction within Gd layers (resolving itself into the non-collinear spin structure), the fcc-type layer stacking, and the presence of the interlayer exchange coupling. As far as the electronic structure is concerned, TlGdZ 2 appear to be topologically-trivial semiconductors. However, a magnetic topological insulator state can be induced by constructing superlattices of TlGdZ 2 with isostructural Tl-based non-magnetic TIs. In this way, we predict that TlGdTe 2 /(TlBiTe 2 ) n potentially hosts a magnetic topological insulator state both intrinsically and under an external magnetic field. Our study not only suggests a new approach for achieving magnetic topological states of matter, but also uncovers a new large family of the TlREZ 2 /(TlXZ 2 ) n superlattices (RE = Pr, Nd, Sm, Gd, Tb, Dy, Ho, Er, and Tm; X = Bi, Sb; Z = S, Se, Te) with tunable magnetic, electronic and topological properties.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

MTI
Magnetic topological insulator FM Ferromagnetic AFM Antiferromagnetic NCAFM Noncollinear antiferromagnetic QAHE quantum anomalous Hall effect TME Topological magnetoelectric effect DFT Density functional theory SOC Spin-orbit coupling