**Advances in Topological Materials**

Editor

**Artem Pronin**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Artem Pronin Universitaet Stuttgart Germany

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Crystals* (ISSN 2073-4352) (available at: https://www.mdpi.com/journal/crystals/special issues/Advances Topological Materials).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-4235-5 (Hbk) ISBN 978-3-0365-4236-2 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editor**

**Artem Pronin** graduated from the Moscow Physics and Engineering Institute, Russia, in 1994. His PhD research was conducted at the General Physics Institute of the Russian Academy of Sciences, Moscow. In 1998, he received a PhD degree in solid-state physics. Afterwards, he held positions at McMaster University, Canada; Leiden University, the Netherlands; the University of Leeds, UK; and the High Magnetic Field Laboratory in Dresden, Germany. Since 2015, he has been leading a group at the 1st Physics Institute of the University of Stuttgart, Germany. The group investigates the optical properties of different topological materials, particularly, of Weyl and Dirac semimetals.

### *Editorial* **Advances in Topological Materials**

**Artem V. Pronin**

1. Physikalisches Institut, Universität Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany; artem.pronin@pi1.physik.uni-stuttgart.de

Materials with electronic bands that possess nontrivial topology have remained a focal point of condensed matter physics since 2005, when topological insulators were theoretically discovered by Kane and Mele [1,2]. In parallel to this remarkable discovery, Haldane and Raghu [3] realized that topological phases are a universal phenomenon of waves in periodic media. Thus, topological concepts can also be applied, for example, to electromagnetic waves in photonic crystals [4], magnons in magnetic materials [5], and sound waves in different periodic structures [6]. This Special Issue of *Crystals* represents a collection of 11 papers devoted to different aspects of experimental and theoretical studies on topological materials.

Five papers from the Special Issue focus on the theory side. Gao and Wang [7] propose a new design for an ideal photonic Weyl metacrystal ("ideal" means that there are no additional states at the Weyl-node energy). The Weyl nodes of this metacrystal are stabilized by the screw rotation symmetry of space group 19. The authors argue that this design might be advantageous for further experimental studies of photonic Weyl materials. Cheng and Gao [8] study a non-interacting Λ/*V*-type dice model composed of three triangular sublattices. By considering certain nearest-neighbor and next-nearestneighbor hopping terms, as well as a quasi-staggered on-site potential, they acquire the full phase diagrams for different energy band fillings. They find abundant topologically nontrivial phases with different Chern numbers and a metallic phase in several regimes. Nikolaev et al. [9] study the influence of uniaxial deformation on the band structure and topological properties of the multifold semimetal CoSi with large topological charges. The **k***·***p** Hamiltonian, which takes the deformation into account, is constructed from symmetry considerations near the Γ and *R* points of the Brillouin zone. The transformation of the multifold band crossings into nodes of other types with different topological charges, their shift in energy and in reciprocal space, and the tilt of the dispersion around the nodes are studied in detail, depending on the direction of uniaxial deformation. Polatkan and Uykur [10] present a theoretical study of the band structure and optical conductivity for another multifold semimetal, PdGa. They identify several characteristic features in the optical conductivity and relate their origin to the band structure. Yaresko and Pronin [11] calculate the *ab*-plane optical conductivity of the Weyl semimetal TaP and compare it to the experimental data. Based on these calculations, they propose an explanation of the strong low-energy peak observed in the experimental spectra: this peak originates from transitions between the almost parallel non-degenerate electronic bands split by spin-orbit coupling.

The other papers in this Special Issue report experimental findings. Dally et al. [12] present small-angle inelastic neutron scattering measurements of Fe3Sn2. Fe3Sn2 has recently been discovered to host room temperature skyrmionic bubbles and is known to have competing magnetic exchange interactions, correlated electron behavior, weak magnetocrystalline anisotropy, and lattice anisotropy. The results of Dally et al. reveal that, at elevated temperatures, there is an absence of significant magnetocrystalline anisotropy and that the system behaves as a nearly ideal isotropic exchange interaction ferromagnet. Hatnean et al. [13] report on the growth of large high-quality Ce-substituted SmB6 crystals via the floating zone method. The topological properties of SmB6 are currently being intensively discussed in relation to Kondo physics. Hence, the investigation of substituted

**Citation:** Pronin, A.V. Advances in Topological Materials. *Crystals* **2021**, *11*, 680. https://doi.org/10.3390/ cryst11060680

Received: 10 June 2021 Accepted: 11 June 2021 Published: 14 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

SmB6 samples is of interest. The structural, magnetic and transport properties of single crystals with different Ce contents are investigated by Hatnean et al. using X-ray diffraction techniques, electrical resistivity and magnetization measurements. The authors find that the substitution of Sm with magnetic Ce does not lead to long-range magnetic ordering.

The remaining experimental reports focus on optics. Shuvaev et al. [14] present sub-terahertz measurements of the quantum anomalous Hall effect (QAHE). In the static regime, the QAHE is observed as a step in Hall resistivity. At optical frequencies, it is transformed into a step in the polarization rotation, with the size of this step being equal to the fine structure constant, α ≈ 1/137. The authors measure the polarization rotation in thin films of the topological insulator (Bi,Sb)2Te3 doped with Cr and observe the expected steps at temperatures below 20 K. However, due to material issues, the size of the steps only reaches up to 20% of the theoretical value (at 1.85 K). At millikelvin temperatures, full-size steps are anticipated. Kamenskyi et al. [15] perform magnetooptical measurements of the topological insulator Bi2Te3 in the terahertz frequency range in magnetic fields up to 10 T. They report on the observation of a cyclotron resonance mode and ascribe it to free bulk carriers. The width of the mode demonstrates a non-monotonous behavior in the magnetic field. The authors propose that the mode width is defined by two competing factors: impurity scattering and electron–phonon scattering, which exhibit opposite behaviors in applied magnetic fields. Another topological insulator, Bi2Te2Se, is investigated by Zhukova et al. [16] by mid- and near-infrared optical measurements. The optical conductivity of Bi2Te2Se is found to be dominated by bulk carriers and shows a linear-in-frequency increase at 0.5 to 0.8 eV. This linearity might be interpreted as a signature of the three-dimensional (bulk) Dirac bands; however, the band structure-based calculations performed by the authors show that transitions between bands with complex dispersions contribute instead to the inter-band optical conductivity at these frequencies and, hence, the observed linearity is accidental. These results warn against oversimplified interpretations of optical conductivity measurements in different Dirac materials. Finally, Schilling et al. [17] investigate the broadband optical conductivity of the two-dimensional Dirac material CaMnBi2. They find that both components of the intraband conductivity follow a universal power law as a function of frequency at low temperatures. This conductivity scaling differs from the standard Drude-like behavior and might point toward quantum criticality in this system.

Overall, this Special Issue represents a few recent developments in the broad and growing field of topological material studies.

**Acknowledgments:** I am grateful to Dancy Yu for her editorial assistance in the production of this Special Issue and to Ece Uykur for proofreading the draft of this paper.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


## *Article* **Infrared Optical Conductivity of Bulk Bi2Te2Se**

#### **Elena S. Zhukova 1,2,\*, Hongbin Zhang 3, Victor P. Martovitskiy 4, Yurii G. Selivanov 4, Boris P. Gorshunov 1,2 and Martin Dressel 1,2,\***


Received: 8 June 2020; Accepted: 24 June 2020; Published: 28 June 2020

**Abstract:** Mid- and near-infrared measurements reveal that the optical conductivity of the three-dimensional topological insulator, Bi2Te2Se, is dominated by bulk carriers and shows a linear-in-frequency increase at 0.5 to 0.8 eV. This linearity might be interpreted as a signature of three-dimensional (bulk) Dirac bands; however, band-structure calculations show that transitions between bands with complex dispersion contribute instead to the inter-band optical conductivity at these frequencies and, hence, the observed linearity is accidental. These results warn against the oversimplified interpretations of optical-conductivity measurements in different Dirac materials.

**Keywords:** topological insulators; optical conductivity; Dirac materials

#### **1. Introduction**

Spin-orbit coupling often leads to the formation of linear bands in solids. Electrons in such bands (the Dirac electrons) manifest themselves in special ways in different experiments [1–5]. One of these manifestations is in their optical response: the contribution of a *d*-dimensional Dirac band to the inter-band optical conductivity, which is calculated to follow a simple power–law frequency dependence [6,7]:

$$
\sigma(\omega) \propto \omega^{\mathrm{d}-2}.\tag{1}
$$

Such optical-conductivity behavior—unusual for conventional materials—has indeed been confirmed for (quasi)-2D electrons in graphene, graphite, and the line-node semimetal ZrSiS, where σ(ω) ≈ *const*(ω) was reported [8–10]. In turn, the 3D Dirac electrons in Dirac and Weyl semimetals, such as ZrTe5, Cd3As2, and TaAs, provide the inter-band optical conductivity to be proportional to frequency, σ(ω) ∝ ω [11–13]. The linearity in σ(ω) over a broad frequency range in a 3D electron system is often considered as a "smoking gun" for Dirac physics. For example, Timusk et al. [14] suggested the presence of 3D Dirac fermions in a number of quasicrystals, based entirely on the observation of a linear σ(ω) in these materials.

Besides, enormous efforts have been made to investigate the symmetry-protected surface states of topological insulators [2,3]. However, the dominant physics of the bulk often obscures the surface properties and hence is generally considered as an obstacle for experiments targeting the surface states. Achieving dissipationless surface spin currents may be of primary importance for potential applications of topological insulators, nevertheless, investigations into bulk electronic properties are essential for understanding the complete picture of the topological-state formation [15].

Our experiments reveal that the bulk optical conductivity of Bi2Te2Se follows a linear frequency dependence in an appreciably broad spectral range. Based on band-structure calculations, we argue that this linearity is not due to transitions within (a) particular 3D linear band(s), but instead a result of contributions from the transitions between the bands with complex dispersion.

#### **2. Materials and Methods**

Bi2Te2Se bulk crystals were synthesized by a modified Bridgman method [16]. Highly purified (99.9999%) elemental starting materials (Bi, Te, and Se) (Chimmed, Moscow, Russia) were loaded in quartz ampules inside an inert-gas glove box in the stoichiometric ratio 2:2:1. The sealed evacuated ampules were kept at 850 ◦C for 24 h with periodic stirring to ensure the homogeneity of the melt, followed by a cooldown to 520 ◦C with a rate of 5 ◦C/h. The crystals were then annealed at 520 ◦C for six days. The typical crystal sizes obtained in this way were in the centimeter range. The crystals were cut into appropriate pieces for X-ray, Hall, and optical measurements (and kept in vacuum until the measurements).

Utilizing an X'Pert Pro Extended MRD X-ray diffractometer (PANalytical, Almelo, the Netherlands) we have confirmed the high structural quality of the crystals, see Figure 1. The free-carrier concentration and mobility were measured in a standard Hall geometry. Indium-soldered contacts were applied to razor-cut Hall bars with typical dimensions of 2 <sup>×</sup> 0.5 <sup>×</sup> 0.2 mm3. For all samples, the conduction was by *n*-type carriers. The properties of the sample, used in our infrared studies, are listed in Table 1.

**Figure 1.** Bi2Te2Se X-ray diffraction pattern. Inset: Rocking curve for the (0 0 15) reflection peak.

**Table 1.** Room-temperature properties of the single-crystalline Bi2Te2Se sample used for the optical measurements. The mobility value is typical for the samples with such electron densities [17].


Optical reflectivity was measured from the (001) plane on freshly cleaved surfaces. The room-temperature experiments were performed in the mid- and near-infrared spectral ranges (600–8000 cm−1, 75 meV–1 eV) with a Bruker Vertex 80v Fourier-transform infrared spectrometer (Bruker Corporation, Billerica, MA, USA). Freshly evaporated gold mirrors served for reference measurements. We used unpolarized light, because Bi2Te2Se possesses C3 rotational symmetry along the [001] direction and hence the (001)-plane response, expressed via a second-rank tensor, such as optical conductivity, is isotropic.

#### **3. Results and Discussion**

In the top panel of Figure 2, we plot the raw reflectivity data recorded at 300 K. The reflectivity is very flat between 4000 and 8000 cm<sup>−</sup>1. In order to obtain the optical conductivity from the reflectivity data, we first tried to fit the measured spectra using a standard Drude–Lorentz procedure [18]. However, we found that such flat reflectivity is impossible to fit in an acceptable way with a physically meaningful number of Lorentzians. In an alternative approach, we used Kuzmenko's variational dielectric function method [19], which produces optical functions with an accuracy equivalent to Kramers–Kronig. For the sake of convenience, the variable part of the dielectric response function was described by a large number of Lorentzians. Justification and details of this approach can be found in [20]. Similar to the Kramers–Kronig analysis, this method gives less accurate results near the edges of the experimental window. Thus, the results below approximately 2000 cm−<sup>1</sup> and above 7000 cm−<sup>1</sup> cannot be considered as accurate.

**Figure 2.** Top panel: [001]-plane reflectivity of Bi2Te2Se at 300 K: measurements (black line) and fit (red line). Bottom panel: bulk optical conductivity (real part) of Bi2Te2Se, as obtained from the reflectivity fit (black straight line) and the inter-band portion of optical conductivity, computed from the band structure of Figure 3 at 0 K (dashed line), as detailed in the text. The thin orange line is to mimic a linear increase in frequency.

The real part of the optical conductivity obtained from this fit is plotted in the bottom panel of Figure 2. The eye-catching feature of the figure is the linear increase in σ(ω) at 4000 to 7000 cm−<sup>1</sup> (~0.5–0.8 eV).

Let us first argue that the observed optical conductivity originates from the bulk of Bi2Te2Se. In Bi2Te2Se, the surface Dirac point lies inside the bulk band gap [21,22] and metallic surface states have been experimentally confirmed [21–27]. Nevertheless, Bi2Te2Se samples usually possess a significant concentration of bulk charges due to the basically unavoidable presence of defects, the so-called self-doping [26–30]. This is also the case for our sample—its bulk carrier concentration is rather large, as shown in Table 1. Furthermore, the skin depth, calculated from the complex optical conductivity,

is above 30 nm at any measurement frequency, while the thickness of the topologically non-trivial surface layer is believed to be around 1 nm [3]. Hence, the response detected by our optical measurements is due to the bulk.

**Figure 3.** Band structure of Bi2Te2Se. Black dashed (red solid) horizontal line indicates the original (shifted) Fermi energy.

Let us also note that, at elevated temperatures, the optical detection of surface carriers in Bi2Te2Se, as well as in similar compounds, such as Bi2Te3 and Bi2Se3, remains so far elusive, while bulk carriers clearly manifest themselves in the optical response of Bi2Te2Se [16,28–30] and related compounds [31–35]. Reijnders et al. have reported on a mixed (surface plus bulk) optical response in Bi2Te2Se for low frequencies at temperatures below some 40 K [30]. However, at room temperature, as well as at frequencies above 2000 cm<sup>−</sup>1, their data are perfectly reconciled with entirely bulk response.

Coming back to the linear σ(ω), it is tempting to interpret it in terms of Equation (1), namely, as a signature of a 3D Dirac band (because our σ(ω) reflects the bulk response). Such a band, however, is not expected to appear in the bulk of Bi2Te2Se [36]. We would like to point out that all the available optical conductivity spectra (ours and those previously reported in [16,28–30]) are rather similar to each other, although the linearity of σ(ω) is most apparent in our data. The deviations between the data sets can be assigned, for example, to the abovementioned difference in the exact Fermi-level position in different samples of Bi2Te2Se. In order to check the origin of the linear frequency increase in σ(ω), we performed band-structure calculations for Bi2Te2Se and then calculated its inter-band optical conductivity.

The band-structure and optical-conductivity calculations were performed using the full potential linear augmented plane-wave method, as implemented in theWIEN2k code [37]. The exchange-correlation functional is parameterized using the GGA approximation [38]. The self-consistent charge-densities and optical-conductivity calculations were done with 400 and 2000 k-points in the whole Brillouin zone, respectively. The results of the calculations are shown in Figures 2 and 3. The obtained band structure is basically identical to the one reported in [36]. In order to be reconciled with the bulk electron concentration (the self-doping problem mentioned above), the Fermi level needs to be shifted upwards, as compared to the undoped situation, as shown in Figure 3. From the figure, it is apparent that there is no truly Dirac band in the bulk of Bi2Te2Se.

The calculated optical conductivity is shown as a dashed line in Figure 2. Taking into account the generally poor reproducibility of the experimental infrared optical conductivity by first-principles calculations (cf., e.g., in [39,40]), the agreement between theory and experiment can be considered as fairly good. Further, we should point out that the computed σ(ω) has no intra-band (free-carrier) contribution. Thus, it is not surprising that the low-frequency experimental σ(ω) is larger than the theoretical line. Additionally, the effect of temperature broadening is absent in the calculations. Such broadening would make the smooth step at around 3000 cm−<sup>1</sup> even broader [10].Taking into account the mentioned issues in the computations of σ(ω) is outside of our capacity and beyond the scope of the paper. The important result of our computations is that the linear σ(ω) is nicely reproduced at 4000 to 6000 cm−<sup>1</sup> (~0.5–0.75 eV). Thus, we can conclude that this linearity comes as a cumulative effect of transitions between the bands, which do not have a simple linear dispersion. We note that recent measurements of BaCoS2 and GdPtBi provide other examples of linear σ(ω) not due to a simple 3D Dirac band [41,42].

#### **4. Conclusions**

We have experimentally found that the bulk optical conductivity of Bi2Te2Se is linear in frequency at 4000 to 7000 cm−<sup>1</sup> (~0.5–0.8 eV). Our computations demonstrate that this linearity is not due to transitions within a 3D Dirac band, but emerges as a cumulative effect of transitions between the bands with complex dispersion. Obviously, similar situations can appear in other systems and, thus, suggestions for Dirac physics based on optical-conductivity measurements have to be made cautiously.

**Author Contributions:** Infrared measurements, E.S.Z. and B.P.G.; sample preparation, V.P.M. and Y.G.S.; data analysis, E.S.Z.; writing—original draft preparation, E.S.Z.; writing—review and editing, H.Z., M.D., and B.P.G.; project administration, B.P.G.; funding acquisition, Y.G.S. M.D., and B.P.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Ministry of Science and Higher Education of the Russian Federation (Program "5 top 100"), by the Russian Science Foundation (grant No. 17-12-01544) and by the Deutsche Forschungsgemeinschaft (DFG) via grants No. DR228/51 and No. ZH559/2-1.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Ideal Photonic Weyl Nodes Stabilized by Screw Rotation Symmetry in Space Group 19**

**Wenlong Gao 1,\* and Yao-Ting Wang 2,3,\***


Received: 20 April 2020; Accepted: 6 July 2020; Published: 12 July 2020

**Abstract:** Topological photonics have developed in recent years since the seminal discoveries of topological insulators in condensed matter physics for electrons. Among the numerous studies, photonic Weyl nodes have been studied very recently due to their intriguing surface Fermi arcs, Chiral zero modes and scattering properties. In this article, we propose a new design of an ideal photonic Weyl node metacrystal, meaning no excessive states are present at the Weyl nodes' frequency. The Weyl node is stabilized by the screw rotation symmetry of space group 19. Group theory analysis is utilized to reveal how the Weyl nodes are spawned from line nodes in a higher symmetry metacrystal of space group 61. The minimum four Weyl nodes' complex for time reversal invariant systems is found, which is a realistic photonic Weyl node metacrystal design compatible with standard printed circuit board techniques and is a complement to the few existing ideal photonic Weyl node designs and could be further utilized in studies of Weyl physics, for instance, Chiral zero modes and scatterings.

**Keywords:** Weyl nodes; screw rotation symmetry; line node; space group 19; space group 61

#### **1. Introduction**

Weyl nodes (WN) are linear band crossings existing in odd dimensions [1–3]. Three dimensional WN have been extensively exploited in recent years for their intriguing physics. As a drain/source of Berry curvatures, WN hold a quantized topological index, Chern number, of ±1. Due to the stability of the topological index, WN are extremely tolerant to disorders and can be gapped only by the coalesce and annihilation between WN with opposite Chern numbers. WN are also renowned for their Chiral Anomaly Magnetoresistance [4–6], Fermi arc surface states [7–15] and nonlinear Hall effects [16–18]. More recently, topological nodes beyond WN [19], for instance, Dirac points [20–26], spin-1, spin-2/3 WN and double WN [27–30] have also gained significant progress in research.

In photonics, WN have been exploited in photonic crystals [31–34], metamaterials [35–37], magnetized plasma [38,39] and in synthetic spaces [40,41]. Despite substantial studies, there have been only a few ideal WN designs. Ideal WN are salient because no excessive modes are present at the WN frequency, which is beneficial for revealing the properties of WN without disturbance coming from excessive bands [42]. In this article, we proposed a printed circuit board (PCB) technique compatible ideal WN metacrystal functioning in the millimeter wavelength regime (~36 GHz). Group theory analysis shows that the WN are stabilized by screw rotation symmetry in space group (SG) 19, and are spawned from ideal line nodes (LN) in SG 61. Fermi arcs connecting WN with opposite Chern numbers are also found.

#### **2. Results**

The metacrystals we propose comprise two double-layer PCB boards in one primitive cell. The primitive cell's dimension is 3 <sup>×</sup> 3 <sup>×</sup> 3 mm3. Each PCB board layer's thickness is 1 mm, hence, thickness of the dielectric layer between the PCB boards is 0.5 mm and has a relative permittivity of 2.2, which is within the common PCB board material Teflon's parameter range. The PCB board's relative permittivity is 2.2 as well. All the metallic layers have the standard 37μm thickness and are assumed to work as perfect electric conductors within the frequency range of interest. The metacrystal is shown in Figure 1a. Each PCB layer consists of double-layer metallic wire-like structures and the layers are electrically connected by metal-coated 0.2 mm radius through holes. The metallic pads connecting the through holes to the wires have radii of 0.4 mm. Note that in the figures, the primitive cell in the x–y direction is denoted by the red dashed lines. For clarity, the two layers are not drawn together, though it should be remembered that in the real structure, the two layers are gapped by a merely 0.5 mm thickness dielectric board. The designed structure belongs to SG 61 (Pbca). It will be shown later in the article that this metacrystal could be reduced to SG 19 (P212121) by introducing a deformation, as shown in Figure 2a, of which the inversion symmetry and glide symmetries are broken, and are later essential for the creation of the WN.

**Figure 1.** (**a**) First and second layer structure of the line node metacrystal. (**b**) Band diagram of the line node structure. (**c**) Detailed band diagram of the line node. (**d**) Brillouin zone and the path for band diagram. (**e**) *Ex* field components of the line nodes eigen modes found at *ky* = <sup>π</sup> <sup>2</sup>*<sup>a</sup>* , *kx* = 0, *kz* = 0 and the frequency ω at 34.5 and 36.9 GHz, respectively.

**Figure 2.** (**a**) First and second layer of the perturbed metacrystal structure. (**b**) Band diagram of the perturbed metacrystal. (**c**) Detailed band diagram of the metacrystal around the Weyl node frequency. The Weyl nodes are found on the Γ*X* and Γ*Y* high symmetry lines, and the band crossings at other momentums are gapped.

#### *2.1. Line Node from SG 61*

Without the deformation, the metacrystal structure belongs to SG 61, whose group representatives are translational symmetry **T**, three glide symmetries *Gx*,*y*,*z*, three screw rotation symmetries *S*2,*x*,*y*,*z*, and inversion symmetry P. The nonsymmorphic group operations are expressed explicitly as:

$$G\_{\mathbf{x}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \rightarrow \left(-\mathbf{x}, \mathbf{y}, \mathbf{z} + \frac{1}{2}\right)$$

$$G\_{\mathbf{y}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \rightarrow \left(\mathbf{x} + \frac{1}{2}, -\mathbf{y}, \mathbf{z}\right)$$

$$G\_{\mathbf{z}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \rightarrow \left(\mathbf{x}, \mathbf{y} + \frac{1}{2}, -\mathbf{z}\right)$$

$$S\_{2\mathbf{x}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \rightarrow \left(\mathbf{x} + \frac{1}{2}, -\mathbf{y}, -\mathbf{z} + \frac{1}{2}\right)$$

$$S\_{2\mathbf{y}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \rightarrow \left(-\mathbf{x}, \mathbf{y} + \frac{1}{2}, -\mathbf{z}\right)$$

$$S\_{2\mathbf{z}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \rightarrow \left(-\mathbf{x} + \frac{1}{2}, -\mathbf{y}, \mathbf{z} + \frac{1}{2}\right)$$

Nonsymmorphic group operations have been extensively used to protect degeneracies at the boundaries of the Brillouin zone beyond the capabilities of point group [43,44]. A band diagram of the metacrystal is given in Figure 1b, whose sweeping path is given in Figure 1c, and presents band crossings along the high symmetry lines Γ*X S*Γ and Γ*Y* that is line node (LN) degeneracy. Noticeably, the band diagram is pseudo-gapped for the LN, meaning no other excessive bands can be found

around the LN frequency at other sites in the Brillouin zone. Furthermore, on the edge of the Brillouin zone (RS, RU and RT high symmetry lines), Dirac nodal lines are found, featuring four-fold degeneracy on the whole edges [45].

#### *2.2. Group Theory Analysis*

Since the LN is due to an accidental degeneracy that resides on the high symmetry line/plane, the two bands consisting of LN degeneracy are expected to belong to different irreducible group representations. On the high symmetry lines Γ*X* and Γ*Y*, the little group reduced to *G*Γ*<sup>X</sup>* = *E*, *S*2*y*, *Gx*, *Gz*, *T* and *G*Γ*<sup>Y</sup>* = *E*, *S*2*x*, *Gy*, *Gz*, *T*}, respectively. Despite the existence of the half unit cell translations within the nonsymmorphic operations, the coset groups *G*Γ*X*/*T*, *G*Γ*Y*/*T* of the little groups are found to be isomorphic to simple point group *C*2*v*, whose character and compatability tables are given below in Tables 1 and 2. Note that on the Γ*X* and Γ*Y* high symmetry lines, the σ*v* mirror operation are σ*y* and σ*x* while σ *<sup>v</sup>* is σ*z*.

**Table 1.** Character table of group *C*2*v*.


**Table 2.** Compatibility table of group *C*2*v*.


Field distributions of the electrical component *ex* = *Exe*−*ikr* on the momentums half way on Γ*X* and Γ*Y* are given in Figure 1e. Note that after multiplying the phase factor *e*−*ikr*, *ex* is the periodic function in the primitive cell. The left two and the right two images belong to the same frequency, respectively. In the figures, the mirror plane and the rotation axis are illustrated by the solid black lines, and the half unit cell translations related to the nonsymmorphic group operations are illustrated by the yellow arrow. In the left two figures in Figure 1e, it can be inspected that the state is an even state under the glide operation, and an odd state under the screw rotation operation from their (0 1 0) and (1 0 0) direction's fields. The group representation is, thus, *B*<sup>1</sup> according to Table 1. Whereas, for the other state consisting the LN, it is odd under the glide operation and even under the screw rotation operation, meaning the group representation is *A*1.

#### *2.3. Weyl Point from SG 19*

What the group theory analysis can show is that after introducing certain structure deformation, whether the band crossing is kept or gapped can be predicted. Introducing the deformation shown in Figure 2a will essentially reduce the original SG 61 to SG 19, since the only remaining symmetries are the three screw rotations. Therefore, on the Γ*X* and Γ*Y* high symmetry lines, the little groups are isomorphic to point group *C***<sup>2</sup>** (Table 2)**.** In the compatibility relation in Table 2, A1 and B1 can be reduced to different irreducible representations *A* and *B*, meaning the band crossings are kept, while gapped when k is away from the high symmetry lines.

Instead, if the original metacrystal is reduced to SG 29 by keeping *Sx*,*<sup>y</sup>* and *Gz*, the isomorphic point group is *C <sup>v</sup>* in the compatibility relation in Table 2. The irreducible representations are *A* for both the bands, meaning the LN is immediately gapped. This is consistent with the fact that Weyl nodes cannot be found on mirror planes since mirror operations can flip the chirality of the Weyl nodes. The band structure of the SG 29 metacrystal is given in Figure 3, showing a full band gap around 36 GHz.

**Figure 3.** Band structure of the SG 29 metacrystal.

#### *2.4. Surface States of the LN and the WN*

It has been well understood that LN and WN hold surface states. LN hold the 'drumhead' surface states and are normally flat bands that could enhance interactions and are promising for high temperature superconductivity. The WN, on the other hand, hold the so-called Fermi arcs that have been exploited in their intriguing transport properties and unconventional quantum Hall effects [14].

To explore the surface states of the metacrystal, we created super cell configurations that are periodic in the x and y direction and are confined with perfect electrical conductors (PEC) with 10 unit cells in the z-direction. The results are illustrated in Figure 4 and the configuration of the projected surface Brillouin is give in Figure 5. For the unperturbed LN metacrystal, surface states that are reminiscences of the drumhead surface states are found, while for the WN metacrystal, Fermi arcs are found around the momentums where LN are found (red dots in Figure 4b).

**Figure 4.** (**a**) Surface states (gray) and projected bulk band structure of the line node metacrystal in SG 61. (**b**) Surface states (gray/red) and projected bulk band structure of the line node metacrystal in SG 61; the topological Fermi arcs connecting the Weyl nodes are in red.

**Figure 5.** (**a**) Equi-frequency contour of the projected Weyl nodes and the Fermi arcs. Red and Magenta lines are on the top and bottom surface, respectively. (**b**,**c**) Excitation of topological fermi arc surface states with and without a square-shaped defect. Source is an x direction line current fixing the *kx* at 0.2π/*a*, marked by the yellow pentagram. The metacrystals are bounded by the perfect electric conductor (PEC) boundary condition. The unit cell is outlined by the solid black line box.

At a fixed frequency (36.5 GHz in Figure 5), Fermi arcs are found to connect WN with opposite Chern numbers (WN1 and WN2 located on the *kx* and *ky* axis, respectively), which are the results of the integral of Berry curvatures on a closed surface that include a WN [1]. The Fermi arcs are illustrated by the magenta and red lines that are residing on the top and the bottom surfaces, respectively. Figure 5b shows the real space field distribution of the topological surface state excited by a line current source that fixes the *kx* at 0.2π/*a*. The unidirectionality of the surface state propagation is consistent to the calculated equi-frequency contour in Figure 5a. Figure 5c shows the same surface state excitation simulation with an extra square-shaped defect. Again, the unidirectionality demonstrates the topological robustness of the fermi arc surfaces.

#### **3. Discussion**

Our design has shown how to obtain WN in SG 19 from LN in SG 61. The pure PCB layered design could also benefit the tunings of the WN. In Figure 6, we show the band diagrams of the WN metacrystal with various z-direction periods. Note that the PCB layers' thicknesses are conserved at 1 mm and only the interlayer dielectric boards' thickness is changed. Intriguingly, locations of the WN in the Brillouin zone can be effectively tuned by the periodicity. Note that when the thickness is reduced to below 4.8 mm, the WN coalesce and annihilate. The highly sensitive locations of the WN to the periodicity could be used to generate giant effective magnetic fields in the WN metacrystal by introducing a gradually changing thickness. Since photons do not respond to actual magnetic fields, effective magnetic fields are crucial for mimicking the physical effects of real magnetic fields, for example, Landau quantization and chiral zero modes [46,47].

**Figure 6.** Band diagrams of the Weyl node for various unit cell z-direction periods. (**a**) The Weyl nodes coalesce and annihilate at 4.6 mm. (**b**) Weyl nodes are in the proximity of Γ point in the Brillouin zone at 4.8 mm. (**c**,**d**) Weyl node at 5 and 5.5 mm period.

**Author Contributions:** Conceptualization, W.G. and Y.-T.W.; software, W.G.; validation, W.G. and Y.-T.W.; formal analysis, W.G. and Y.-T.W.; writing—original draft preparation, W.G.; writing—review and editing, W.G. and Y.-T.W.; supervision, W.G. and Y.-T.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 724306); Engineering and Physical Science Research Council (EP/T002654/1); Leverhulme Trust (Topologically protected flexural waves in thin elastic plates).

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

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Bulk Cyclotron Resonance in the Topological Insulator Bi2Te3**

**Dmytro L. Kamenskyi 1,2,\*, Artem V. Pronin 3,\*, Hadj M. Benia 4, Victor P. Martovitskii 5, Kirill S. Pervakov <sup>5</sup> and Yurii G. Selivanov <sup>5</sup>**


Received: 28 July 2020; Accepted: 16 August 2020; Published: 20 August 2020

**Abstract:** We investigated magneto-optical response of undoped Bi2Te3 films in the terahertz frequency range (0.3–5.1 THz, 10–170 cm<sup>−</sup>1) in magnetic fields up to 10 T. The optical transmission, measured in the Faraday geometry, is dominated by a broad Lorentzian-shaped mode, whose central frequency linearly increases with applied field. In zero field, the Lorentzian is centered at zero frequency, representing hence the free-carrier Drude response. We interpret the mode as a cyclotron resonance (CR) of free carriers in Bi2Te3. Because the mode's frequency position follows a linear magnetic-field dependence and because undoped Bi2Te3 is known to possess appreciable number of bulk carriers, we associate the mode with a bulk CR. In addition, the cyclotron mass obtained from our measurements fits well the literature data on the bulk effective mass in Bi2Te3. Interestingly, the width of the CR mode demonstrates a behavior non-monotonous in field. We propose that the CR width is defined by two competing factors: impurity scattering, which rate decreases in increasing field, and electron-phonon scattering, which exhibits the opposite behavior.

**Keywords:** topological insulators; cyclotron resonance; Dirac materials

#### **1. Introduction**

From the theory point of view, a three-dimensional (3D) topological insulator (TI) possesses insulating bulk and conducting surfaces, the conduction channels at surfaces being spin-polarized [1–4]. Since the spin polarization can potentially be utilized in spintronic devices, topological insulators have attracted a lot of attention in the past years [5,6]. In practice, the real samples of 3D topological insulators often conduct not only on their surfaces, but also in the bulk. Considerable efforts have been made to understand and separate the properties of surface and bulk charge carriers. These properties can particularly be studied via different spectroscopic techniques, such as angle-resolved photoemission spectroscopy (ARPES) or optical and magneto-optical spectroscopy. The optical conductivity and cyclotron resonance (CR) of a number of 3D TI materials have been reported in the literature. Perhaps the most studied family of such TIs is the bismuth selenide: bismuth telluride series, Bi2(Te1 <sup>−</sup> *<sup>x</sup>*Se*x*)3, which also includes the undoped members, Bi2Te3 and Bi2Se3 [7–9]. In this study, we concentrate on Bi2Te3. Namely, we investigate experimentally the CR in this compound. Surprisingly, the published cyclotron-resonance measurements performed on this well-studied TI produce rather diverging

results [10–13] with the absorption features being generally of rather complex shapes. One of the reasons for such diversity might be the sample-dependent variation between the surface and bulk contributions, which, in turn, greatly depend on the exact position of the Fermi level.

Unlike in the majority of previous reports [11–13], the CR absorption observed in our study can be well described by a single Lorentzian-shaped mode (which is rather consistent with the earliest study on this issue from 1999 [10]). We believe the absorption we detect is of bulk origin. Thus, our findings might be useful for the proper interpretations of the CR modes in doped Bi2Te3, where the balance between the surface and bulk-states contributions can be shifted towards the former, but the bulk still cannot be completely ignored.

#### **2. Materials and Methods**

We grew thin layers of Bi2Te3 on (111)-oriented BaF2 substrates by molecular beam epitaxy [14]. For the growth, we used binary Bi2Te3 and elemental Te. This is different from the standard practice, when elemental (Bi and Te) sources are utilized with the typical flux ratio of Te/Bi being about 10 to 20. Using Bi2Te3 and Te allowed us to reduce this ratio to the values below 1 and, hence, to precisely control the stoichiometry of the growing layer. Along with the employed "ramp up" growth procedure [15], these two approaches successfully suppress twin formation in the growing films. X-ray diffraction (XRD) φ scans about the [0 0 1] axis on the asymmetric (1 0 10) reflection revealed only 120◦-periodic peaks and confirmed that the films obtained by this method are either single-domain or have a very small twin volume fraction (3–7% for the films with 1 cm2 area) with the *c* axis of Bi2Te3 being perpendicular to the substrate surface. To the best of our knowledge, this thin-film growth method is unique.

In order to prevent possible influence of atmospheric oxygen and water, we have developed a method to cover the TI films in situ with optically-friendly protecting layers of BaF2 [16]. We have found that 30–50 nm of BaF2 provide the optimal protection. Our measurements have shown that the BaF2 cap layers affect neither crystal-structure parameters nor optical properties at the frequencies of interest.

The sample used in this study was thoroughly characterized by XRD, scanning electron microscopy (SEM), atomic force microscopy (AFM), and angle-resolved photoemission spectroscopy (ARPES). The results of these investigations, presented in the Supplementary Materials, confirm high structural and morphological quality of the film and show that the film possesses the topological surface electronic states as well as the states in the bulk conduction band.

For optical measurements, we utilized the infrared optical setup available at the High Field Magnet Laboratory in Nijmegen [17]. This setup consists of a commercial Fourier-transform infrared (FTIR) spectrometer (Bruker IFS113v) (Bruker Optik GmbH, Ettlingen, Germany) combined with a continuous-field 33-Tesla Bitter magnet. A detailed description of this setup could be found elsewhere [18]. The measurements were performed in the Faraday geometry [19] at 2 K. A mercury lamp was used as a radiation source. The far-infrared radiation was detected using a custom-made silicon bolometer operating at 1.4 K. The FTIR spectra were recorded in a number of magnetic fields from 0 to 30 T. The optical data were collected between 10 and 170 cm−<sup>1</sup> (300–5100 GHz), using a 200-μm Mylar beamsplitter and a scanning velocity of 50 kHz. At each field, at least 100 scans were averaged. As will be seen below, the data obtained in the fields above 10 T cannot be used in our analysis because of a low signal-to-noise ratio. Thus, in this study we concentrate on the measurements performed between 0 and 10 T.

#### **3. Results and Discussion**

In Figure 1, we show raw transmission data measured though the 115-nm-thick Bi2Te3 film on a 0.49-mm-thick BaF2 substrate (cf. S2 in the Supplemental Material) in magnetic fields *B* up to 10 T. We note that all the measurements reported in this study are performed on a single sample. As seen from Figure 2, the substrate has no detectable field dependence. Hence, all the field-inducted changes come from the film. We note that the BaF2 substrate has intense phonon modes at roughly 50 and

140 cm−<sup>1</sup> [20]. Thus, accurate measurements around these frequencies are impossible. The spectra of Figure 1 are dominated by a single broad mode, which position shifts to higher frequencies in increasing field. The spectra can be fitted by a single Lorentzian, as exemplified in Figure 3 for 2, 4, and 6 T. In zero field, the Lorentzian central frequency is zero, i.e., the observed absorption mode is due to free carriers (Drude conductivity). The field evolution of the mode can be traced in Figure 1: with increasing field, the mode shifts upwards and eventually goes above 100 cm<sup>−</sup>1, i.e., in the range, where the signal-to-noise ratio is worsened by the spectrometer noise and the phonons in the substrate (this prevents a meaningful spectra analysis in higher fields). Still, the shift of the mode in the applied magnetic field is apparent and can straightforwardly be interpreted as a magnetic-field-induced free-carrier localization or, in other words, a cyclotron-resonance absorption.

**Figure 1.** Raw transmission spectra of Bi2Te3 films on BaF2 substrates as obtained in magnetic fields of up to 10 T. The areas with low signal due to either the substrate phonons or spectrometer electronic noise are shaded. The signal-to-noise ratio is best at around 80 cm−<sup>1</sup> and becomes appreciably lower as frequency increases (see also Figure 2), preventing thus any meaningful measurements of the cyclotron resonance (CR) mode at the fields higher than 10 T.

**Figure 2.** Frequency-dependent transmission of a bare BaF2 substrate at 10 and 15 T normalized to its zero-field spectrum. The noise at high frequencies is due the experimental setup. The Figure is meant to demonstrate: (i) the absence of any field-induced changes in the optical spectra of BaF2 and (ii) the frequency limits of the setup used.

The Lorentzian-fit results for this CR absorption mode in the fields from 0 to 10 T are shown in in Figure 4. One can see that the central frequency of the absorption line is linear in field (left panel). This immediately signals that the electronic band(s) responsible for the observed absorption have a quadratic dispersion relation. For linear electronic bands, the field dependence of the CR lines is supposed to have a square root dependence on applied field [21]. Thus, following the Occam's razor principle, we conclude that the mode is due to bulk (i.e., not linear, not topological) electronic bands. This conclusion is in full agreement, e.g., with ARPES [9] and quantum-oscillations [22] measurements, which show that the Fermi level in undoped Bi2Te3 crosses the bulk conduction band and hence there exists a large bulk Fermi surface.

**Figure 3.** Examples of Lorentzian fits of the transmission spectra from Figure 1 for a few magnetic-field strengths as indicated. The raw experimental data are smoothed, using a Savitzki–Golay method [23]. Note that the spectra for 4 and 6 T are shifted upwards for clarity.

**Figure 4.** CR-line central frequency (left frame) and full width at half maximum (FWHM) (right frame) versus applied magnetic field. The dashed line is a fit with *m*\* = 0.1*m*e. The red solid line is a guide for the eye.

We note that weak modes due to the surface conduction channels may exist on top of the dominating bulk abortion, but within our accuracy they cannot be resolved.

The linear field dependence of the central CR frequency, ω0, can be fitted with the standard parabolic-band expression, connecting the slope of ω0(*B*) and the carrier cyclotron mass, *m*\*, ω<sup>0</sup> = *eB*/*m*\**c* (CGS units are used, *e* is the elementary charge, *c* is the speed of light). This fit is shown in the left panel of Figure 4 as a straight line and provides *m*\* = 0.1*m*<sup>e</sup> (*m*<sup>e</sup> is the free-electron mass). This value is in very good agreement with the available literature data on the bulk effective mass in Bi2Te3 [24]: *m*\* = 0.109*m*<sup>e</sup> for the response perpendicular to the *c*-axis, which we do probe in our transmission experiment with unpolarized light. This match provides another confirmation for the correctness of our interpretation. We would like to note here that the complete agreement between the calculated electronic band structure of Bi2Te3 and the entire body of the available experimental work is still to be achieved, as emphasized in a recent review [25].

Finally, we turn to the width of the absorption band. As one can see from the right panel of Figure 4, the full width at half maximum (FWHM) of the band demonstrates a non-monotonous field dependence: in low fields, it decreases with increasing *B* and then, starting at approximately 5 T, the FWHM starts growing with the applied field. The initial decrease of FWHM can be naturally explained by the decreasing cyclotron-orbit radius with increasing *B* and the consequent decrease of impurity scattering. The reason for the CR mode broadening in *B* > 5 T is not entirely clear. We propose that this could be due to the increased electron-phonon scattering. In higher fields, the CR mode approaches the frequencies, where phonon density grows (roughly, above 40 cm<sup>−</sup>1; cf. the left panel of Figure 4 and [26], where the phonon density for Bi2Te3 was calculated) and hence the rate of the electron-phonon scattering starts to increase, leading to the observed total broadening of the CR line according to the Matthiessen rule.

#### **4. Conclusions**

We have investigated the magneto-optical response of undoped Bi2Te3 films at terahertz frequencies and in magnetic fields of up to 10 T. We observed an intense CR line, which can be fitted with a single Lorentz oscillator. The central frequency of the CR increases linearly with applied field, signaling the bulk origin of this resonance. In addition, we found the "in-plane" cyclotron mass, *m*\* = 0.1*m*e, which matches well the literature data for bulk Bi2Te3. The width of the CR mode demonstrates a behavior non-monotonous in field. We propose that the CR width is defined by two competing factors: impurity scattering, which rate decreases in increasing field, and electron-phonon scattering, which rate demonstrates the opposite behavior. We believe our findings can be exploited in future measurements of the surface-states CR in Bi2Te3 to disentangle the bulk and surface contributions.

#### **Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4352/10/9/722/s1.

**Author Contributions:** Far-infrared optical measurements and data analysis, D.L.K. and A.V.P.; ARPES measurements, H.M.B.; sample preparation and characterization, V.P.M., K.S.P., and Y.G.S.; writing—original draft preparation, D.L.K.; writing—review and editing, A.V.P. and Y.G.S.; funding acquisition, D.L.K., A.V.P., and Y.G.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported by the PRIME program of the German Academic Exchange Service (DAAD) with funds from the German Federal Ministry of Education and Research (BMBF) and by the Russian Foundation for Basic Research (grant No. 20-02-00989). H.M.B. acknowledges funding from the German Research Foundation (DFG) via Project No. BE 5190/1-1. We also acknowledge the support of HFML-RU/NWO, member of the European Magnetic Field Laboratory (EMFL).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Crystal Growth by the Floating Zone Method of Ce-Substituted Crystals of the Topological Kondo Insulator SmB6**

#### **Monica Ciomaga Hatnean \*, Talha Ahmad, Marc Walker, Martin R. Lees and Geetha Balakrishnan**

Department of Physics, University of Warwick, Coventry CV4 7AL, UK; S.T.S.Ahmad@warwick.ac.uk (T.A.); M.Walker@warwick.ac.uk (M.W.); m.r.lees@warwick.ac.uk (M.R.L.); G.Balakrishnan@warwick.ac.uk (G.B.) **\*** Correspondence: M.Ciomaga-Hatnean@warwick.ac.uk

Received: 19 August 2020; Accepted: 14 September 2020 ; Published: 17 September 2020

**Abstract:** SmB6 is a mixed valence topological Kondo insulator. To investigate the effect of substituting Sm with magnetic Ce ions on the physical properties of samarium hexaboride, Ce-substituted SmB6 crystals were grown by the floating zone method for the first time as large, good quality single crystal boules. The crystal growth conditions are reported. Structural, magnetic and transport properties of single crystals of Sm1−*x*Ce*x*B6 (*x* = 0.05, 0.10 and 0.20) were investigated using X-ray diffraction techniques, electrical resistivity and magnetisation measurements. Phase composition analysis of the powder X-ray diffraction data collected on the as-grown boules revealed that the main phase was that of the parent compound, SmB6. Substitution of Sm ions with magnetic Ce ions does not lead to long-range magnetic ordering in the Sm1−*x*Ce*x*B6 crystals. The substitution with 5% Ce and above suppresses the cross-over from bulk conductivity at high temperatures to surface-only conductivity at low temperatures.

**Keywords:** crystal growth; optical floating zone method; SmB6; Sm1−*x*Ce*x*B6; topological insulator; kondo insulator

#### **1. Introduction**

Extensive investigations of the physical properties and excitations in the rare earth (RE) hexaboride compounds, REB6, have been carried out over the past decades. These strongly correlated electron systems display an array of interesting magnetic and electrical properties, such as superconductivity (YB6 [1–3]), intricate antiferromagnetic ordered phases owing to the displacement of rare earth ions within the rigid framework formed by the boron ions (GdB6 [4–6]), complex antiferromagnetic phases with Kondo-like characteristics (CeB6 [7–10]), semimetallic behaviour correlated with the transition to an unusual ferromagnetic state (EuB6 [11–13]), typical metallic behaviour (LaB6 [14–16]) or an exotic Kondo-like topological insulating state (SmB6 [17–19]). Amongst the rare earth hexaboride compounds, cerium and samarium hexaborides have puzzled experimentalists and theoreticians alike, for a long time, in view of their intriguing physical properties. SmB6 and CeB6 are isostructural, crystallising in the same cubic CsCl-type structure (*Pm*3*m* space group) [20–22]. Sm and Ce ions replace the Cs ion, whilst the B6 cubo-octohedral clusters take the place of the Cl ions at the corners of the cube. Nevertheless, the similarities between samarium and cerium hexaborides stop at the structural level, as they display very unusual, but different physical properties.

SmB6 has long been known to be a Kondo insulator [23,24]; in recent times, new theoretical and experimental studies demonstrated that samarium hexaboride is a topological Kondo insulator (TKI) exhibiting topological surface properties [18,25–31], although this remains open to further investigation [32]. SmB6 is one of the most investigated Kondo insulators, mainly due to its exciting low temperature transport behaviour. As the temperatures decreases, an energy gap arises due to the interaction of the strongly correlated *f*-electrons and the conducting *d*-electrons, leading to an exponential increase in the electrical resistance of SmB6 [33–35]. Unexpectedly, upon further cooling, the resistance of SmB6 does not continue rising, as would be the case for a conventional insulator, but instead the resistance saturates at a finite value, below 5 K. This plateau in the resistivity has been attributed to a transition from a bulk conductivity characteristic of the high temperature region to a surface-dominated conductivity with bulk insulation at low temperature [36]. SmB6 is a mixed valence system that does not order magnetically, despite exhibiting antiferromagnetic correlations [33,37–40]. The Sm3<sup>+</sup>:Sm2<sup>+</sup> ratio was determined to be independent of the temperature, and equal to approximately 0.6∼0.7:0.4∼0.3 [37,41]. Nevertheless, recent studies have shown that, upon the application of an external pressure, the Sm3<sup>+</sup> configuration can be stabilised for sufficient time to allow long-range magnetic ordering of the samarium ions [42,43].

CeB6 is known to have a typical dense Kondo compound behaviour and a complex magnetic phase diagram [44–47]. Cerium hexaboride exhibits Kondo-like behaviour and has a Kondo temperature of *T*<sup>K</sup> = 19 K. Upon cooling, CeB6 undergoes two magnetic ordering transitions: the first to a state in which antiferroquadrupolar and field-induced octupolar order coexist, below *T*<sup>Q</sup> = 3.2 K, and then to an antiferromagnetic ordering of the Ce dipoles, below *T*<sup>N</sup> = 2.3 K. Moreover, a subsequent study reported a new transition, of unknown origin, at *T*<sup>2</sup> = 1.6 K [45].

Recent progress, e.g., the discovery of the coexistence of an unusual metallic surface state and an insulating bulk state in SmB6 [19,48] and the observation of the long-range-ordered multipolar phases in CeB6 [47], has generated new interest in these materials. One route towards the investigation of the exotic metallic surface state arising in SmB6 and understanding of its topological nature, is through chemical substitution in this TKI with other rare earth ions. Recently, studies have been carried out on Eu, Gd, La, Y and Yb-substituted SmB6 [18,21,49–54]. High levels of substitution of non-magnetic ions (above 30%), and substitutions with small amounts of magnetic ions, were found to destroy the saturation seen in the low temperature resistivity of pure SmB6. It would therefore be interesting to investigate the effect that the substitution with a magnetic rare earth, such as Ce, in samarium hexaboride has on the robustness of the topological surface state of this TKI. Such an investigation is of course best carried out on high quality single crystals. In the present work, we investigated single crystals of Sm1−*x*Ce*x*B6, with a focus on studying the effects that the substitution of the magnetic Ce ion have on the structural and physical properties of SmB6. The physical properties of Ce-substituted SmB6 samples have previously been investigated; however, this has only been done on polycrystalline samples and flux grown crystals [21,49,54]. Crystals of pure cerium and samarium hexaboride have previously been grown using the floating zone (FZ) technique [55–57]; however, Ce-substituted SmB6 compounds have only been grown in crystal form using the flux method [54]. SmB6 crystals grown using Al flux could suffer from contamination by the flux affecting some of the physical properties of the crystals, thereby making it difficult to study the intrinsic properties of pure samarium hexaboride [58]. We have successfully grown, for the first time, crystal boules of Sm1−*x*Ce*x*B6 by the FZ method, which yields large, good quality crystals, free from flux or crucible contamination. The crystals obtained are especially suitable for the investigation of how the substitution with magnetic ions affects the surface and bulk behaviour of this interesting TKI.

#### **2. Materials and Methods**

Crystal boules of Sm1−*x*Ce*x*B6 (*x* = 0.05, 0.10 and 0.20) were grown by the floating zone technique [57] using a CSI FZ-T-12000-X\_VI-VP four-mirror xenon arc lamp (3 kW) optical image furnace (Crystal Systems Incorporated, Yamanashi, Japan). The crystal quality was checked using a backscattering X-ray Photonic-Science Laue camera system (Photonic-Science, St Leonards-on-Sea, UK). Single crystal samples were aligned for selected experiments, and rectangular prism-shaped samples with [001], [1-10] and [110] directions perpendicular to the faces of the prism were cut from the Sm1−*x*Ce*x*B6 crystal boules.

Phase composition analysis was carried out using a Panalytical X-Pert Pro MPD diffractometer (Malvern Panalytical Ltd, Malvern, UK) with Cu K*α*<sup>1</sup> radiation (*λ*K*α*<sup>1</sup> = 1.5406 Å). The diffraction patterns were collected at room temperature over an angular range of 10 to 110◦ in 2*θ* with a step size in the scattering angle 2*θ* of 0.013◦ and at various scanning times. The analysis of the X-ray patterns was performed using the Fullprof software suite [59].

Chemical composition of the crystal boules was investigated by energy dispersive X-ray spectroscopy (EDX) using a Zeiss SUPRA 55-VP scanning electron microscope (Carl Zeiss GmbH, Jena, Germany). LaB6 was used as a standard for the EDX measurements. X-ray photoelectron spectroscopy (XPS) analysis was also carried out in order to determine the elemental composition and the valence of the Sm ions. The samples were attached to electrically-conductive carbon tape, mounted on to a sample bar and loaded into a Kratos Axis Ultra DLD (Kratos Analytical Ltd, Manchester, UK) spectrometer (base vacuum of ∼<sup>2</sup> × <sup>10</sup>−<sup>10</sup> mbar). The measurements were performed using a monochromated Al K*α* X-ray source, at room temperature and at a take-off angle of 90◦ with respect to the surface parallel. The data were analysed in the CasaXPS package (Casa Software Ltd, Teignmouth, UK), using Shirley backgrounds and mixed Gaussian-Lorentzian (Voigt) line-shapes and asymmetry parameters, where appropriate.

Magnetic susceptibility measurements were performed with a Quantum Design Magnetic Property Measurement System (Quantum Design Incorporated, San Diego, USA) on rectangular-prism-shaped Sm1−*x*Ce*x*B6 samples with an applied field parallel to the [100] (tetragonal), [110] (rhombic) and [111] (trigonal) crystallographic directions. The samples were cooled to 1.8 K in zero field and then the susceptibility as a function of temperature up to 300 K was measured on warming and then cooling with an applied field of *H* = 500 Oe.

Alternating current (*ac*) resistivity measurements were performed using a Quantum Design Physical Property Measurement System on bar shaped samples of the Sm1−*x*Ce*x*B6 single crystals using the standard four-probe technique. Silver wire contacts were attached with silver paint, in a linear configuration, to the surfaces of the samples. The resistivity measurements were made from 2 to 300 K on both cooling and warming in zero applied field with an ac current of 1 mA at a frequency of 113 Hz.

#### **3. Results and Discussion**

#### *3.1. Crystal Growth*

Stoichiometric ratios of high purity SmB6 (99.9%, American Elements UK, Manchester, UK) and CeB6 (99.5%, Cerac Incorporated, Milwaukee, USA) powders were mixed together by ball milling for over 15 h, to prepare 5%, 10% and 20% Ce-substituted SmB6 polycrystalline samples. The resulting materials were then isostatically pressed into rods (typically 5–7 mm diameter and 40–50 mm long) and sintered in an alumina boat, at 1450 ◦C in a flow of argon gas for 12 h. Before the sintering process, the furnace was evacuated to give a vacuum of ∼10−<sup>5</sup> mbar. The resulting polycrystalline feed rods were used for the crystal growth. A binder (polyvinyl alcohol or polyvinyl butyral) was mixed with the powders in some cases to facilitate the formation of the rods.

Crystals of Sm1−*x*Ce*x*B6 (*x* = 0.05, 0.10 and 0.20) were successfully grown by the floating zone method. The growths were carried out in an argon atmosphere at a pressure of ∼3 bar, using a growth rate of 18 mm/h. The feed and the seed rods were counter-rotated at ∼15–25 rpm. Initially, a crystal boule of SmB6 was used as a seed. Once good quality crystals were obtained, Sm1−*x*Ce*x*B6 crystal seeds were used for subsequent growths. A dark grey coloured deposition on the quartz tube surrounding the feed and seed rods was observed for all the growths, indicating the evaporation of boron during the growth process.

The Sm1−*x*Ce*x*B6 boules were typically 4–5 mm in diameter and measured approximately 45–50 mm in length. All the crystals obtained developed facets as they grew and two very strong facets were present on almost the entire lengths of most of the grown crystals. Figures 1a–c show photographs of Sm1−*x*Ce*x*B6 crystals grown in argon atmosphere at a growth speed of 18 mm/h. The quality of the grown boules was investigated by X-ray Laue diffraction, and Laue photographs were taken along the length of the boule, on the faceted sides (see Figure 1). The Laue patterns were identical along the whole length of the crystal boules.

**Figure 1.** Crystal boules of (**a**) Sm0.95Ce0.05B6, (**b**) Sm0.90Ce0.10B6 and (**c**) Sm0.80Ce0.20B6, prepared by the floating-zone method in argon atmosphere at a growth rate of 18 mm/h. X-ray Laue back reflection photographs show the [001] orientation of aligned Sm1−*x*Ce*x*B6 samples used for the physical properties measurements.

#### *3.2. Structural and Composition Analysis*

Structural and phase purity analysis was carried out using powder X-ray diffraction measurements on small pieces of the Sm1−*x*Ce*x*B6 crystals selected from close to the end of each crystal boule. Figures 2a–c show the patterns for *x* = 0.05, 0.10 and 0.20, and profile matching (goodness of fit, GOF = 1.35, 1.51 and 1.92, respectively) to the cubic *Pm*3*m* space group [20] indicates that in each case the main phase is Sm1−*x*Ce*x*B6, with no significant impurity phases present. One peak that does not belong to the *Pm*3*m* cubic structure can be observed at ∼26.6◦ in the powder X-ray profiles of each of the Sm1−*x*Ce*x*B6 crystals grown. The impurity was identified to be a hexagonal (*P*63/*mmc*) SmBO3 phase [60]. Lattice parameters calculated from the profile matching were determined to be 4.1351(2) Å, 4.1384(2) Å and 4.1393(2) Å, respectively, for Sm0.95Ce0.05B6, Sm0.90Ce0.10B6 and Sm0.80Ce0.20B6 (see Table 1). The values are in agreement with those reported in previous studies on Sm1−*x*Ce*x*B6 polycrystalline samples [21].

**Figure 2.** Powder X-ray diffraction patterns of Sm1−*x*Ce*x*B6 with (**a**) *<sup>x</sup>* = 0.05, (**b**) *<sup>x</sup>* = 0.10 and (**c**) *x* = 0.20) for samples taken from the crystal boules. The experimental profile (red closed circles) and a full profile matching refinement (black solid line) made using the *Pm*3*m* cubic structure are shown, with the difference given by the blue solid line. The orange coloured symbols \* indicate the impurity peaks belonging to SmBO3 impurity phases. (**d**) Evolution of the lattice parameter, *a*, as a function of the concentration, *x*, of the Ce-substituent for Sm1−*x*Ce*x*B6. The experimental values obtained in the present work (red open circles) are also given in Table 1. The previously reported values (red, black and orange closed circles) of the crystallographic parameters for the Sm1−*x*Ce*x*B6−*<sup>y</sup>* series [21] are given for completeness.

Figure 2d shows the dependence of the lattice constant on the concentration of Ce for the Sm1−*x*Ce*x*B6 samples. The composition dependence of the cubic parameter, *a*, does not obey Vegard's law [61], for the Sm1−*x*Ce*x*B6 series. The anomalously large positive deviation observed in Figure 2d can be attributed to the mixed valence of samarium ions [37,38,62]. As the concentration of the Ce-substituent changes from *x* = 0 to 1, the Ce3<sup>+</sup> ions replace the Sm3<sup>+</sup> ions preferentially, whereas the concentration of Sm2<sup>+</sup> ions remains constant [21,38,50]. The effective ionic radius [63,64] of Ce3<sup>+</sup> (1.01 Å) is larger than the ionic radius of Sm3<sup>+</sup> (0.958 Å); thus, the substitution of samarium with cerium ions results initially in a lattice expansion (up to *x* ∼ 0.6). Further substitution of samarium with cerium is followed by a subtle lattice contraction, which is attributed to the replacement of the larger Sm2<sup>+</sup> ions (1.15 Å) with Ce3+. A similar effect on the lattice constant has been observed in the case of Gd and La-substituted SmB6 [38,50,65].


**Table 1.** Lattice parameters calculated from profile matching the powder X-ray diffraction patterns of the Sm1−*x*Ce*x*B6 (*<sup>x</sup>* = 0.05, 0.10 and 0.20) crystals to the *Pm*3*<sup>m</sup>* cubic structure. The previously reported structural parameters quoted for other members of the Sm1−*x*Ce*x*B6−*<sup>y</sup>* series [21] are included for completeness.

Composition analysis of the crystals of Sm1−*x*Ce*x*B6 was carried out by EDX to determine the concentrations of Ce in each crystal. The results, given in Table 2, show that the ratios for Sm:Ce are similar to the expected chemical compositions for the crystals, relative to the starting compositions of the polycrystalline materials (5%, 10% and 20% Ce-substituted SmB6 samples).

**Table 2.** Chemical composition and valence of the Sm ions determined by EDX and XPS for the Sm1−*x*Ce*x*B6 crystal boules grown. The data collected on a pure SmB6 crystal are included for completeness. The XPS measurements were carried out on a piece of an as-grown SmB6 crystal boule and on a sample cleaved (in-situ) from the as-grown SmB6 crystal fragment.


Core level XPS spectra were recorded using a pass energy of 20 eV (resolution ∼0.4 eV) on an area of 300 μm × 700 μm of the Sm1−*x*Ce*x*B6 crystals and used to study the electronic states of Sm 4*d*, Ce 3*d*3/2 and Ce 3*d*5/2 levels, shown in Figure 3. The Sm 4*d* XPS spectrum (see Figure 3a) is composed of one singlet, at 123.5 eV, and one multiplet, at 134.1 eV, separated by approximately 10.6 eV. The Sm2<sup>+</sup> (4 *f* <sup>6</sup> ground-state) feature appears near 129 eV (Sm 4*d* photoelectron line position), which is in agreement with previously published XPS studies on pure SmB6 [66,67]. The Sm3<sup>+</sup> (4 *f* 5) multiplet appears at a higher binding energy, well separated from the 2+ peak. The contributions of the two features to the XPS spectra were used to determine the valence of the Sm ions. The results, given in Table 2, reveal that Sm1−*x*Ce*x*B6 are mixed valence systems, similar to the parent compound SmB6 [37,38]. The average Sm valence values of the Sm1−*x*Ce*x*B6 boules are slightly larger than the values determined previously for pure SmB6 [67–70], due to surface oxidation effects (an increased

concentration of Sm3<sup>+</sup> to the detriment of the Sm2<sup>+</sup> ions). Previous XPS results reported an increased average Sm valence and a B/Sm ratio lower than the nominal stoichiometric value of 6:1 when the SmB6 crystals were exposed to ambient conditions [67]. To confirm this hypothesis, XPS spectra were collected on two SmB6 crystal samples, an as-grown and a cleaved crystal fragment. The average samarium valence is ∼2.8 for the as-grown crystal fragment of SmB6. In the case of SmB6 cleaved in-situ from the as-grown crystal, the Sm valence is 2.7, corresponding to a Sm3+:Sm2<sup>+</sup> ratio of approximately ∼0.7:∼0.3, which is in agreement with previous results [37,41].

**Figure 3.** (**a**) Sm 4*d* XPS spectrum and (**b**) Ce 3*d*3/2,5/2 XPS spectra collected for the Sm1−*x*Ce*x*B6 (*x* = 0.05, 0.10, and 0.20) crystal boules.

The Ce 3*d* spectrum, shown in Figure 3b, is comprised of two multiplets, at 885.8 eV and 904.8 eV, corresponding to the spin-orbit split 3*d*5/2 and 3*d*3/2 core levels. The spin-orbit splitting is approximately 19 eV, with the complex electronic structure of different Ce oxidation states yielding useful spectral features which can be used to distinguish Ce3<sup>+</sup> and Ce4<sup>+</sup>. In our data, each component of the Ce 3d XPS spectrum is dominated by two features. The absence of a third component at 916 eV, characteristic of the Ce4<sup>+</sup> (4 *<sup>f</sup>* 0) [71,72], indicates that the Ce ion is in the 3+ state in the Sm1−*x*Ce*x*B6 samples. The analysis of the XPS results, given in Table 2, allowed us to estimate the amount of Ce-substituent in the Sm1−*x*Ce*x*B6 crystal boules. A comparison of the Ce concentrations determined from the XPS spectra and those estimated from the EDX compositional analysis is provided in Table 2.

#### *3.3. Magnetisation*

Zero-field-cooled warming (ZFCW) and field-cooled cooling (FCC) magnetisation versus temperature curves were collected on pieces of the Sm1−*x*Ce*x*B6 (*x* = 0.05, 0.10 and 0.20) single crystals with an applied field of 500 Oe along three different crystallographic directions ([001], [110] and [111]). The temperature dependence of the *dc* magnetic susceptibility, *χ* (*T*), is shown in Figure 4a. The magnetic susceptibility measured along the different crystallographic directions for all three Sm1−*x*Ce*x*B6 crystals decreased on warming from 1.8 K to room temperature. In addition, for each Sm1−*x*Ce*x*B6 composition, the magnetic susceptibilities collected with field applied along the three different high-symmetry directions all overlap to within experimental error across the whole temperature range.

The temperature dependent magnetic susceptibility of the Sm1−*x*Ce*x*B6 crystals was compared with data collected on a pure SmB6 crystal grown by the floating zone method [57]. In the temperature range 300 to 60 K, the magnetic susceptibility of Ce-substituted and pure SmB6 crystals show a similar behaviour, i.e., a gradual increase of *χ* (*T*) with decreasing temperature. Below 60 K, the Sm1−*x*Ce*x*B6 crystals exhibit a more rapid increase in susceptibility, down to 1.8 K. In contrast, the susceptibility data of pure SmB6 crystals contain a broad maximum centred around 50–60 K, characteristic of a

Kondo insulator, before a more gradual upturn at lower temperatures. Moreover, in the temperature range 1.8–60 K, the magnetic susceptibility of Sm1−*x*Ce*x*B6 crystals increases sharply with increasing Ce content. The change in the magnetic response of both Ce-substituted and pure SmB6 crystals below 60 K coincides with the increase observed in the resistivity (see Figure 5).

**Figure 4.** (**a**) Temperature dependence of the *dc* magnetic susceptibility, *χ* versus *T*, in the temperature range 1.8–100 K for the Sm1−*x*Ce*x*B6 (*<sup>x</sup>* = 0, 0.05, 0.10 and 0.20) crystals, with a magnetic field applied along the [001] (black), [110] (red) and [111] (orange) crystallographic directions. The previously reported susceptibility data for a SmB6 crystal [57] are given for comparison. The inset shows *χ* versus *T*, on a logarithmic scale, in the temperature range 1.8–300 K. (**b**) Temperature dependence of the reciprocal of the *dc* susceptibility, *<sup>χ</sup>*−<sup>1</sup> versus *<sup>T</sup>*, of Sm1−*x*Ce*x*B6 for a field applied along the [001] direction. The inset shows the normalised magnetic susceptibilities of Sm1−*x*Ce*x*B6 samples, with a magnetic field applied along the [001] direction. The *χ*/*χ* (10 K) versus *T* data increase rapidly at low temperatures, but with no signature of long-range magnetic order, for all Ce concentrations.

**Figure 5.** Temperature dependence of the bulk *ac* resistivity, *ρ* versus *T*, in the temperature range 1.8–300 K for the Sm1−*x*Ce*x*B6 (*x* = 0, 0.05, 0.10 and 0.20) crystals. The previously reported resistivity data for a SmB6 crystal [57] are given for comparison.

Attempts to fit the temperature-dependent reciprocal magnetic susceptibilities, *χ*−<sup>1</sup> (*T*) (see Figure 4b), in the temperature range 100–300 K reveal that pure SmB6 and the Sm1−*x*Ce*x*B6 materials all appear to follow a Curie–Weiss law. The effective moment, *μ*eff, per formula unit at 300 K varies from 2.4(1)*μ*<sup>B</sup> for *x* = 0.00 to 2.5(1)*μ*<sup>B</sup> for *x* = 0.20. The form of *χ* (*T*) for the Sm1−*x*Ce*x*B6 crystals is qualitatively similar than data reported for aluminium flux grown Ce-substituted SmB6 single crystals, although the effective moments in our samples are substantially lower, especially for lower Ce concentrations [54]. The *χ* (*T*) data are consistent with magnetic response expected for a mixture of 4 *f* <sup>1</sup> Ce3<sup>+</sup> ions - 2.54*μ*B/Ce3<sup>+</sup> and divalent and trivalent Sm ions in a variety of magnetic and nonmagnetic electronic configurations, 4 *f* <sup>6</sup> , 4 *f* <sup>5</sup>*d*<sup>1</sup> , and 4 *f* <sup>5</sup> [68,73,74].

A previous study reported that substituting Sm with another magnetic rare earth ion, such as Gd3+, in large concentrations (≥40%), leads to antiferromagnetic ordering at low temperature due to coupling between the Gd sites [18]. This is predicted by the existence of a saturation plateau in the normalised magnetic susceptibility data of 40% Gd-substituted SmB6. In contrast, the magnetisation curves for our Sm1−*x*Ce*x*B6 crystals exhibit a rapid increase at low temperatures, but with no evidence for the onset of long-range magnetic order down to 2 K, as shown in the inset of Figure 4b. For the Ce3<sup>+</sup> concentrations used in our work, the magnetic data suggest that the Ce ions are distributed randomly in the lattice.

#### *3.4. Resistivity*

Alternating current resistivity versus temperature, *ρ* (*T*), measurements were made on bar shaped samples cut from the Sm1−*x*Ce*x*B6 (*x* = 0.05, 0.10 and 0.20) single crystals. The *ac* resistivity data are shown in Figure 5 for temperatures between 1.8 and 300 K. These resistivity data are compared with data for a pure SmB6 crystal grown by the FZ method and reported in our previous work [57]. At 300 K, the Sm1−*x*Ce*x*B6 samples all have resistivity values similar to SmB6 and *ρ* (300 K) increases with *x*. Below 300 K, SmB6 exhibits a continuous increase in the bulk electrical resistivity. In contrast, the *ρ* (*T*) data for the Ce-substituted samples exhibit a broad maximum centred around 150 K, followed by an increasingly prominent minimum at ∼50 K. On further cooling below 50 K, the resistivity of SmB6 increases by four orders of magnitude, whereas the resistivity of the Ce-substituted samples increases by only a single order of magnitude or less. Nevertheless, the Sm1−*x*Ce*x*B6 samples still have resistivities larger than pure CeB6, over the entire temperature range studied. The resistivity of CeB6 is approximately 10−<sup>5</sup> <sup>Ω</sup>-cm from 2 to 300 K [75], whereas for the Sm1−*x*Ce*x*B6 samples it is 10−<sup>3</sup> <sup>Ω</sup>-cm or higher over the same temperature range, for the *x* = 0.20 sample.

In contrast to the saturation plateau seen in the resistivity of SmB6 at lowest temperatures, *ρ* (*T*) for the Sm1−*x*Ce*x*B6 samples increases monotonically with decreasing temperature below 10 K. These results are in agreement with the transport measurements performed on aluminium flux grown Ce-substituted SmB6 single crystals [54]. There is an evolution from the TKI behaviour of pure SmB6 to a dense Kondo system with low temperature spin ordering of CeB6 [7,8,75]. The data suggest that it is bulk conductivity, modified by crystalline electric field and Kondo effects alongside phonon scattering, that determines the form of the *ρ* (*T*) curves for these Sm1−*x*Ce*x*B6 samples over the entire temperature range studied. A more quantitative description of the transport properties of the Sm1−*x*Ce*x*B6 crystals, including extensive measurements in a magnetic field, will be presented elsewhere [76].

#### **4. Conclusions**

Crystal boules of Sm1−*x*Ce*x*B6 (*x* = 0.05, 0.10, and 0.20) compounds were grown, for the first time, by the FZ technique. Investigation of the crystals using X-ray diffraction techniques revealed that the Ce-substituent is successfully incorporated in the SmB6 structure and that the structural distortions due to the substitution of Sm with Ce follow a similar trend to the one reported for polycrystalline samples of Ce-substituted SmB6. EDX and XPS results confirm that the Ce concentration is close to the nominal stoichiometric values of *x* = 0.05, 0.10 and 0.20. Analysis of the average Sm valence data determined by XPS on Sm1−*x*Ce*x*B6 and pure SmB6 samples showed that the results are extremely dependent on the quality of the surface studied; i.e., an increase in the Sm valence is observed when the surface is exposed to ambient conditions. Magnetic property measurements show that our Sm1−*x*Ce*x*B6 crystals exhibit no sign of long-range magnetic ordering, at substitution concentrations below 20%. Temperature dependent resistivity measurements revealed that a 5% (and above) substitution with Ce suppresses the crossover from bulk to surface conductivity seen in pure SmB6 as the temperature is reduced. Detailed low temperature magneto-transport measurements are now being carried out to investigate the bulk and surface properties of the Sm1−*x*Ce*x*B6 crystals.

**Author Contributions:** M.C.H. and T.A. performed the crystal growths; characterisation measurements were carried out by M.C.H. and T.A. with M.W. and M.R.L.; M.C.H., T.A., M.W., M.R.L. and G.B. analysed the data; M.C.H. and T.A. drafted the paper; M.R.L. and G.B. reviewed the manuscript; G.B. secured the funding and managed the project. All authors have read and agreed to the published version of the manuscript.

**Funding:** Financial support was provided by EPSRC, UK, through grant EP/T005963/1.

**Acknowledgments:** The authors thank S. York for the EDX compositional analysis, and A. Julian and T. E. Orton for valuable technical support.

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Optical Response of Chiral Multifold Semimetal PdGa**

**Sascha Polatkan \* and Ece Uykur**

Physikalisches Institut, Universität Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany; ece.uykur@pi1.physik.uni-stuttgart.de

**\*** Correspondence: sascha.polatkan@pi1.uni-stuttgart.de

**Abstract:** We present a theoretical study of the band structure and optical conductivity for the chiral multifold semimetal PdGa. We identify several characteristic features in the optical conductivity and provide their origins within the band structure. As experimental optical studies for the mentioned compound have not been reported, we contrast our results with the related compounds, RhSi and CoSi. We believe that the presented hallmarks will provide guidance to future experimental works.

**Keywords:** topology; chirality; multifold semimetal; optics; DFT

#### **1. Introduction**

First hints of topological characters of materials were found in crystals, where avoided band crossings caused a separation of band energies accompanied by a mixture of band characters. Brought into contact with a different material, that hosts a regular order of bands, without band mixing, this enforces surface states at the interface, which would host linear, intersecting bands with opposite spin characters. Later, it was found that this can be extended to be a bulk property, when the Weyl semimetal was born [1], hosting a three-dimensional bulk realization of linear, spin non-degenerate bands. Their equivalency to the surface states of topological insulators is immediatelly obvious, but can further be quantized by the Chern number *χ*, for which the Weyl semimetal is the lowest integer realization with *χ* = ±1 [2–4]. Weyl semimetals were promising for applications in spintronics [5–7], in optoelectronics [8–11], or even in chemistry [12–16]. Of essence, here, was either the breaking of the time-reversal symmetry, or the inversion symmetry in combination with sufficiently strong spin-orbit coupling (SOC). The first confirmed Weyl-semimetal is TaAs [17–19], which lacks inversion symmetry.

This can be extended to crystals, which host, besides broken inversion symmetry, also a lack of mirror symmetries. Among the 230 space groups, the 65 Sohncke groups fall under this condition and can provide a *chiral* crystal structure. Note that the 65 Sohnke groups are not necessarily chiral space groups, but can provide chiral crystal structures. In total, there are only 22 chiral space groups (11 enantiomorphic pairs). These 22 chiral space groups are contained within the 65 Sohnke groups. The group of PdGa, *P*213 is not a chiral space group, but a Sohnke group [20,21].

It has recently been suggested, that in chiral crystals a new type of fermionic state can be realized, extending the pool of topological quasiparticles by the so called multifold fermion [22]. They differ from the Weyl fermion in that they are hosted by conical band intersection with Chern numbers |*χ*| > 1. To this end, chiral crystal structures belonging to the RhSi family, space group *P*213, number 198, have garnered attention as angle resolved photoemission spectroscopy (ARPES) measurements provided strong indications that these materials host surface states, which are maximally extended in *k*-space [23–28]. The lack of inversion and mirror symmetries, in combination with SOC, leads to the splitting of bands around the high symmetry points Γ and *R*, giving rise to a non-collinear spin arrangement with Chern numbers *χ* = ±4 [29] with maximally extended surface states, ranging from the center to the edge of the Brillouin zone [25,27,30]. In such semimetals, the

**Citation:** Polatkan, S.; Uykur, E. Optical Response of Chiral Multifold Semi-metal PdGa. *Crystals* **2021**, *11*, 80. https://doi.org/10.3390/ cryst11020080

Academic Editor: Yuri Kivshar Received: 16 December 2020 Accepted: 18 January 2021 Published: 21 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional clai-ms in published maps and institutio-nal affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

quantized circular photogalvanic effect has been predicted, which, by effect of topological states, constitutes a photocurrent, quantized in units of material-independent fundamental constants [31–34].

In this paper, we theoretically investigated the material PdGa, which belongs to the same family. So far, PdGa has not been as extensively investigated as its silicon-based counterparts. In contrast to CoSi and RhSi, the Fermi surface of PdGa is considerably larger and is expected to host a larger number of free carriers, perhaps making the experimental observation of the low energy dynamics challenging, as in the case of RhSi [35]. Hence, the theoretical approach would be beneficiary for comparison and provide some guide for the interpretation of the high energy part of the measured spectra.

Furthermore, the optical conductivity provides a useful tool to identify fingerprints of Dirac fermions [36–40], which should express themselves as linear signatures *σ*(*ω*)<sup>1</sup> ∝ *ω* [41,42]. In the silicides, they are thought to occur as low energy excitations near the Γ-point, as the Fermi level almost coincides with the multifold chiral intersection. Even though the Fermi level in PdGa is much higher, we nonetheless find a nearly linear ridge in the optical spectrum and identify its origin, as discussed in the results.

The combination with external stimuli, such as magnetic fields [43–45] or external pressure [46,47], which provide clean, non-invasive tuning mechanisms, makes the optical conductivity especially desirable in the search for topological materials.

We provide theoretical estimates for intra- and interband optical conductivities, which we contrast to and interpret along with the published experimental results of PdGa's sister compounds CoSi [48] and RhSi [35], as no experimental data for PdGa are available.

#### **2. Results and Discussion**

We performed DFT calculations using Wien2k's full-potential linearized augmented plane wave (LAPW) methods with the Perdew–Burke–Ernzerhof (PBE) exchange correlation, accounting for the semimetallic nature of the chiral multifold compound PdGa [49,50]. The lattice parameters were adopted from Ref. [29]. Figure 1 shows the chiral atomic structure and the Fermi surface of PdGa. The Fermi surface is comparatively large, with reference to RhSi and CoSi, indicating a much stronger intraband response. In Figure 1d the Fermi surfaces of the separate contributing bands are shown, sorted by energy from left to right. Within the cube-like structure, a set of eight droplet-shaped Fermi surfaces are hidden, positioned within the corners of the cube.

The band structures of PdGa were calculated on an 18 × 18 × 18 *k*-mesh with and without SOC, as shown in Figure 2a–c. The calculations were converged within 14 cycles down to the charge 10−<sup>5</sup> *e*, with *e* being the electron charge. Core leakage was well within acceptable levels: the core charges for Pd and Ga integrated to 29.998 *e* and 17.999 *e*, respectively. These core electrons arise from choosing a generous energy interval of −9.0 Ry for outer electrons. We chose the parameter *RK*Max = 7.0 and excluded relativistic linear orbitals (RLOs). Due to SOC in combination with the lack of inversion symmetry, bands with opposite spin split throughout the Brillouin zone (BZ) under the action of the Dresselhaus effect. Exempt from the splitting is the Γ-point, at the BZ center, and points or paths on the BZ boundaries. The Γ-point hosts a four-fold degenerate spin 3/2 fermion, also known as the Rarita–Schwinger fermion [51], and the *R*-point two three-fold degenerate spin 1 fermions, each with the maximum possible Chern number of |*χ*| = 4 [29,30]. A detailed view of bandstructures around the Γ-point is given in Figure 2b,c without and with SOC, respectively.

The calculations summarized in Figures 3 and 4 were performed on a denser 32 × 32 × 32 *k*-mesh with SOC, using the optic module [52]. Figure 3 shows the density of states (DOS) of PdGa, as well as the atomic densities of Pd and Ga separately. It can be clearly seen that the majority of carriers at the Fermi surface, as well as electrons involved in optical transitions far beyond the visible spectrum, are entirely contributed by Pd. Significant contributions to the DOS from Ga arise at energies only as low as −15 eV, and beyond.

**Figure 1.** (**a**) The crystal structure of PdGa. The Pd and Ga atoms are arranged chirally along the *c*-axis, with a distinct handedness, giving rise to chiral properties in the band structure and optical interactions. (**b**) Fermi surface of PdGa. (**c**) Brillouin zone of PdGa representing the high symmetry points used in the band structure plots. (**d**) Contributions of the different bands to the Fermi surface, sorted from lowest to highest energy going from left to right.

**Figure 2.** Band structure of PdGa (**a**) with spin-orbit coupling (SOC). Bands that cross the Fermi surface are labeled in pairs *A*<sup>±</sup> to *D*±, with respect to their energy. The Fermi energy is positioned at *EF* = 0 eV. The label ± refer to spin pairs, which split away from high symmetry points. This splitting, due to lack of inversion and mirror symmetries, gives rise to a 4-fold intersection at Γ and 6-fold intersection at *R*, both with Chern numbers of magnitude |*χ*| = 4. (**b**) Magnified view of the band structure around the Γ-point without SOC. (**c**) Magnified view around Γ with SOC.

**Figure 3.** Density of states of PdGa and the atoms Pd and Ga separately. The Gd contribution is negligible, and has been magnified by a factor of 10. Most bands around the Fermi energy are thus contributed by the Pd atoms, underlining the relevance of chirality among carriers and optical transition. The Fermi level is positioned at *EF* = 0 eV as indicated by the blue dotted line.

**Figure 4.** (**a**) The calculated real part of the interband optical conductivity *σ*1(*ω*) with weight analysis. The features were sorted according to the labelling in Figure 2a. (**b**) Intra-, interband, and total optical conductivity *σ*1(*ω*) of PdGa for parameters Γ = 2.5 meV and *ω<sup>p</sup>* = 4.43 eV.

The interband contribution to the optical conductivity of PdGa, without scattering, is shown in Figure 4a. It can be obtained from the imaginary part of the relative permittivity *ε*2, which is the standard output of Wien2k, via the convention (CGS units)

$$
\sigma\_1 = \frac{\omega \varepsilon\_2}{4\pi}.\tag{1}
$$

Several distinct features can be seen. Two features of special interest have been marked <sup>1</sup> , at around 150–300 meV, and <sup>2</sup> , which shows a distinct peak at around 600 meV. The former originates from bands near the Γ-point, see again Figure 2c, corresponding to transitions from the Fermi surface droplets in Figure 1d to the corners of the cube-shaped Fermi surface. The latter is of particular interest as it seems to reproduce a known feature in CoSi around 560 meV [48] very well and is also in the energy range of a broader feature, likely of similar origin, in RhSi, around 750 meV [35], both measured at *T* = 10 K. The sharp peak at 600 meV can be attributed to transitions between the parallel bands dispersing from Γ to *R* and similar transitions between *M* and *R*. Since for CoSi and RhSi the corresponding bands between Γ and *R* are not as parallel, we predict that the optical transition for PdGa will be distinctly sharper. Note that the interband conductivities were not broadened in

our calculations. In the case of RhSi, the mentioned feature merges with a rising edge the in optical conductivity around 1 eV [35], which may indicate what could be identified with the lower edge of feature <sup>3</sup> in our calculations, which arises from a combination of transitions *A* → *C*, *D* and *B* → *D*.

As mentioned earlier, another interesting aspect of the the calculated *σ*1(*ω*), is the nearly linear ridge extending over a rather large energy range, from roughly 0.9 eV to 1.5 eV. Linear features in the optical conductivity are readily identified as signatures of Dirac fermions [36–40]. We can assign it to the transition *B* → *C*. In an experimental setting, this feature may mistakenly be interpreted as transitions between bands extending from the lower lying multifold point, intersecting Γ at around −1.0 eV, to bands above, which we can safely exclude. Note further that the contributing bands contain accidental near-degeneracies, akin to gapped Dirac points, between Γ − *X* and Γ − *R*. Away from the high-symmetry lines, these features may develop into real conical intersections, implying that the discussed transition may very well contain character of Dirac fermions. This warrants closer inspection in future works and may motivate experimental investigations.

In Figure 4b the intra- and interband contributions are plotted alongside the total optical conductivity, *σ*1(*ω*). For the intraband contribution we make the reasonable assumption for the dc conductivity, *<sup>σ</sup>*<sup>0</sup> = <sup>1</sup> · <sup>10</sup><sup>6</sup> <sup>Ω</sup>−1cm−<sup>1</sup> supported by on-site dc measurements, and use the theoretically obtained value for the plasma frequency, *ω<sup>p</sup>* = 4.43 eV (note that this is the unscreened plasma frequency). Both values are significantly larger than the experimentally obtained parameters for CoSi and RhSi [35,48], which is consistent with the much larger Fermi surface of PdGa. These values, under assumption of a single Drude contribution for the intraband transitions, yield a scattering rate of Γ = 2.5 meV, which is expressed by a very sharp Drude peak, seen in Figure 4b. All discussed interband features remain clearly visible after inclusion of the intraband contribution. Whether such a sharp Drude peak is experimentally reproducible remains to be seen. Accounting for interband broadening, the feature <sup>2</sup> might be brought to closer coincidence with the CoSi peak, as well as the much broader RhSi feature. Note that the low energy feature <sup>1</sup> is inherent only to PdGa. This is, again, due to fact that bands, originating from the Γ-point in PdGa, disperse in parallel towards *R*, while in CoSi and RhSi these bands take the form of a double flat band overlaid with a Dirac cone, contributing an energetically much broader joint DOS, smearing out the transition.

#### **3. Conclusions**

In summary, we performed DFT calculations on one of the chiral multifold semimetals, PdGa, and estimated the optical conductivity in a broad frequency range, providing a detailed picture the band structure and a guide for interpretation of future optical experiments on this currently popular compound. A plethora of tuning mechanisms, such as application of pressure, magnetic fields, gating or doping may be applied to effectively alter the band structure, which our work should provide a helpful reference to. Optical transitions were assigned to the specific bands based on the band structure of the compound. Several common features were found in the related compounds, CoSi and RhSi, based on the experimental reports: A sharp, prominent mid-infrared absorption, which originates from the parallel bands between the Γ and *R*-points and the successive linear-in-frequency increase. In addition, PdGa seems to possess a low lying optical transition at around 150–300 meV that is predicted to be absent in the sister compounds. Indeed, it has not been reported for either CoSi or for RhSi in previous experimental studies [35,48].

Furthermore, a linear-in-frequency section of the optical conductivity, often taken as a hint towards the presence of Dirac or Weyl fermions, around 0.9 eV to 1.5 eV has been identified, excluding the origin as transitions between the multifold chiral points intersecting Γ at around 0.5 eV and 1.0 eV. Nearly degenerate bands between Γ − *X* and Γ − *M* may hint at accidental Dirac-like degeneracies near the high-symmetry lines, warranting further investigation.

**Author Contributions:** All authors contributed to the calculations and the writing of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was supported by the Deutsche Forschungsgesellschaft (DFG) via DR228/51-3. E.U. acknowledges the European Social Fund and the Baden-Württemberg Stiftung for the financial support of this research project by the Eliteprogramme.

**Acknowledgments:** We thank Alexander A. Tsirlin for insightful discussions.

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

#### **References**


## *Article* **Effect of Deformation on Topological Properties of Cobalt Monosilicide**

**Sergey Nikolaev †, Dmitry Pshenay-Severin \*,†, Yuri Ivanov † and Alexander Burkov †**

Ioffe Institute, 26 Politekhnicheskaya, 194021 St Petersburg, Russia; sergey.nikolaev.w@mail.ru (S.N.); yu.ivanov@mail.ioffe.ru (Y.I.); a.burkov@mail.ioffe.ru (A.B.)

**\*** Correspondence: d.pshenay@mail.ru

† These authors contributed equally to this work.

**Abstract:** Recently, it was shown that materials with certain crystal structures can exhibit multifold band crossings with large topological charges. CoSi is one such material that belongs to non-centrosymmetric space group P213 (#198) and posseses multifold band crossing points with a topological charge of 4. The change of crystal symmetry, e.g., by means of external stress, can lift the degeneracy and change its topological properties. In the present work, the influence of uniaxial deformation on the band structure and topological properties of CoSi is investigated on the base of ab initio calculations. The **k** · **p** Hamiltonian taking into account deformation is constructed on the base of symmetry consideration near the Γ and *R* points both with and without spin-orbit coupling. The transformation of multifold band crossings into nodes of other types with different topological charges, their shift both in energy and in reciprocal space and the tilt of dispersion around nodes are studied in detail depending on the direction of uniaxial deformation.

**Keywords:** topological semimetal; cobalt monosilicide; mechanical deformation

**Citation:** Nikolaev S.; Pshenay-Severin D.; Ivanov Y.; Burkov A. Effect of deformation on topological properties of cobalt monosilicide. *Crystals* **2021**, *11*, 143. https://doi.org/10.3390/cryst11020143

Academic Editor: Artem Pronin Received: 30 December 2020 Accepted: 26 January 2021 Published: 29 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Cobalt monosilicide crystallizes in the cubic noncentrosymmetric space group #198 (P213). The unit cell and the Brillouin zone of CoSi are shown in Figure 1a,b. The band structure, magnetic, optical, transport and, in particular, thermoelectric properties of CoSi have been extensively studied [1–15]. Initially, the interpretation of experimental results was based on a simple semimetallic band structure model with small energy overlap of parabolic valence and conduction bands [1,2]. With the development of firstprinciple density functional theory (DFT) methods, a more realistic CoSi band structure has emerged [3–8]. Band structures calculated with and without the account of spin-orbit coupling (SOC) are plotted in Figure 1c,d, respectively. Earlier calculations [3,4] without the account of SOC revealed the presence of multiple band crossings at the Γ and *R* points of the Brillouin zone, but they did not consider the topology of the band structure. The symmetry analysis allowed to predict the existence of chiral fermions and multifold band crossings with high topological charges in crystals belonging to several space groups (including space group #198) in the presence of time-reversal symmetry [16,17]. In CoSi, multifold linear band crossing and spin texture was initially investigated around the Γ point, based on firstprinciple fully-relativistic calculations [5]. Later, detailed studies of band structure topology were made for CoSi [6,8] and for isostructural RhSi [7]. Effective **k** · **p** Hamiltonians around the time-reversal invariant momentum (TRIM) points were written down in Ref. [17] for *R* point and in Ref. [8] for Γ point. It was shown [6,8] that the topological charges at the Γ and *R* points are equal to ±4 and there are four surface Fermi arcs, connecting projections of these points on the surface Brillouin zone. Because spin-orbit coupling in CoSi is not strong, the Chern numbers were also calculated without SOC [6,18]. It was shown that multifold nodes have large topological charges of ±2 even without SOC (see Figure 1d

for illustration). The existence of multifold fermions and surface Fermi arcs in CoSi was recently confirmed by angle resolved photoemission spectroscopy (ARPES) [18–20].

**Figure 1.** The unit cell (**a**), the Brillouin zone (**b**) and the band structure of CoSi, calculated with (**c**) and without (**d**) the account of spin-orbit coupling. Insets in (**c**,**d**) show band structure around multi-fold band crossings at the Γ and *R* points (numbers are topological charges).

The theoretical study of the band structure beyond DFT, taking into account dynamic on-site correlations of d-electrons, revealed that, in contrast to FeSi, the electronic states in CoSi are only moderately influenced by electronic correlations [21]. Band broadening in CoSi is small in the range of ±0.3 eV near the Fermi level and decreases with the temperature. Thus, DFT description of CoSi band structure should give quite accurate results, that is confirmed by ARPES experiments [18–20] and by the better agreement of calculated lattice constants and elastic modules with experimental results for CoSi, compared to other monosilicides of the elements of the 4th period [22].

New information on the band structure of CoSi prompted to study the manifestation of its non-trivial topology and provided a base for correct interpretation of experimental results on conventional transport properties of the compound. For example, the account of real band structure and energy dependent relaxation time allowed to adequately explain the concentration dependencies of thermoelectric and galvanomagnetic properties of Co1−*x*Fe*x*Si and Co1−*x*Ni*x*Si alloys [10–12]. Recently, the effects of nonstoichiometry of CoSi-based materials on thermoelectric [9] and magnetic [14] properties were studied. In particular, in samples with the excess of Co, magnetically ordered states with helical and skyrmionic spin structures were observed near room temperature [14]. Quantum oscillations of thermopower with a beating pattern were observed in high-quality CoSi crystals [13]. They were successfully interpreted by the coexistence of two close Fermi surfaces in agreement with DFT results for the band structure. The influence of chiral fermions and charge density waves on magnetic field dependent electrical transport was studied in [23]. The experimental and theoretical investigation of optical conductivity of CoSi revealed various exotic multifold quasiparticles [15]. Moreover, low-frequency part of optical conductivity spectrum confirmed the existence of previously experimentally unobserved four-fold spin-3/2 node at the Γ point [15].

As the features of the band structure topology are due to particular crystal symmetry of CoSi, it is interesting to investigate the evolution of these properties when the symmetry changes. Such changes can appear due to mechanical stress, for example, in thin film devices or experimental setups, and can be important for device operation or interpretation of experimental data. In addition to the change of symmetry, mechanical deformation, in principle, can lead to the opening of a gap in the energy spectrum and the disappearance of the topological nodes. The possibility of using CMOS-compatible CoSi thin films for thermoelectric and sensor applications were considered recently in Ref. [24]. The stability of CoSi under hydrostatic pressure was theoretically investigated in Ref. [25], where it was predicted that the transition to CsCl structure (*Pm*3¯*m*) take place at hight pressure of 270 GPa. In the present work, we theoreticaly investigate another possibility—the change of band structure under uniaxial strain. In contrast to isotropic strain, uniaxial deformation changes the crystal symmetry even at low pressure. We considered deformation in [100], [110] and [111] directions. Based on symmetry analysis, the **k** · **p** Hamiltonian, taking into account deformation, was constructed for both the Γ and *R* points. Combining ab initio calculations, analytical model and symmetry considerations, the band splitting at the Γ and *R* points, the types of nodes arising from multifold band crossings and their energy and *k*-space positions were carefully studied both with and without SOC.

#### **2. Method of Calculation**

DFT calculations were performed in an integrated suite of Open-Source computer codes for electronic-structure calculations—Quantum ESPRESSO (QE) [26], using fully relativistic optimized normconserving Vanderbilt pseudopotentials (ONCV) [27]. The plane wave cut-off energy was 80 eV. The calculations were performed on 8 × 8 × 8 Monkhorst–Pack(MP) grid with the optimized lattice parameter *a*<sup>0</sup> = 4.438 Å. Four atomic positions of each of atomic species in the unit cell of undeformed CoSi are (*xA*, *xA*, *xA*), (−*xA* + 1/2, −*xA*, *xA* + 1/2), (−*xA*, *xA* + 1/2, −*xA* + 1/2), and (*xA* + 1/2, −*xA* + 1/2, −*xA*). Their optimization gives *x*Co = 0.144, *x*Si = 0.843.

Under uniaxial deformation, we set the unit cell parameters based on corresponding strain tensor and performed the relaxation of atomic positions, that allowed to determine the space group of deformed crystal.

For detailed study of the band structure, we performed Wannier interpolation using Wannier90 [28]. The position of nodes, topological charges and Fermi arcs were calculated using WannierTools [29] software package.

In order to analyze low-energy excitations around the nodes at Γ and *R* points, we constructed **<sup>k</sup>** · **<sup>p</sup>** Hamiltonian *<sup>H</sup>*<sup>ˆ</sup> in the presence of deformation from symmetry considerations. This allowed to independently verify the position of nodes and topological charges. As the effects of strain was assumed to be small, we considered only zeroth order in *k* terms in Hamiltonian proportional to strain tensor [30,31]: *Dijij*, where *ij* are strain tensor components and *Dij* are deformation potential parameters. The independent terms can be identified, applying symmetry operations of the considered space group (P213, #198), as it was made for the construction of Hamiltonian without stain [8,16,17]. We took into account that *ij* is transformed under symmetry operations as a product of wave vector components *kikj*, and used irreducible representations of space (double space) groups from Bilbao Crystallographic Server [32] for the case without (with) spin-orbit coupling. Since the spin-orbit coupling is also small in CoSi, we did not consider terms in *H*ˆ that depend on both strain and spin-orbit interaction. As will be seen from what follows, this approximation is sufficient.

The form of obtained Hamiltonians and their parameters are given in Appendix A. In the equations we used eV as units of energy. The wave vector components *ki*, *i* = 1, 2, 3 (in crystal coordinates) are measured in fractions of the reciprocal lattice vectors. The wave vector components *ki*, *i* = *x*, *y*, *z* (in Cartesian coordinates) are measured in the units of 2*π*/*a*<sup>0</sup> = 1.416 Å<sup>−</sup>1. If not stated otherwise, the latter units were used in band structure plots.

The deformation potential parameters were obtained using shifts of energy levels from ab initio calculations of undeformed and deformed crystal. The deformation potential parameters at the diagonal elements of strain tensor *ii* in **k** · **p** Hamiltonian determine the absolute shift of energy levels upon deformation. For CoSi, as metallic material, the absolute shift of energy level *<sup>n</sup>* due to deformation can be calculated as Δ*<sup>n</sup>* = ( (*d*) *<sup>n</sup>* <sup>−</sup> (*d*) *<sup>F</sup>* ) − ( (*u*) *<sup>n</sup>* <sup>−</sup> (*u*) *<sup>F</sup>* ), where (*d*(*u*)) *<sup>n</sup>* are energy levels in deformed (undeformed) crystal, and (*u*(*d*)) *<sup>F</sup>* are corresponding Fermi levels (see Appendix B for details).

#### **3. Results without SOC**

The band structure of CoSi without SOC features a triply-degenerate energy level at the Γ point close to the Fermi level. It is plotted in the inset of the Figure 1d and in the Figure 2 with dotted lines. The wave functions are transformed according to the three-dimensional single-valued representation Γ<sup>4</sup> of the little group of the Γ point of P213 (#198). The low-energy excitations around this point can be considered as effective spin-1 quasiparticles [33]. The topological charges of lower and upper linear branches are −2 and 2 respectively, while nearly flat band has zero charge.

**Figure 2.** The splitting of energy levels around the Γ point under uniaxial strain along [100] direction. Dotted lines represent the spectrum of undeformed crystal, solid (dashed) lines represent the spectrum in the case of compressive (tensile) strain. The absolute value of strain is |*e*| = 0.01. The wave vectors are measured in 2*π*/*a*<sup>0</sup> units.

Let us consider the simplest deformation of a crystal along the crystallographic direction [100]. When stretched along this direction, the spatial symmetry of group P213 (#198) is lowered to P212121 (#19), and essential remaining symmetry elements are: {*C*2*x*<sup>|</sup> <sup>1</sup> <sup>2</sup> <sup>0</sup><sup>1</sup> 2 }, {*C*2*y*<sup>|</sup> <sup>1</sup> 2 1 <sup>2</sup> <sup>0</sup>}, {*C*2*z*|0<sup>1</sup> 2 1 <sup>2</sup> }. At the Γ point, this group has only one-dimensional irreducible representations and, using character theory, we can expect that the three-dimensional representation Γ<sup>4</sup> of the space group P213 (#198) splits into three one-dimensional representations

of P212121 (#19) as: Γ<sup>4</sup> → Γ<sup>2</sup> + Γ<sup>3</sup> + Γ4. As Γ2(3,4) are real single-valued representations, they are not combined due to time-reversal symmetry (TRS) [32].

Thus, without taking into account the spin-orbit interaction, the triply-degenerate level at the Γ point is split into three with different energies. The low-energy band structure around the Γ point in CoSi under 1% uniaxial deformation in [100] direction is shown in the Figure 2 with solid lines for compressive strain and with dashed lines for tensile strain.

As the time-reversal symmetry is preserved in the considered cases, the whole energy spectrum around the TRIM point is symmetric with respect to the change of the sign of the wave vector **k** → −**k**. In addition, our **k** · **p** Hamiltonian is linear both in **k** and in components of deformation tensor ˆ. Hence, the low-energy band structure for stretched crystal can be obtained from the band structure of compressed crystal by changing the sign of energy (**k**, −ˆ) = −(**k**, ˆ), that can be seen by comparing solid and dashed curves in the Figure 2. It should be noted that this conclusion applies to all considered cased with the exception of eight-band Hamiltonian around the *R* point including SOC, as the latter contains terms independent of both wave vector and deformation tensor.

In compressed crystal there are nodes shifted in *kz* and *ky* directions upwards and downwards in energy relative to unstrained case, as shown in Figure 2. The topological charges of both nodes are equal to ±1. In stretched crystal, the sign of the energy changes, and the nodes swap. In both cases, we have two doubly degenerate (spinless) nodes at the same energy at the positions ±*knz* with the total topological charge of ±2 and similar nodes at positions ±*kny*. In the case, when only <sup>11</sup> = *e* is not zero, the node positions for small deformation can be obtained from eigenvalues of **k** · **p** Hamiltonian, and are equal to *knz* = (*D*<sup>1</sup> − *<sup>D</sup>*2)(*D*<sup>3</sup> − *<sup>D</sup>*2)*e*/*<sup>v</sup>* and *kny* = (*D*<sup>3</sup> − *<sup>D</sup>*1)(*D*<sup>3</sup> − *<sup>D</sup>*2)*e*/*v*. In these expressions *v* is the Fermi velocity at the Γ point in unstrained crystal and *Di* are the deformation potential parameters, defined in Appendix A after Equations (A1) and (A2).

It can be seen also, that the dispersion around doubly degenerate nodal points is tilted. As was shown in the Ref. [34], the general form of the Hamiltonian for Weyl point is the following: *H*(**k**) = ∑*i*,*<sup>j</sup> kiAijσj*, where *Aij* is a 3 × 4 matrix of coefficients and *σ<sup>j</sup>* are the 2 × 2 unit matrix and the three Pauli matrices for *j* = 0 and *j* = 1, 2, 3 respectively. The spectrum can be written as <sup>±</sup>(**k**) = *<sup>T</sup>*(**k**) <sup>±</sup> *<sup>U</sup>*(**k**), where *<sup>T</sup>*(**k**) = <sup>∑</sup><sup>3</sup> *<sup>i</sup>*=<sup>1</sup> *kiAi*<sup>0</sup> and *<sup>U</sup>*(**k**) = ∑3 *j*=1 - ∑3 *<sup>i</sup>*=<sup>1</sup> *kiAij*<sup>2</sup> . The nodal point is of the type II if there is a direction, in which *T*(**k**) > *U*(**k**). In the present case, it can be shown that in linear approximation *T*(**k**) = *U*(**k**) in *k*<sup>010</sup> or *k*<sup>001</sup> directions independently of the magnitude of strain *e*. Thus, under strain the nodal points are at the border of transition from type I to type II nodal points. However, it should be emphasized that, in contrast to ordinary type-II Weyl fermions, fermion states discussed here are spin degenerate.

The shift of nodes in reciprocal space implies the modification of the shape of surface Fermi arcs, that should emanate from the projections of the nodal points on the (100) surface Brillouin zone. In CoSi, it appeared that due to large extension of the Fermi arcs between the projections of the Γ and *R* point, their general shape changes quite moderately (at the scale of the full Brillouin zone) compared to undeformed case, presented in Refs. [6,8]. Therefore here we illustrate their variation around the Γ point only for several selected cases. The Figure 3 shows the Fermi arcs in the (100) surface Brillouin zone for the case of compressive deformation in the [100] direction for *e* = −0.01. For better visualisation, Fermi level was shifted to the energy of nodal points. The position of nodes, obtained by **k** · **p** calculations, are plotted in the figure by black asterisks. It can be seen, that two nodes are shifted along *kz* direction, and their positions correlate with the sources of two surface Fermi arcs.

When deforming in [110] direction, off-diagonal elements of the strain tensor begin to play a role. The spatial symmetry of deformed cobalt monosilicide is described by the space group P21 (#4) with the only symmetry element {*C*2*z*|00<sup>1</sup> <sup>2</sup> } (except lattice translations). Therefore, without taking into account the spin-orbit interaction, the three-dimensional representation Γ<sup>4</sup> of the group *P*213(#198) splits into three one-dimensional representations

of *P*21(#4) as: Γ<sup>4</sup> → Γ<sup>1</sup> + 2Γ2. In this way, the degeneracy at the Γ point is completely lifted, as in the case of the deformation along the main unit cell directions. In the most simple case, compatible with considered symmetry, where the only non-zero components of stain tensor are <sup>12</sup> = <sup>21</sup> = *e*/2, two effective spin-1/2 nodes symmetrically diverge from Γ point along [110] ([11¯0]) crystallographic direction and shift to lower energies in the case of compressive (tensile) strain *e* < 0 (*e* > 0). Similar nodes appear at higher energies but they are shifted along [11¯0] ([110]) directions for *e* < 0 (*e* > 0). In more general case, when <sup>11</sup> = <sup>22</sup> = <sup>12</sup> = *e*/2, that corresponds to the absence of deformation in directions normal to [110], the nodes split along the line, that is rotated by a small angle (about *φ* = 5◦ at *e* = 0.01) from [110] ([11¯0]) axis. Let's denote these directions by *k φ* <sup>110</sup> (*k φ* 110¯ ). The low-energy band structure for these directions and topological charges are given in the Figure 4. The situation is somewhat similar to [100] case (see, Figure 2), but the shift of nodes in both energy and k-space are larger. In this case we obtain again effective spin-1/2 nodes with the tilt intermediate between tilts of the type-I and type-II nodes. The total topological charge of each pair of the nodes is ±2.

**Figure 3.** The details of surface Fermi arcs around the Γ point under uniaxial strain along [100] direction (*e* = −0.01) without the account of spin-orbit coupling (SOC). Asterisks depict the positions of nodal points.

**Figure 4.** The splitting of energy levels around the Γ point under uniaxial stain along [110] direction. Dotted lines represent the spectrum of undeformed crystal, solid (dashed) lines represent the spectrum in the case of compressive (tensile) strain. The absolute value of strain is |*e*| = 0.01.

The deformation along the [111] crystallographic direction needs special consideration for CoSi. When deformed in [111] direction, the symmetry of cobalt monosilicide is reduced to *R*3(#146) space group with the single 3-fold rotation axis. Character theory suggests that without SOC, the three-dimensional representation Γ<sup>4</sup> of *P*213(#198) space group splits into three one-dimensional representations of *R*3(#146) as: Γ<sup>4</sup> → Γ<sup>1</sup> + (Γ<sup>2</sup> + Γ3). Since the representations Γ<sup>2</sup> and Γ<sup>3</sup> are mutually conjugate, they should be combined due to time reversal symmetry. The triply degenerate level is split into two (nondegenerate and doubly degenerate) levels, in contrast to other types of deformation in which the degeneracy is completely lifted at the Γ point. In this case, strain tensor was taken in the form *ij* = *e*/3 for all *i*, *j*. In deformed crystal, one unusual node is located at the Γ point and two other nodes are displaced along the [111] direction at the positions ±*kn*,111 with *kn*,111 = *D*4*e*/*v*(1 + *e*) (see the Figure 5). In this case we obtained one node with topological charge ±2 at the Γ point and two nodes with charges ±1 at ±*kn*,111 points. The spectrum for compressive and tensile strains are again can be obtained by the energy sign change. Thus the number of nodes between the two upper bands depends on the sign of deformation *e*. Another way to see this result is to plot Fermi arcs in the (001) surface Brillouin zone for compressive and tensile strain (see Figure 6). Consider a compressed crystal. Below the Fermi level, there are two nodes shifted in [111] direction. Their projections on (001) plane are sources of two Fermi arcs (left panel). The starting points of the arcs are shifted towards the projections of the nodes (black asterisks) but do not coincide with them exactly, since the Fermi level is located above these nodes. In the case of expansion, the node with topological charge of 2 is below the Fermi level, and both Fermi arcs start from the Γ point (right panel).

**Figure 5.** The splitting of energy levels around the Γ point under uniaxial stain along [111] direction. Dotted lines represent the spectrum of undeformed crystal, solid (dashed) lines represent the spectrum in the case of compressive (tensile) strain. The absolute value of strain is |*e*| = 0.01.

**Figure 6.** The details of surface Fermi arcs around the Γ point under uniaxial compression (**left** panel) and extension (**right** panel) along [111] direction (|*e*| = 0.01) without the account of SOC. Asterisks depict the position of nodal points.

The fermions around the node with topological charge of 2 at the Γ point are similar to that of quadratic double-Weyl fermions in SrSi2, described in Ref. [35], but they arise for

another reason. In SrSi2 without SOC, on the four-fold rotation axes of the crystal, there are two-fold band crossings with linear dispersion. The bands are spin degenerate. The account of spin-orbit coupling adds terms independent of wave vector to the Hamiltonian that leads to a partial lifting of the degeneracy and to a change of linear dispersion to quadratic in certain directions. In CoSi, the k-independent terms in *H*ˆ appear due to strain even without SOC. Counting energy relative to the charge-2 node in the Figure 5, the dispersion in [111] direction remains linear: Δ = ±*k*111*v*(1 + *e*). While, in perpendicular direction, it becomes quadratic, Δ = *D*4*e*/2 − sign(*e*) (*D*4*e*/2)<sup>2</sup> + (*k*⊥,111*v*)<sup>2</sup> or flat. Thus, in CoSi under [111] deformation, there is a strain induced transition from spin-1 quasiparticles to quadratic effective-spin-1/2 fermions.

Let us now consider low-energy band structure around the *R* point without spin-orbit coupling. In this case in unstrained CoSi, the energy level at the *R* point is four-fold degenerate not considering spin (see Figure 1d). Low energy excitations around the *R* point are double spin-1/2 fermions [33] with the Chern number −2 [6]. The wave functions are transformed according to the direct sum of single-valued mutually conjugated complex two-dimensional representations *R*<sup>1</sup> and *R*3, combined due to TRS. The switching to *P*212121 (#19) space group under [100] uniaxial deformation does not lead to the energy level splitting, because the representation *R*<sup>1</sup> + *R*<sup>3</sup> is transformed into the direct sum of the pseudoreal representations *R*<sup>1</sup> + *R*<sup>1</sup> of *P*212121. Thus the node at the *R* point remains intact.

Under uniaxial [110] stress, the space group *P*213 is reduced to *P*21 group and the *R* point goes into the *E* point, which is also located at the vertex of the deformed Brillouin zone. The representation *R*<sup>1</sup> + *R*<sup>3</sup> is transformed into 2(*E*<sup>1</sup> + *E*2). The *E*<sup>1</sup> and *E*<sup>2</sup> representations of *P*21 are one-dimensional and mutually conjugated, therefore they are combined due to TRS. Thus, four-fold degenerate energy level splits into two doubly degenerated levels (see Figure 7). The nodes are split along *kz* direction and are situated at *knz* = ± *D*<sup>2</sup> <sup>2</sup> + *<sup>D</sup>*<sup>2</sup> <sup>3</sup>*e*/2*v*, where deformation potential constants are given in Appendix A after Equation (A3). The topological charges of these two nodes are ∓1. Two crossing points at the *R* point (*E* point of P212121 space group) in the Figure 7 are not nodes as the states are degenerate on the (001) surface of the deformed Brillouin zone.

**Figure 7.** The splitting of energy levels around the *R* point under uniaxial stain along [110] direction (*E* point, P21 space group). Dotted lines represent the spectrum of undeformed crystal, solid lines represent the spectrum in the case of compressive strain. The value of strain is *e* = −0.01.

In the case of [111] deformation, the four-fold degenerate level at the *R* point splits into three levels (two nondegenerate, and one doubly degenerate) at the corner (*T* point) of deformed Brillouin zone (see Figure 8). Representations transform according to the expression *R*<sup>1</sup> + *R*<sup>3</sup> → 2*T*<sup>1</sup> + (*T*<sup>2</sup> + *T*3). All representations at the *T* point are one-dimensional, but *T*<sup>2</sup> and *T*<sup>3</sup> are mutually conjugated and should be combined due to TRS. Similarly to the case of Γ point, doubly degenerate node at the *T* point has topological charge of ∓2 and has quadratic dispersion in the direction perpendicular to [111]. There are also nodes

with tilted dispersion, shifted in the [111] direction by *kn*,111 = ±(2*D*<sup>2</sup> <sup>2</sup> − *<sup>D</sup>*<sup>2</sup> <sup>3</sup>)*e*/2 <sup>√</sup>3*D*2*v*. Their topological charge is ∓1.

**Figure 8.** The splitting of energy levels around the *R* point under uniaxial stain along [111] direction (*T* point of *R*3 space group). Dotted lines represent the spectrum of undeformed crystal, solid (dashed) lines represent the spectrum in the case of compressive (tensile) strain. The absolute value of strain is |*e*| = 0.01.

#### **4. Results with SOC**

Taking into account the spin-orbital coupling in the crystal without deformation, the 6-fold degenerate level at the Γ point splits into a doublet and 4-fold degenerate levels. Their wave functions are transformed according to the Γ<sup>5</sup> irreducible representation and mutually conjugated Γ<sup>6</sup> and Γ<sup>7</sup> irreducible representations, combined due to timereversal symmetry [8]. The low-energy excitations around the 4-fold degenerate node at the Γ point are spin-3/2 fermions [6,33] with topological charge of 4. Its dispersion is shown in the Figure 1c, and is also plotted with dotted lines in the Figure 9. The doublet does not move from the Γ point due to deformation, so we consider in more details the evolution of fourfold node at the Γ point under the influence of deformation along the main crystallographic directions.

**Figure 9.** The splitting of energy levels around the Γ point under uniaxial stain along [100] direction. Dotted lines represent the spectrum of undeformed crystal, solid lines represent the spectrum in the case of compressive strain with *e* = −0.01.

Under the deformation along [100] the fourfold degenerate level splits into two twofold degenerate levels, both corresponding to Γ<sup>5</sup> representation of the P212121 (#19) group, as this two-dimensional representation is time-reversal invariant. The low-energy Hamiltonian for these levels is given by Equations (A4) and (A5) of Appendix A. When the

only non-zero component of deformation is <sup>11</sup> = *e*, the gap between two levels at the Γ point is equal to 2|*D*3|*e*. In general case, it is equal to 2|*D*3|*e*(1 + *νP*), where *ν<sup>P</sup>* is the Poisson ratio. The low-energy spectrum is shown in the Figure 9. The doublets at the Γ point form Weyl nodes with unit topological charge. Four nodes are shifted from the Γ point along the diagonals of the *ky* − *kz* plane and each of them has a unit topological charge. In addition, two Weyl nodes shifted from the Γ point in the ±*ky* directions appear between two lower branches of the spectrum, and two similar Weyl nodes shifted in the ±*kz* directions (not shown in Figure 9) appear between two upper branches of the spectrum. The tensile and compressive strain spectra differ from each other by a change in the sign of the energy.

Without deformation, there are four surface Fermi arcs, emanating from the projection of the nodal point at Γ, as the topological charge is 4. The shift of the node positions for [100] deformation (*e* = −0.01) is accompanied by a change in the surface Fermi arcs, shown in the Figure 10, for the case when the Fermi level coincides with the nodes located along the diagonals of the *ky* − *kz* plane. The projections of the nodes are marked with asterisks. Their positions, calculated by the **k** · **p** method, quite well coincide with the sources of four Fermi arcs.

**Figure 10.** The details of surface Fermi arcs around the Γ point under uniaxial strain along [100] direction (*e* = −0.01) with the account of SOC. Asterisks depict the positions of nodal points.

Considering the deformation along [110] axis, we found similar splitting of four-fold degenerate level at the Γ point into 2 doublets: Γ<sup>6</sup> + Γ<sup>7</sup> → 2(Γ<sup>3</sup> + Γ4), where irreducible representations Γ3(4) of the space group *P*21(#4) are one-dimensional and combined together due to TRS. The energy gap at the Γ point is equal to *e D*<sup>2</sup> <sup>2</sup> + *<sup>D</sup>*<sup>2</sup> <sup>3</sup>. When the only non-zero components of stain tensor are <sup>12</sup> = <sup>21</sup> = *e*/2, the four-fold degenerate node splits into 4 nodes, moving along crystallographic directions [100], [100], [010], [010] by the distance *D*2*e*/2 <sup>√</sup>*b*<sup>2</sup> <sup>−</sup> *<sup>a</sup>*2. Due to distortion, small deviation from corresponding Cartesian axis appears by an angle of *φ* = arctan(*e*/2), which is equal, for example, to 0.3◦ at *e* = 0.01. In the case of the absence of deformations in the directions, normal to [110], when <sup>11</sup> = <sup>22</sup> = <sup>12</sup> = *e*/2, the shift of energy levels due to volume change leads to additional deviation of nodes from Cartesian axes, which is about *φ* = 6◦ at *e* = 0.01. Similar to the case of [100] deformations, there are also two Weyl nodes below and two Weyl nodes above the node in unstrained crystal, but here they are shifted close to the diagonals of *kx* − *ky* plane. To within a change of the directions of node shifts, the low-energy spectrum around the Γ point for this case is qualitatively very similar to the case of [100] deformations (see the Figure 9).

When the stress is applied along [111] axis, we obtain similar 4-fold level slitting into two doublets at the Γ point with Γ<sup>6</sup> + Γ<sup>7</sup> → (Γ<sup>4</sup> + Γ4)+(Γ<sup>5</sup> + Γ6), where all irreducible representations Γ4(5,6) of the space group *R*3 (#146) are one-dimensional. Γ<sup>4</sup> is real and it is doubled due to TRS, while Γ5(6) are complex conjugated and they are combined together due to TRS. These two band crossings at the Γ point are shown in the Figure 11. The lower one is a Weyl node. The node at higher energy has linear dispersion in the [111] direction and nonlinear dispersion in the perpendicular plane (this branch is not shown in the figure). It looks like the triple-Weyl node [36–38] with topological charge ± 3, but it is located at the TRIM point. In addition, near the Γ point, there are two groups of nodes connected by timereversal symmetry. Each of these groups consist of 4 nodes. One of them *kn*<sup>1</sup> is shifted from the <sup>Γ</sup> point along the [111] axis by a distance *<sup>k</sup>* <sup>=</sup> *<sup>D</sup>*2*e*(*<sup>a</sup>* <sup>+</sup> <sup>√</sup>*a*<sup>2</sup> <sup>+</sup> <sup>2</sup>*b*2)/ <sup>√</sup>3*b*2, which is about 0.007 at |*e*| = 0.01 (see Figure 11). In the case of compressive (tensile) strain the energy of the node shifts downwards (upwards) relative to the node at the Γ point in unstrained crystal. More detailed calculations showed that there are another 3 nodes in each group: one of them *kn*<sup>2</sup> is shifted from *kn*<sup>1</sup> into [112] direction and positions of two other nodes *kn*3(4) can be obtained using 2*π*/3 rotation around [111] axis. The distance between nodes is rather small, about 0.0003 at *e* = −0.01. The calculation of topological charge showed that three equivalent nodes *kn*2(3,4) have topological charge of 1 each, and the node *kn*<sup>1</sup> has a charge of −1. Thus total topological charge of each group is 2, giving 4 for both groups together. The dispersion around the nodes is rather complex and is given in the Figure 12. It can be seen that the electronic velocities are very different in different directions, and the nodes are tilted. In addition to these groups of nodes, there are also two tilted Weyl nodes in the directions [111] and [1¯1¯1¯]. The lowest band crossing in Figure 11 is one of these node.

**Figure 11.** The splitting of energy levels around the Γ point under uniaxial stain along [111] direction. Dotted lines represent the spectrum of undeformed crystal, solid lines represent the spectrum in the case of compressive strain for *e* = −0.01.

**Figure 12.** The dispersion around the nodes shifted from the Γ point under uniaxial compressive strain along [111] direction. The dispersions along [111] (solid lines), [110] (dot-dashed lines) and [112] (towards the node *kn*2, dashed lines) directions are plotted for *e* = −0.01.

Let us now consider the effect of deformation on the low-energy spectrum around the *R* point in the presence of spin-orbit coupling. In unstrained crystal, the 8-fold degenerate level at the *R* point splits into 6-fold degenerate level *R*<sup>7</sup> + *R*<sup>7</sup> (double spin-1 quasiparticles) and doublet *R*<sup>5</sup> + *R*<sup>6</sup> due to SOC. The doublet is not split by the strain, so we will consider only 6-fold degenerate level (see the Figure 1c). In the case of [100] deformation, it is split into 3 doublets with the wave functions, transforming according to the representations *R*<sup>2</sup> + *R*2, *R*<sup>3</sup> + *R*<sup>3</sup> and *R*<sup>4</sup> + *R*<sup>4</sup> of the P212121 (#19) group. The change of the dispersion around the *R* point is shown in the Figure 13. Two nodes at lower energy are shifted along positive and negative *kx* (*kz*) directions under compressive (tensile) strain. There is a crossing of four energy branches at each of these nodal points. They are a tilted double spin-1/2 nodes [33] with Chern numbers of ∓2. In the case of eight-band Hamiltonian, the equality (**k**, −ˆ) = −(**k**, ˆ) does not hold exactly. But in the case of [100] deformation the spectrum for six considered bands approximately follows this rule. This implies that there are the two additional double spin-1/2 nodes in *kz* (*kx*) directions for compressive (tensile) strain (see Figure 13).

**Figure 13.** The splitting of energy levels around the *R* point under uniaxial stain along [100] direction. Dotted lines represent the spectrum of undeformed crystal, solid (dashed) lines represent the spectrum in the case of compressive (tensile) strain of |*e*| = 0.01.

The deformation in [110] direction again leads to the splitting of 6-fold degenerate level into 3 doublets at the vertex of deformed Brillouin zone (*E* point). Wave functions of two doublets are transformed according to the *E*<sup>3</sup> + *E*<sup>3</sup> representation and one doublet according to the *E*<sup>4</sup> + *E*<sup>4</sup> representation, where both *E*3(4) are real one-dimensional. These doublets do not form Weyl nodes, as the energy branches, starting from them are degenerate at the edges of the Brillouin zone parallel to (001) plane. Four simple Weyl nodes are formed near the *E* point. They shifted mainly into *kz* direction to the points with coordinates ±*kn*, where *kn* = (±0.001, ∓0.001, 0.014) at *e* = −0.01 (see Figure 14). In the case of tensile strain the shift appeared to be almost the same. As the degeneracy of the bands is completely lifted under this deformation, the topological charge of each of the 4 nodes is ∓1.

The case of [111] deformation is similar to the two previously considered cases in the sense that one obtains 3 doublets instead of 6-fold degenerate level at the vertex of the Brillouin zone (2*R*<sup>7</sup> → (*T*<sup>4</sup> + *T*4) + 2(*T*<sup>5</sup> + *T*6), where *T*<sup>4</sup> is real one-dimensional representations, *T*<sup>5</sup> and *T*6, are complex conjugated one dimensional representations of the little group of the *T* point of the *R*3(#146) space group). Each of these doublets corresponds to a Weyl node. The two band crossings at the *T* point are the conventional Weyl nodes with Chern number of ∓1. The third band crossing looks like the triple-Weyl node with the topological charge of ∓3. It is similar to triple-Weyl node at the Γ point. In addition, near the *T* point, there are two groups of nodes connected by time-reversal symmetry, as in the case of the Γ point. Two of them are simple Weyl nodes shifted into [111] and [1¯1¯1¯] directions, but they have topological charge 1. Around each of them, there are three

conventional Weyl nodes with topological charges −1. The total charge of these eight nodes is −4. For example, in the case of compressive strain with *e* = −0.01, one of the nodes at the [111] direction have coordinates *kn*<sup>1</sup> = (0.00197, 0.00197, 0.00197) relative to the *T* point. One of its satellites have coordinates *kn*<sup>2</sup> = (0.00169, 0.00196, 0.00227), and the coordinates of another two satellites *kn*3(4) can be obtained by cyclic permutations. The dispersion around the *T* point towards nodes *kn*<sup>1</sup> and *kn*<sup>2</sup> is plotted in the Figure 15. It is almost linear. At the same time, the dispersion along the line connecting the central node with one of its satellites (Figure 16) looks like a result of the crossing of two nonlinear bands.

**Figure 14.** The splitting of energy levels around the *R* point (*E* point, P21 space group) under uniaxial stain along [110] direction. Dotted lines represent the spectrum of undeformed crystal, solid lines represent the spectrum in the case of compressive strain of *e* = −0.01.

**Figure 15.** The splitting of energy levels around the *R* point under uniaxial stain along [111] direction (*T* point, *R*3 space group). Dotted lines represent the spectrum of undeformed crystal, solid lines represent the spectrum in the case of compressive strain of *e* = −0.01.

**Figure 16.** The electronic dispersion along the direction from central node *kn*<sup>1</sup> towards to one of its satellites under uniaxial stain along [111] direction at *e* = −0.01.

#### **5. Conclusions**

In the present work the influence of deformation on the band structure and topological properties of CoSi was studied using both ab initio calculations, and symmetry considerations. The symmetry prescribed **k** · **p** Hamiltonians at the Γ and *R* TRIM points taking into account deformation were written down for the cases with and without SOC. It was shown that in almost all considered cases, the degeneracy is partially lifted at the TRIM points. The only exception is the fourfold degenerate level at the *R* point (without SOC) under [100] strain. A lowering in symmetry leads to the appearance of a significant number of different band crossings with topological charges from ±1 to ±3 around the TRIM points. The nodes often have a tilted dispersion.

The unusual results were obtained upon deformation of CoSi along the [111] direction. Without spin-orbit coupling, the doubly degenerate nodes with quadratic dispersion in the plane orthogonal to the [111] direction appear at the Γ and *T* points of the deformed Brillouin zone. These band crossings have Chern numbers of ±2 and resemble the wellknown double-Weyl nodes, but they are spin degenerate. Calculation with account of SOC revealed doubly degenerate nodes with the topological charges of ±3 at the TRIM points. These band crossings are located on the threefold rotation axis and are analogous to triple-Weyl nodes.

The band structure with SOC around the *R* point under [100] strain exhibits another example of the change of node type. The double spin-1 node with topological charge of 4 splits into pairs of double spin-1/2 nodes with topological charges of 2 per node. Thus, using mechanical deformation, the transition between different types of topological nodes can be realized in the same material.

A lowering of the crystal symmetry under strain also leads to a modification of the surface Fermi arcs shape. A change in the sign of the deformation and the Fermi level position switches the ends of the Fermi arcs from one group of nodes to another. However, the number of Fermi arcs always remains equal to two without taking into account SOC and four with SOC.

As a byproduct of low-energy Hamiltonian fitting, the absolute deformation potential parameters were obtained for considered energy states, and the work function of CoSi was calculated (4.55 eV), which correlates with available experimental data.

**Author Contributions:** Conceptualization, Y.I. and A.B.; methodology, Y.I.; investigation, S.N. and D.P.-S.; writing—original draft preparation, S.N. and D.P.-S.; writing—review and editing, Y.I. and A.B.; supervision, A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** The study was supported by the Russian Foundation for basic Reseach, project 18-52- 80005 (BRICS).

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

#### **Appendix A. k** · **p Hamiltonians and Their Parameters**

In this section, the form of **k** · **p** Hamiltonians and their parameters will be given. We used eV units for energy, and the dimensionless wave vector components *ki*, *i* = 1, 2, 3 are measured in fractions of the reciprocal lattice vectors.

Without spin-orbital coupling linear in wave vector part of Hamiltonian at the Γ point is given by the following equation:

$$H\_{\Gamma1} = \begin{pmatrix} 0 & ivk\_3 & -ivk\_2 \\ -ivk\_3 & 0 & ivk\_1 \\ ivk\_2 & -ivk\_1 & 0 \end{pmatrix} \tag{A1}$$

where *v* = 1.73 eV. The node at the Γ point lies 3.6 meV above the Fermi level.

The perturbation Hamiltonian in the linear approximation in the deformation tensor *ik* and in the zero approximation in the wave vector has the following form:

$$H\_{\Gamma2} = \begin{pmatrix} S\_1(\mathfrak{k}) & D\_4 \mathfrak{e}\_{12} & D\_4 \mathfrak{e}\_{13} \\ D\_4 \mathfrak{e}\_{12} & S\_2(\mathfrak{k}) & D\_4 \mathfrak{e}\_{23} \\ D\_4 \mathfrak{e}\_{13} & D\_4 \mathfrak{e}\_{23} & S\_3(\mathfrak{k}) \end{pmatrix} \tag{A2}$$

where *S*1(ˆ) = *D*<sup>1</sup><sup>11</sup> + *D*<sup>2</sup><sup>22</sup> + *D*<sup>3</sup>33, *S*2(ˆ) = *D*<sup>3</sup><sup>11</sup> + *D*<sup>1</sup><sup>22</sup> + *D*<sup>2</sup>33, *S*3(ˆ) = *D*<sup>2</sup><sup>11</sup> + *D*<sup>3</sup><sup>22</sup> + *D*<sup>1</sup>33, *D*<sup>1</sup> = −0.319 eV, *D*<sup>2</sup> = −0.400 eV, *D*<sup>3</sup> = 0.470 eV and *D*<sup>4</sup> = 2.479 eV.

At the *R* point the Hamiltonian, linear in k-vector, was given in Ref. [16]. It can be represented as *HR*<sup>1</sup> <sup>=</sup> *<sup>v</sup>*1ˆ <sup>⊗</sup> (*<sup>σ</sup>* · **<sup>k</sup>**) with *<sup>v</sup>* = 1.28 eV. The node at the *<sup>R</sup>* point lies 0.211 eV below the Fermi level. The perturbation due to elastic strain reads:

$$H\_{R2} = \begin{pmatrix} D\_1 \text{Tr}\mathfrak{k} - D\_2 \epsilon\_{12} & -D\_2(\epsilon\_{23} - i\epsilon\_{13}) & iD\_3\epsilon\_{12} & D\_3(\upsilon\_6^\* \varepsilon\_{23} + i\upsilon\_6 \varepsilon\_{13}) \\ -D\_2(\epsilon\_{23} + i\epsilon\_{13}) & D\_1 \text{Tr}\mathfrak{k} + D\_2 \epsilon\_{12} & D\_3(\upsilon\_6^\* \varepsilon\_{23} - i\upsilon\_6 \varepsilon\_{13}) & -iD\_3 \epsilon\_{12} \\ -i D\_3^\* \epsilon\_{12} & D\_3^\*(\upsilon\_6 \varepsilon\_{23} + i\upsilon\_6^\* \varepsilon\_{13}) & D\_1 \text{Tr}\mathfrak{k} + D\_2 \epsilon\_{12} & D\_2(\epsilon\_{23} - i\varepsilon\_{13}) \\ D\_3^\*(\upsilon\_6 \epsilon\_{23} - i\upsilon\_6^\* \varepsilon\_{13}) & iD\_3^\* \epsilon\_{12} & D\_2(\epsilon\_{23} + i\varepsilon\_{13}) & D\_1 \text{Tr}\mathfrak{k} - D\_2 \epsilon\_{12} \end{pmatrix} \tag{A3}$$

where Tr<sup>ˆ</sup> is a trace of strain tensor, *<sup>ν</sup>*<sup>6</sup> <sup>=</sup> *<sup>e</sup>i<sup>π</sup>*/6, *<sup>D</sup>*<sup>1</sup> = 0.264 eV, *<sup>D</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>1.29 eV and *D*<sup>3</sup> = 3.27 eV.

Taking into account spin-orbit coupling, linear Hamiltonian at the Γ point was written down in Ref. [8] and has the following form:

$$H\_{\Gamma1}^{(\text{SOC})} = \begin{pmatrix} ak\_3 & a(k\_1 - ik\_2) & b(\nu\_3 k\_1 - \nu\_6 k\_2) & bk\_3 \\ a(k\_1 + ik\_2) & -ak\_3 & bk\_3 & -b(\nu\_3 k\_1 + \nu\_6 k\_2) \\ b^\*(\nu\_3^\* k\_1 - \nu\_6^\* k\_2) & b^\*k\_3 & -ak\_3 & -a(k\_1 + ik\_2) \\ b^\*k\_3 & -b^\*(\nu\_3^\* k\_1 + \nu\_6^\* k\_2) & -a(k\_1 - ik\_2) & ak\_3 \end{pmatrix},\tag{A4}$$

where *ν*<sup>3</sup> = *ei<sup>π</sup>*/3, *a* = 0.56 eV and *b* = 1.19 eV [8]. The 4-fold degenerate node position is 21 meV above the Fermi level, and the Weyl cone at the Γ point is shifted down due to SOC by 54 meV relative to this node.

The perturbation due to deformation reads:

$$H\_{\Gamma2}^{(\text{SOC})} = \begin{pmatrix} D\_1 \text{Tr}\mathfrak{e} + D\_2 \mathfrak{e}\_{12} & D\_2(\mathfrak{e}\_{23} - i\mathfrak{e}\_{13}) & 0 & D\_3 \Sigma(\mathfrak{e}) \\ D\_2(\mathfrak{e}\_{23} + i\mathfrak{e}\_{13}) & D\_1 \text{Tr}\mathfrak{e} - D\_2 \mathfrak{e}\_{12} & -D\_3 \Sigma(\mathfrak{e}) & 0 \\ 0 & -D\_3^\* \Sigma^\*(\mathfrak{e}) & D\_1 \text{Tr}\mathfrak{e} + D\_2 \mathfrak{e}\_{12} & D\_2(\mathfrak{e}\_{23} + i\mathfrak{e}\_{13}) \\ D\_3^\* \Sigma^\*(\mathfrak{e}) & 0 & D\_2(\mathfrak{e}\_{23} - i\mathfrak{e}\_{13}) & D\_1 \text{Tr}\mathfrak{e} - D\_2 \mathfrak{e}\_{12} \end{pmatrix},\tag{A5}$$

where <sup>Σ</sup>(<sup>ˆ</sup>) = <sup>11</sup> − *<sup>ν</sup>*<sup>3</sup><sup>22</sup> + *<sup>ν</sup>*<sup>2</sup> <sup>3</sup> 33, *D*<sup>1</sup> = −0.085 eV, *D*<sup>2</sup> = 1.40 eV. Parameter *D*<sup>3</sup> is complex. Eigenvalues at *k* = 0 does not depend on its phase, but it affects the spectrum for nonzero *<sup>k</sup>* values. The fitting gives *<sup>D</sup>*<sup>3</sup> <sup>≈</sup> 0.233*e*−*i<sup>π</sup>*/6eV.

At the *R* point the Hamiltonian for 6-fold degenerate node including SOC was given in Ref. [17]. After uniaxial deformation this node splits into three doubly-degenerate nodes. Under deformation in [100] direction, the shift of energy levels at the *R* point is linear in deformation. If the deformation is applied in [111] direction, only for small *e* < 0.004, the shift of energy levels can be considered as linear, and the deformation along [110] axis leads to nonlinear shift of the two pairs of energy levels (see Figure A1). It was shown in Ref. [16] in the framework of a simple model that a linear Hamiltonian that includes SOC in the zeroth order with respect to the wave vector and takes into account all eight bands leads to the correct nonlinear band dispersion near the *R* point (see, e.g., Figure 3b–c in Ref. [8]). So, we consider both nodes together and obtain 8 × 8 Hamiltonian at the *R* point. The zero-order Hamiltonian has only nonzero matrix elements (*H*(*SOC*) *<sup>R</sup>*<sup>0</sup> )*ii* = −Δ, *i* = 7, 8, which describe energy shift of doublet downwards in energy due to SOC. Including SOC, the position of the 6-fold degenerate node is 0.202 eV below *<sup>F</sup>*, while the band splitting Δ = 32 meV. The zero- and linear-order in *k* parts together reads:

$$H^{(\text{SOC})}\_{\text{R01}} = \begin{pmatrix} 0 & a\_1 k\_3 & -a\_1^\* k\_2 & 0 & a\_2 k\_3 & -a\_2 k\_2 & a\_3 k\_1 & a\_4 k\_1 \\ a\_1^\* k\_3 & 0 & a\_1 k\_1 & a\_2 k\_3 & 0 & a\_2 k\_1 & v\_3 a\_3 k\_2 & v\_3^\* a\_4 k\_2 \\ -a\_1 k\_2 & a\_1^\* k\_1 & 0 & -a\_2 k\_2 & a\_2 k\_1 & 0 & -v\_3^\* a\_3 k\_3 & -v\_3 a\_4 k\_3 \\ 0 & a\_2^\* k\_3 & -a\_2^\* k\_2 & 0 & -a\_1^\* k\_3 & a\_1 k\_2 & a\_4^\* k\_1 & -a\_3^\* k\_1 \\ a\_2^\* k\_3 & 0 & a\_2^\* k\_1 & -a\_1 k\_3 & 0 & -a\_1^\* k\_1 & v\_3 a\_4^\* k\_2 & -v\_3^\* a\_4^\* k\_2 \\ -a\_2^\* k\_2 & a\_2^\* k\_1 & 0 & a\_1^\* k\_2 & -a\_1 k\_1 & 0 & -v\_3^\* a\_4^\* k\_3 & v\_3 a\_3^\* k\_3 \\ a\_3^\* k\_1 & v\_3^\* i\_3^\* k\_2 & -v\_3 a\_3^\* k\_3 & a\_1 k\_1 & v\_3^\* a\_2 k\_2 & -v\_3 a\_4 k\_3 & -\Delta & 0 \\ a\_4^\* k\_1 & v\_3 a\_4^\* k\_2 & -v\_3^\* a\_4^\* k\_3 & -a\_3 k\_1 & -v\_3 a\_3 k\_3 & v\_3^\* a\_3 k\_3 & 0 & -\Delta \end{pmatrix} \tag{A6}$$

where parameters *ai* are complex. Their values were obtained by fitting to electron spectrum around *R* point. They are not unique, but they were checked to give correct values of topological charges. These values are *a*<sup>1</sup> = (0.342 − 0.686*i*) eV, *a*<sup>2</sup> = (1.043 + 0.060*i*) eV, *a*<sup>3</sup> = (−0.459 − 0.657*i*) eV and *a*<sup>4</sup> = (0.181 − 0.100*i*) eV.

**Figure A1.** The splitting of energy levels at the *R* point with uniaxial stress of magnitude *e* in [110] (**left** panel) and in [111] (**right** panel) directions.

The perturbation due to deformation is:

*H*(*SOC*) *<sup>R</sup>*<sup>2</sup> = ⎛ ⎜⎜⎜⎜⎜⎝ *S*1(ˆ) *D*<sup>5</sup><sup>12</sup> −*D*<sup>∗</sup> <sup>5</sup> <sup>13</sup> 0 *D*<sup>6</sup><sup>12</sup> *D*<sup>6</sup><sup>13</sup> *D*<sup>7</sup><sup>23</sup> *D*<sup>8</sup><sup>23</sup> *D*∗ <sup>5</sup> <sup>12</sup> *S*2(ˆ) *D*<sup>5</sup><sup>23</sup> −*D*<sup>6</sup><sup>12</sup> 0 *D*<sup>6</sup><sup>23</sup> *D*<sup>7</sup><sup>13</sup>*ν*<sup>3</sup> *D*<sup>8</sup><sup>13</sup>*ν*<sup>∗</sup> 3 −*D*<sup>5</sup><sup>13</sup> *D*<sup>∗</sup> <sup>5</sup> <sup>23</sup> *S*3(ˆ) −*D*<sup>6</sup><sup>13</sup> −*D*<sup>6</sup><sup>23</sup> 0 −*D*<sup>7</sup><sup>12</sup>*ν*<sup>∗</sup> <sup>3</sup> −*D*<sup>8</sup><sup>12</sup>*ν*<sup>3</sup> 0 −*D*<sup>∗</sup> <sup>6</sup> <sup>12</sup> −*D*<sup>∗</sup> <sup>6</sup> <sup>13</sup> *S*1(ˆ) *D*<sup>∗</sup> <sup>5</sup> <sup>12</sup> −*D*<sup>5</sup><sup>13</sup> −*D*<sup>∗</sup> <sup>8</sup> <sup>23</sup> *D*<sup>∗</sup> <sup>7</sup> <sup>23</sup> *D*∗ <sup>6</sup> <sup>12</sup> 0 −*D*<sup>∗</sup> <sup>6</sup> <sup>23</sup> *D*<sup>5</sup><sup>12</sup> *S*2(ˆ) *D*<sup>∗</sup> <sup>5</sup> <sup>23</sup> −*D*<sup>∗</sup> <sup>8</sup> <sup>13</sup>*ν*<sup>3</sup> *D*<sup>∗</sup> <sup>7</sup> <sup>13</sup>*ν*<sup>∗</sup> 3 *D*∗ <sup>6</sup> <sup>13</sup> *D*<sup>∗</sup> <sup>6</sup> <sup>23</sup> 0 −*D*<sup>∗</sup> <sup>5</sup> <sup>13</sup> *D*<sup>5</sup><sup>23</sup> *S*3(ˆ) *D*<sup>∗</sup> <sup>8</sup> <sup>12</sup>*ν*<sup>∗</sup> <sup>3</sup> −*D*<sup>∗</sup> <sup>7</sup> <sup>12</sup>*ν*<sup>3</sup> *D*∗ <sup>7</sup> <sup>23</sup> *D*<sup>∗</sup> <sup>7</sup> <sup>13</sup>*ν*<sup>∗</sup> <sup>3</sup> −*D*<sup>∗</sup> <sup>7</sup> <sup>12</sup>*ν*<sup>3</sup> −*D*<sup>8</sup><sup>23</sup> −*D*<sup>8</sup><sup>13</sup>*ν*<sup>∗</sup> <sup>3</sup> *D*<sup>8</sup><sup>12</sup>*ν*<sup>3</sup> *D*4Trˆ 0 *D*∗ <sup>8</sup> <sup>23</sup> *D*<sup>∗</sup> <sup>8</sup> <sup>13</sup>*ν*<sup>3</sup> −*D*<sup>∗</sup> <sup>8</sup> <sup>12</sup>*ν*<sup>∗</sup> <sup>3</sup> *D*<sup>7</sup><sup>23</sup> *D*<sup>7</sup><sup>13</sup>*ν*<sup>3</sup> −*D*<sup>7</sup><sup>12</sup>*ν*<sup>∗</sup> <sup>3</sup> 0 *D*4Trˆ ⎞ ⎟⎟⎟⎟⎟⎠ , (A7)

where *D*<sup>1</sup> = 0.318 eV, *D*<sup>2</sup> = 0.211 eV, *D*<sup>3</sup> = 0.261 eV, *D*<sup>4</sup> = 0.270 eV are real and determine the shift of each pair of energy levels under [100] deformation. The splitting of levels under the compressive or tensile strain *e* in [110] direction also can be obtained analytically. Two doublets linearly shift with deformation *e*

$$\epsilon\_R = \left( 2D\_1 + D\_2 + D\_3 \pm \sqrt{(D\_2 - D\_3)^2 + 4(|D\_5|^2 + |D\_6|^2)} \right) \varepsilon / 4,\tag{A8}$$

from which |*D*5|<sup>2</sup> + |*D*6|<sup>2</sup> = 3.5 eV can be obtained. But another four levels shift nonlinearly

$$\varepsilon\_{R} = \left( (D\_{2} + D\_{3} + 2D\_{4})\varepsilon - 2\Delta \pm \sqrt{((D\_{2} + D\_{3} - 2D\_{4})\varepsilon + 2\Delta)^{2} + 4\varepsilon^{2}(|D\gamma|^{2} + |D\varsigma|^{2})} \right) / 4. \tag{A9}$$


#### **Appendix B. Absolute Deformation Potentials and Work Function of CoSi**

In order to obtain absolute shift of energy level *<sup>n</sup>* after deformation, it is necessary to have common reference energy in deformed and undeformed crystals or to determine the shift of the reference due to deformation. For example, energy can be measured

from macroscopic average of effective self-consistent potential *Ve*. Then absolute shift of energy level *<sup>n</sup>* due to deformation is equal to Δ*<sup>n</sup>* = ( (*d*) *<sup>n</sup>* <sup>−</sup> *<sup>V</sup>*(*d*) *<sup>e</sup>* ) <sup>−</sup> ( (*u*) *<sup>n</sup>* <sup>−</sup> *<sup>V</sup>*(*u*) *<sup>e</sup>* ) + <sup>Δ</sup>*Ve*, where Δ*Ve* is the reference energy offset due to strain and superscript *d*(*u*) corresponds to bulk calculation for deformed (undeformed) crystal.

The vacuum level can be used as a common reference energy, but it is not accessible in bulk DFT calculation. Hence we apply approach similar to that used for work function calculations. The superlattice configuration was considered with alternating layers of material and vacuum gaps. Average effective potential inside material layer *Ve* was calculated relative to its value inside vacuum gap (vacuum energy level). Then, the change of average effective potential Δ*Ve* due to deformation can be calculated as a deference between the values obtained from separate calculations for strained and unstrained layers.

Alternative approach was used in Ref. [39], where another superlattice method was proposed in order to obtain reference energy offset. The superlattice was formed from layers, extended or compressed along the direction of superlattice axis, and the layers were undeformed in the plane. In all-electron calculations, performed in Ref. [39], localized core levels, used as an energy reference, can be associated with each of the layers. The difference in their energy positions in the limit of thick layers allowed to obtain the reference energy offset due to deformation. Similar approach was used in Ref. [40], where pseudopotential calculations were used and, instead of core levels, macroscopic average effective potential in deformed *<sup>V</sup>*(*DL*) *<sup>e</sup>* and undeformed *<sup>V</sup>*(*UL*) *<sup>e</sup>* layers was used to determine the change of the energy reference due to deformation <sup>Δ</sup>*Ve* <sup>=</sup> *<sup>V</sup>*(*DL*) *<sup>e</sup>* <sup>−</sup> *<sup>V</sup>*(*UL*) *<sup>e</sup>* . We also used the latter approach and made similar calculations for superlattice of strained/unstrained layers of CoSi for [100], [110] and [111] directions. We checked the convergence of Δ*Ve* with respect to the layer thickness. The accuracy of 1–2 meV was reached for the layer thickness of 10*a*0.

Inside thick metallic layers, thicker then screening length, the difference (*Ve* − *<sup>F</sup>*) is determined only by its bulk properties, and the superlattice made of strained/unstrained layers should have common Fermi level *<sup>F</sup>*. Hence, the same Δ*Ve* can be obtained from bulk calculations for deformed and undeformed crystal <sup>Δ</sup>*Ve* = (*V*(*d*) *<sup>e</sup>* <sup>−</sup> (*d*) *<sup>F</sup>* ) <sup>−</sup> (*V*(*u*) *<sup>e</sup>* <sup>−</sup> (*u*) *<sup>F</sup>* ). Then, the absolute shift of energy level *<sup>n</sup>* due to deformation can be calculated as Δ*<sup>n</sup>* = ( (*d*) *<sup>n</sup>* <sup>−</sup> (*d*) *<sup>F</sup>* ) − ( (*u*) *<sup>n</sup>* <sup>−</sup> (*u*) *<sup>F</sup>* ).

All three considered approaches should give the same results for metallic material. Although CoSi is considered as semimetallic, the comparison gave the same results for Δ*Ve* to within 1–2 meV. In addition, we obtained work function for CoSi, equal to 4.55 eV which compares favourably with experimental values of 4.47–4.54 eV [41].

#### **References**


## *Article* **Faraday Rotation Due to Quantum Anomalous Hall Effect in Cr-Doped (Bi,Sb)2Te3**

**Alexey Shuvaev 1,\*, Lei Pan 2, Peng Zhang 2, Kang L. Wang <sup>2</sup> and Andrei Pimenov <sup>1</sup>**


**\*** Correspondence: alexey.shuvaev@tuwien.ac.at

**Abstract:** Quantum anomalous Hall effect (QAHE) represents a quantized version of the classical anomalous Hall effect. In the latter case the magnetization takes over the role of magnetic field and induces nonzero off-diagonal elements in the conductivity matrix. In magnetic topological insulators with the band inversion the QAHE can be reached due to quantized conduction channel at the sample edge if the Fermi energy is tuned into the surface magnetic gap. In the static regime the QAHE is seen as a zero-field step in the Hall resistivity. At optical frequencies this step is transformed into a quantized value of the polarization rotation approaching the fine structure constant *<sup>α</sup>* <sup>=</sup> *<sup>e</sup>*2/2*ε*0*hc* <sup>≈</sup> 1/137. However, due to material issues the steps reach the predicted values at millikelvin temperatures only. In this work we investigate the Faraday polarization rotation in thin films of Cr-doped topological insulator and in the sub-terahertz frequency range. Well defined polarization rotation steps can be observed in transmittance in Faraday geometry. At temperatures down to *T* = 1.85 K the value of the rotation reached about 20% of the fine structure constant and disappeared completely for *T* > 20 K.

**Keywords:** quantum anomalous Hall effect; Faraday rotation; topological insulators; terahertz spectroscopy

#### **1. Introduction**

Topological insulators [1,2] are materials with insulating bulk but revealing conducting surface states. These states possess different helicity thus making them symmetry protected against non-magnetic scattering processes. In two-dimensional (2D) systems and in external magnetic fields the quantized off-diagonal conductivity is observed proportional to an integer times the conductivity quantum *e*2/*h*. However, the application of external magnetic fields can be avoided in magnetically-doped topological insulators [3], see Refs. [4,5] for reviews. In that case the coupling between magnetic moments of the dopants must be strong enough to obtain a magnetically ordered state with finite static magnetization [6]. In addition, the Fermi level must be shifted into the surface magnetic gap thus leading to a single quantized conducting channel at the edge. QAHE in magneticallydoped topological insulators can be seen as a last station starting from the classical Hall resistivity via time reversal breaking and quantization effects [7].

The predicted QAHE has been first observed [8] in Cr-doped (Bi,Sb)2Te3 and at millikelvin temperatures. In these, as well as in several similar experiments [9,10], the temperature range at that well-quantized conduction states are observed is substantially lower than the transition temperature to a magnetically ordered state of typically 20 K. This discrepancy is normally attributed to the sample issues like presence of additional dissipative channels [4,5]. More recently, QAHE in several magnetic topological insulators [11–13] could be achieved at temperatures close to 1 K.

**Citation:** Shuvaev, A.; Pan, L.; Zhang, P.; Wang, K.L.; Pimenov, A. Faraday Rotation Due to Quantum Anomalous Hall Effect in Cr-Doped (Bi,Sb)2Te3. *Crystals* **2021**, *11*, 154. https://doi.org/10.3390/cryst11020154

Academic Editor: Artem Pronin Received: 24 December 2020 Accepted: 30 January 2021 Published: 3 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional clai-ms in published maps and institutio-nal affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In conducting materials a step in the Hall resistivity should normally lead to a step in optical properties. We note here an existing discrepancy in interpretation of the quantum Hall effect in statics and dynamics of the 2D systems [14]. In the former case the quantization is explained by counting the quantum conducting channels at the edge [2,15]. In the optical case, the edges are normally excluded due to contact-free technique. Therefore, to explain the quantization of the optical data, one have to return to the quantization of the bulk conduction [16,17]. In the dynamic regime, one of the typical observation in twodimensional systems is the polarization rotation of linearly-polarized light in transmission geometry, see Figure 1.

**Figure 1.** Scheme of the optical experiment. Linearly polarized incident light is transformed into the elliptical polarization after the sample, and is characterized by Faraday rotation angle *θ* and ellipticity *η*. The analyzer in front of the detector projects the ellipse either into the same direction as the incident beam leading to parallel transmittance *txx* or to perpendicular direction leading to the crossed transmittance *txy*. An external magnetic field is applied parallel to the propagation direction (Faraday geometry).

Full expressions for the polarization rotation include the influence of the substrate and are given elsewhere [18–20]. However, two important approximations substantially simplify the interpretation of the data. First, the influence of the substrate can be removed if the field-dependent experiment is done in the maximum of the Fabry-Pérot resonances of the substrate (see Figure 6 ). At such frequencies the expressions for the Faraday rotation reduce to the substrate-free result [18,21]. Second, in most cases the thin film approximation can be used assuming that the influence of the conducting film is small: (*σxxZ*0, *σxyZ*0)  1. Here *σxx* is the 2D diagonal conductivity, *σxy* is the 2D Hall conductivity and *Z*<sup>0</sup> ≈ 377 Ω is the impedance of free space.

In the following we assume that the incident radiation is linearly polarized with the *ac* electric field along the *x*-axis and is propagating along the *z*-axis. In the thin film approximation the transmittance amplitudes in the parallel *txx* and perpendicular *txy* channels are given by [18]

$$t\_{\rm xx} \approx 1 - \sigma\_{\rm xx} Z\_0 / 2 \approx 1 \qquad \text{and} \qquad t\_{\rm xy} \approx \sigma\_{\rm xy} Z\_0 / 2 \,, \tag{1}$$

respectively. In present experiments, *txx* and *txy* are measured putting the analyzer parallel and perpendicular to the polarization of the incident beam. The phase shift (or optical thickness) of both signals are obtained using the Mach-Zehnder interferometer arrangement, see Methods Section.

Further on, especially for the samples of the present work, the scattering time of the charge carriers is rather small, thus the frequency dependent terms *ωτ* in conductivities *σxx* and *σxy* can be neglected. This can be derived from the fact that we do not observe any signs of the cyclotron resonance in the present range of frequencies and magnetic fields. This indicates that the resonance terms in the Drude conductivity with typical width *ωτ* are negligible. In the same approximation the (Faraday) rotation angle *θ* can be written as [18–20]:

$$
\theta \approx \tan(2\theta)/2 \approx \theta \mathcal{R}(t\_{xy}/t\_{xx}) \approx t\_{xy} \approx \sigma\_{xy} Z\_0/2 \,. \tag{2}
$$

In the quantum regime we expect only a single conduction channel with *σxy* = *e*2/*h* leading to

$$
\theta \approx \sigma\_{xy} Z\_0 / 2 = \mathfrak{a} \,, \tag{3}
$$

where *<sup>α</sup>* = *<sup>e</sup>*2*Z*0/2*<sup>h</sup>* = *<sup>e</sup>*2/2*ε*0*hc* ≈ 7.3 × <sup>10</sup>−<sup>3</sup> rad is the fine structure constant.

#### **2. Results and Discussion**

Figure 2 shows typical magnetic field dependence of the transmittance in crossed polarizers geometry that is most sensitive to weak polarization rotations. In these experiments the amplitude of the signal corresponds to the transmittance amplitude |*txy*| and the optical thickness is related to the phase shift of the transmittance. The absolute values of the optical thickness are mainly determined by the thickness and refractive index of the substrate *ϕ* ≈ *nsds* (see Section 4). The magnetic field-induced changes can be attributed to the film properties that are basically determined by *σxy* in this geometry in agreement with Equation (1). The data reveal a clear step at zero magnetic fields with a hysteresis of about 0.09 T. Compared to the *dc* data shown in Figure 5 below, the transmittance is not affected by the contact resistivity and, therefore, provides more direct information on the sample conductivity.

After a calibration to absolute values, the complex polarization rotation angle *θ* + *iη* can be calculated either using the simplified Equations (1) and (2) or via the exact procedure [18–20].

**Figure 2.** Magnetic field scans of the transmittance in (Cr0.12Bi0.26Sb0.62)2Te3 film and in crossed polarizers geometry *txy* <sup>=</sup> <sup>|</sup>*txy*|*eiϕ*. The external field is applied parallel to the propagation direction (Faraday geometry, see Figure 1). The parameters of the experiment are given in the plot. Bottom panel: amplitude of the crossed signal. Top panel: relative optical thickness (phase shift) of the sample.

Complex polarization rotation angles at the lowest temperature of our experiments (*T* = 1.85 K) and at various frequencies are shown in Figure 3. We observe that in the frequency range of the present experiment the rotation angle is approximately a real number, as the ellipticity corresponds roughly to the noise level of the spectrometer. Similarly to the raw transmittance data in Figure 2, the Faraday rotation angle shows a clear step-like function across zero magnetic field. The inset in Figure 3 shows the absolute values of the rotation angle step at zero magnetic field and as a function of frequency.

**Figure 3.** Complex polarization rotation angle *θ* + *iη* in (Cr0.12Bi0.26Sb0.62)2Te3 at *T* = 1.85 K and at different frequencies. The frequencies were selected at the maxima of the Fabry-Pérot interferences to suppress the effect of the substrate in the spectra. Bottom panel: Faraday rotation angle *θ*, top panel: ellipticity *η*. The inset shows the absolute values of the rotation angle due to quantum anomalous Hall effect (QAHE) in the units of the fine structure constant at zero magnetic field and as a function of frequency. Straight dashed line is to guide the eye.

We conclude that in the frequency range of the present experiment the Faraday angle is roughly frequency independent at the value *θ* ≈ 0.2*α* and the variation of the data corresponds to the uncertainties of the experiment. In order to get more arguments on the absolute values of the step across the zero field, we investigated the temperature dependence of the Faraday rotation. These results are shown in Figure 4. As also seen in the frequency-dependent rotation angles, Figure 3, the ellipticity in our data is close to zero within the experimental uncertainties. The Faraday rotation, as shown in the bottom panel and in the inset to Figure 4, decreases with increasing temperature. In our experiments the step disappears around 20 K that agrees reasonably well with Curie temperature estimated [9] as *T*<sup>C</sup> ≈ 30 K. The Faraday step at our lowest temperatures is Δ*θ*(0) ≈ 1.3 mrad ≈ 0.18*α*. This value is substantially smaller than 1.0*α* expected within simple arguments. However, it is still possible that scattering processes suppressing Δ*θ* will freeze out at millikelvin temperatures. We conclude that at temperatures down to 1.85 K additional dissipative channels like residual carriers from bulk bands or scattering by impurities [4,5,22] are still present, which impedes the dissipationless character of the chiral states and suppresses the universally quantized values of the Hall resistance and of the Faraday rotation.

**Figure 4.** Complex polarization rotation angle *θ* + *iη* in (Cr0.12Bi0.26Sb0.62)2Te3 at *ν* = 203 GHz and at different temperatures. Bottom panel: Faraday rotation angle *θ*, top panel: ellipticity *η*. The inset shows the absolute values of the step of the rotation angle in the units of the fine structure constant at zero magnetic field and as a function of temperature.

As the transmission experiments are done in the Faraday geometry and the sample is magnetic, the magneto-optical Faraday effect [23] cannot be a priori neglected. In fact, in the small angle approximation, the rotation angles due to the off-diagonal conductivity *σxy* and due to the static magnetization *M*<sup>0</sup> are simply added. To estimate the value of the last effect, we use Equation (8) of the Methods Section. The value of the static magnetization in our samples *<sup>μ</sup>*0*M*<sup>0</sup> ≈ 0.9 × <sup>10</sup>−<sup>2</sup> T has been measured in Ref. [9]. It agrees well with an estimate assuming fully ordered moments of Cr3<sup>+</sup> ions. Putting the numbers into Equation (8) we finally get:

$$
\theta\_m \sim 10^{-8} \text{ rad} \ll \alpha \text{ .} \tag{4}
$$

We see that in most cases dealing with QAHE the classical Faraday effect can be neglected.

Finally, we compare the static and dynamic results in our sample. In agreement with Equations (1) and (2) direct correspondence between both properties may be expected. As discussed in the Methods Section, the resistivity of indium contacts was too high thus distorting the magnetic field dependencies of the diagonal and Hall resistivity. Reasonable step-like Hall resistivity data could be obtained for *T* ≥ 5 K only, see Figure 5. Although the Hall data were distorted, we still could estimate the steps across zero field (inset to Figure 5) and compare them with the dynamic data in Figures 3 and 4. We see that both steps disappear at temperatures close to 20 K. The absolute values of the Hall resistivity in Figure <sup>5</sup> correspond to <sup>Δ</sup>*Rxy* ∼ 0.5*h*/*e*<sup>2</sup> in our lowest temperature of 1.85 K that deviates from the values of Faraday rotation ∼0.2*α* observed in the transmittance data. We recall, however, that the resistivity was strongly affected by highly resistive contacts, although the rotation is measured by a contact-free technique.

**Figure 5.** Magnetoresistance data in (Cr0.12Bi0.26Sb0.62)2Te3 film at various temperatures. Top: diagonal resistivity; bottom: Hall resistivity. The inset shows the absolute values of the step in the Hall resistance at zero magnetic field and as a function of temperature.

#### **3. Conclusions**

In this work we investigated the polarization rotation of the sub-terahertz light in thin films of Cr-doped topological insulator (Cr0.12Bi0.26Sb0.62)2Te3. The optical data are compared to the step in the quantum Hall conductivity measured by static technique. Well defined polarization rotation steps can be observed in transmittance at different frequencies and temperatures. At the lowest temperature of *T* = 1.85 K the value of the rotation angle reached about 20% of the fine structure constant and disappeared completely for *T* > 20 K. We estimate that pure magnetic contribution to the Faraday rotation can be neglected in the present case.

#### **4. Materials and Methods**

Single-crystalline (Cr0.12Bi0.26Sb0.62)2Te3 films on insulating (111) GaAs substrates were grown by molecular beam epitaxy [9,10,24]. Both the Cr doping level (12%) and the (Bi/Sb) ratio (0.3/0.7) were optimized so that the Fermi level positions of the as-grown samples were close to the charge neutrality point. The growth was monitored by reflection high-energy electron diffraction and the films with a thickness of 6 quintuple layers (∼6 nm) were obtained. After the film growth, a 2 nm Al was evaporated to passivate the films. During the growth procedure the back side of the sample was fully covered with indium film that was nontransparent for the terahertz radiation. Therefore, prior to the optical experiments this film was removed by polishing. To measure the static resistivity, indium contacts were made at the corners of the hexagon-like sample at soldering temperature of 560 K. Unfortunately, this procedure did not provide good contacts, thus reasonable static Hall resistivity could be measured at *T* ≥ 5 K only (see Figure 5).

Terahertz transmittance experiments at frequencies 0.1 THz < *ν* < 1.0 THz were carried out in a quasioptical arrangement [18,25] which allows measurements of the amplitude and phase shift of the electromagnetic radiation in a geometry with controlled polarization. The spectrometer utilizes linearly polarized monochromatic radiation which is provided by backward-wave oscillators, and a He-cooled bolometer as a detector. The amplitude |*t*| and the phase shift *ϕ* of the radiation transmitted through the sample are measured by using the Mach-Zehnder interferometer setup. Static magnetic field up to 7 Tesla is applied to the sample using a split-coil superconducting magnet with mylar windows. The polarization state of the transmitted beam is determined by measuring the amplitude and phase shift of the radiation both with parallel and crossed polarizer and analyzer. This procedure provides the complex values of *txx* and *txy*, respectively (see Figure 1).

Figure 6 shows the transmittance spectra of the (Cr0.12Bi0.26Sb0.62)2Te3 thin film in the frequency range of the present experiment. Due to Fabry-Pérot resonances within the substrate a clear periodic modulation is seen in the spectra. The frequency positions of the maxima correspond to a resonance relation 2*nsdsν* = *m*, where *m* is an integer, *ds* = 0.478 mm is the sample thickness and *ns* = 3.02 is the refractive index of the substrate. As mentioned in the Introduction section, doing magnetic field-dependent experiments at the maxima of the resonances lead to the Faraday rotation angle that is close to that of the free-standing film, thus strongly simplifying the interpretation of the data. We stress, however, that exact expressions given in detail elsewhere [18,21] have been used to calculate the angle of the polarization rotation. The transmittance maxima in the frequency range 120–500 GHz are close to unity supporting the approximation of a weakly conducting sample. In the frequency range close to 1 THz the absolute values of the transmittance are by about 20 % less than unity. We attribute this effect to a slight non-parallel surfaces of the substrate that appeared after polishing of the backside of the sample. This effect is expected to produce an amplitude correction proportional to the ratio *δ* · *ns*/(*λ*/*D*). Here *δ* · *ns* is the deviation angle *δ* enhanced by the substrate refractive index *ns*, and *λ*/*D* is the diffraction angle as a ratio of the radiation wavelength and the sample aperture. At high frequencies the wavelength becomes smaller, thus enhancing the effect.

Finally , we estimate the value of the magnetooptical Faraday effect on the polarization rotation in terahertz experiments. In calculations below, we neglect the influence of the substrate as the measurements are done in the maxima of the Fabry-Pérot interferences. In addition, we assume isotropic electromagnetic susceptibilities and the normal incidence. Then, in a thin sample approximation, i.e. for *εd*/*λ*  1; *μ*±*d*/*λ*  1, the boundary conditions can be written in an extended manner that includes the sample as part of the surface [26,27]. The transmittance of eigenmodes for a magnetic thin film can then be written as

$$
\mu\_{\pm} \approx 1 - \frac{i\pi d}{\lambda} (\varepsilon + \mu\_{\pm}) \,. \tag{5}
$$

Here *d* is the sample thickness, *ε* is the permittivity, *μ*<sup>±</sup> is the permeability for two circular polarizations, and *λ* is the radiation wavelength. We recall that in the present geometry circularly polarized waves are the eigenmodes of the system and that Equation (5) is closely similar to Equation (1) and to the purely magnetic case in Ref. [28].

We apply now the definitions given in Equation (2) to obtain the magnetic part of the polarization rotation via 2*txx* = *t*<sup>+</sup> + *t*<sup>−</sup> and 2*itxy* = *t*<sup>+</sup> − *t*−:

$$
\theta\_m \approx \frac{i\pi d}{\lambda} (\chi\_+ - \chi\_-) \,. \tag{6}
$$

Here *χ*<sup>±</sup> = (*μ*<sup>±</sup> − 1) are magnetic susceptibilities. Similar expression for the polarization rotation including magnetoelectric susceptibilities has been obtained in Ref. [27]. For a ferromagnetic material *χ*<sup>±</sup> can be written as [28,29]:

$$\chi\_{\pm} = \frac{\gamma M\_0}{\omega\_0 \mp \omega + i \lg \omega},\tag{7}$$

where *ω*<sup>0</sup> = *γ*|*H* − *M*0| is the ferromagnetic resonance frequency in the Faraday geometry, *M*<sup>0</sup> is the static magnetization, *γ* is the gyromagnetic ratio, *g* is the Gilbert damping

parameter and *H* is the external magnetic field (here we avoid using usual notation *α* for the Gilbert damping).

The ferromagnetic resonance frequency can be estimated as *ω*<sup>0</sup> = *γ*|*H* − *M*0| ∼ *γH* ∼ 3 − 10 GHz for fields below *μ*0*H* ∼ 0.3 T. Therefore, the useful approximation in the present case is *ω*0, *gω ω* leading to a simple expression for the magnetic Faraday angle:

$$
\theta\_m = \gamma \mathcal{M}\_0 \frac{d}{2\mathcal{L}}\,. \tag{8}
$$

Finally, it should be noted that Equation (6) differs substantially from the expression *θ <sup>m</sup>* = *<sup>π</sup><sup>d</sup> <sup>λ</sup>* (*n*<sup>+</sup> <sup>−</sup> *<sup>n</sup>*−) used in the classical books by a factor of <sup>√</sup>*<sup>ε</sup>* <sup>∼</sup> 10 for Bi2Te3 [30]. The latter case is derived for a thick sample and it neglects the influence of the surfaces that are dominating in the thin-film geometry.

**Figure 6.** Transmittance spectrum of the sample used in this work in zero magnetic field and at temperature *T* = 1.85 K. This spectrum was measured in the parallel polarizers geometry, *txx*, and shows a series of Fabry-Pérot resonances due to reflections on the substrate surfaces.

**Author Contributions:** Terahertz experiments, A.S.; Sample preparation, L.P., P.Z. and K.L.W.; writing of the manuscript A.P.; Project administration, A.P. and A.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by Austrian Science Funds (Grants No. W-1243, No. P27098-N27, and No. I3456-N27). Open Access was funded by Austrian Science Funds (FWF).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

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

#### **References**


## *Article* **Isotropic Nature of the Metallic Kagome Ferromagnet Fe3Sn2 at High Temperatures**

**Rebecca L. Dally 1,\*, Daniel Phelan 2, Nicholas Bishop 3, Nirmal J. Ghimire 3,4 and Jeffrey W. Lynn <sup>1</sup>**


nbishop3@masonlive.gmu.edu (N.B.); nghimire@gmu.edu (N.J.G.)


**Abstract:** Anisotropy and competing exchange interactions have emerged as two central ingredients needed for centrosymmetric materials to exhibit topological spin textures. Fe3Sn2 is thought to have these ingredients as well, as it has recently been discovered to host room temperature skyrmionic bubbles with an accompanying topological Hall effect. We present small-angle inelastic neutron scattering measurements that unambiguously show that Fe3Sn2 is an isotropic ferromagnet below *TC* ≈ 660 K to at least 480 K—the lower temperature threshold of our experimental configuration. Fe3Sn2 is known to have competing magnetic exchange interactions, correlated electron behavior, weak magnetocrystalline anisotropy, and lattice (spatial) anisotropy; all of these features are thought to play a role in stabilizing skyrmions in centrosymmetric systems. Our results reveal that at the elevated temperatures measured, there is an absence of significant magnetocrystalline anisotropy and that the system behaves as a nearly ideal isotropic exchange interaction ferromagnet, with a spin stiffness *D*(*T* = 480 K) = 168 meV Å2, which extrapolates to a ground state spin stiffness *D*(*T* = 0 K) = 231 meV Å2.

**Keywords:** inelastic neutron scattering; topological materials; anomalous Hall effect; isotropic ferromagnet; kagome; frustrated magnetism; skyrmion; magnetization

#### **1. Introduction**

The two-dimensional kagome lattice lends itself to hosting a variety of phenomena depending on the chemical species occupying the network of corner-sharing triangles. For example, the tight-binding model for itinerant electrons leads to an electronic spectrum with a flat band and two Dirac crossings at the symmetry protected *K* and *K* corner points of the hexagonal Brillouin zone. Chemical tuning can drive the Fermi level to meet the Dirac points (a Dirac semimetal) to realize chiral massless charge carriers such as that in graphene [1,2]. The prediction of the flat band—on the extreme opposite from a Dirac band—is the result of destructive interference of Bloch waves from the lattice geometry. Consequently, this nontrivial flat band can exhibit interesting physics such as flat-band ferromagnetism and a finite Chern number. Experimentally, FeSn was shown to host both flat bands and Dirac fermions [3] due to the isolated Fe kagome layers rendering it a nearly perfect realization of 2D kagome physics. Fe3Sn2 is similar in structure, but features isolated breathing kagome bilayers, as shown in Figure 1a–c. Interestingly, the bilayers and breathing structure were still theorized to have a band structure with similar features. Instead of one Dirac crossing at each *K* and *K* point, there are two which are symmetric about each point [4], and the fermions are both spin-polarized due to the breaking of time-reversal symmetry and massive due to the opening of a gap from spin-orbit coupling. The combination of these effects gives rise to a non-zero Berry curvature which is consistent

**Citation:** Dally, R.L.; Phelan, D.; Bishop, N.; Ghimire, N.J.; Lynn, J.W. Isotropic nature of the metallic kagome ferromagnet Fe3Sn2 at high temperatures. *Crystals* **2021**, *11*, 307. https://doi.org/10.3390/cryst11030307

Academic Editor: Artem Pronin

Received: 7 March 2021 Accepted: 18 March 2021 Published: 20 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

with the quadratic relationship of the anomalous Hall resistivity with the longitudinal resistivity [5], implying the intrinsic Karplus and Luttinger mechanism [6] is responsible for the large anomalous Hall effect. It was also recently shown that nearly flat bands near the Fermi surface exist, which may contribute to the observed high-temperature ferromagnetism [7].

**Figure 1.** Crystal structure and characterization of Fe3Sn2. The crystallographic space group is *R*3¯*m* with reported lattice parameters *a* = *b* = 5.344 Å and *c* = 19.845 Å [8]. (**a**) View along the **a**-axis. The solid black line represents the unit cell and Fe atoms are shown as the smaller green circles and Sn atoms are shown as the larger blue circles. A single Fe site is offset from any high symmetry position (*x*/*a* = 0.4949, *y*/*b* = 0.5051, and *z*/*c* = 0.1134) leading to two different Fe-Fe bond lengths in the *ab*-plane (the so-called "breathing" kagome). Two Fe-Fe bond lengths are shown, where the shorter bond is in orange, and the longer bond is in green. (**b**) A Sn-only layer viewed along the **c**-axis, where the Sn atoms are arranged on a honeycomb lattice. (**c**) An Fe-Sn layer viewed along the **c**-axis, showing the breathing kagome lattice made up of Fe atoms. The axes labels for (**c**) are the same as in (**b**), and the parallelogram outlined by a solid black line for both panels represents the unit cell. (**d**) Magnetization measurements taken at 770 K (solid lines) and 600 K (dashed lines). The samples show no signs of coercivity as the sweep down in field (blue lines) coincides with the sweep up (orange lines) in field. (**e**) Zero-field cooled (ZFC) magnetic susceptibility measurement taken in a 0.1 T applied magnetic field. The derivative (right axis) clearly shows the ferromagnetic transition at *TC* ≈ 665 K. (**f**) Neutron powder diffraction data taken above the magnetic transition at 680 K. The data demonstrate the structure is consistent with that reported. The upper set of red tic marks denote Fe3Sn2 Bragg peak positions and the lower set denote Al Bragg peak positions coming from the sample canister. A few small impurity peaks were observed but not identified, and these are marked by the **\*** symbol.

Metallic kagome ferromagnets clearly exhibit elegant physics; however, they are elusive with only two reported: Co3Sn2S2, a semimetal [9], and the aforementioned Fe3Sn2. As alluded to thus far, the electronic structure has signatures of non-trivial topology, but recently, Fe3Sn2 has garnered growing attention for the discovery of topologically nontrivial spin textures. The observation of room temperature skyrmion bubbles [10] quickly led to reports of nanostructured skyrmionic devices [11–14] and studies of the associated properties such as the topological Hall effect [15–17] and skyrmion thermopower [18]. The space group of Fe3Sn2 is the centrosymmetric *R*3¯*m*, meaning the mechanism for skyrmion bubble formation is not due to the conventional breaking of crystalline inversion symmetry with Dzyaloshinshkii-Moriya interactions found in conventional B20 skyrmion systems. Instead, topological magnetic structures in centrosymmetric systems are due to the presence of anisotropy and/or frustration. The underlying source of each which is needed to stabilize skyrmions has become widely studied in recent years. One common model is the triangular lattice with frustrated Heisenberg antiferromagnetic exchange

interactions. [19] Although the magnetic frustration alone was shown to lead to a skyrmion phase, either lattice and/or spin anisotropy [20], particularly easy-axis anisotropy [21,22], was shown to be helpful in stabilizing the skyrmions.

From a frustrated magnetism perspective, Fe3Sn2 has been of interest for quite some time. Original neutron powder diffraction measurements indicated collinear ferromagnetic order below the onset of magnetism at *TC* ≈ 660 K with moments oriented along the *c*axis [23]. A spin-reorientation transition to the *ab*-plane starting below 250 K was identified and later measurements implied a slightly non-collinear structure was more likely starting below 300 K [24] and that the spin-reorientation transition was actually first-order in nature and occurs at ≈150 K [25,26]. The non-collinearity is thought to be due to frustrated magnetic exchange, and would also explain the large anomalous Hall effect [5,27] and possibly some of the temperature regimes where the topological Hall effect is observed if the scaler spin chirality is finite. Bulk magnetic measurements have also shown Fe3Sn2 to be an extremely soft ferromagnet at all temperatures with no coercivity, implying any easy-axis magnetic anisotropy must be very weak.

Here, we present our small-angle inelastic neutron scattering study of the magnetic excitations in Fe3Sn2 between 480 K and 660 K and unambiguously show that no significant spin wave gap is observed within experimental uncertainties between these temperatures. Below 480 K, the spin stiffness parameter, *D*(*T*), becomes too large and the spin wave full-width-at-half-maximum in energy, Γ(*q*), has narrowed to the point that the excitations move outside our measurement window. However, our results show that down to at least 480 K, Fe3Sn2 behaves as an ideal isotropic ferromagnet, and any onset of significant magnetic anisotropy that may contribute to the topological spin textures must develop below this point.

#### **2. Materials and Methods**

Polycrystalline samples of Fe3Sn2 were synthesized by solid state reaction. Stoichiometric amounts of Fe powder (Alfa Aesar 99+%) and Sn powder (Alfa Aesar 99.995%) were mixed and pelletized. The pellet was sealed in a fused silica ampoule under vacuum. The sealed ampoule was heated to 800 ◦C at the rate of 1 ◦C/hour and was kept at 800 ◦C for 1 week. After 1 week, the ampule at 800 ◦C was quenched into ice water. The pellet was reground, re-pelletized, and sealed into the fused silica ampoule under vacuum and was annealed at 800 ◦C for 1 week.

Magnetic measurements were performed on a piece of pressed pellet of Fe3Sn2 powder employing a Quantum Design MPMS3 magnetometer with an oven heater stick between 300 K and 756 K.

Neutron powder diffraction (NPD) measurements were taken using the triple-axis spectrometer, BT-7, at the NIST Center for Neutron Research [28]. A 17 g sample of polycrystalline Fe3Sn2 was sealed in a cylindrical aluminum canister, which was mounted inside a closed cycle refrigerator. Data were collected in two-axis mode using a position sensitive detector and wavelength of 2.359 Å. Söller collimators of 50 − 40 *R* were used before and after the sample, respectively (where *R* indicates radial), and pyrolytic graphite (PG) filters were employed both in the reactor beam and after the sample to suppress higher order wavelength contributions. Data were refined using the Rietveld method and the program, FullProf [29].

Inelastic neutron scattering data were also taken using BT-7 and the same 17 g sample as in NPD. Two different small-angle inelastic neutron scattering configurations were used in order to obtain data over a wider temperature range. For higher temperatures (630 K to 660 K), PG(002) monochromator crystals with vertical focusing and PG(002) analzyer crystals were used, and constant-Q scans were taken with a fixed incident energy of 13.7 meV. Söller collimators of 10 − 10 − 10 − 25 were used before and after the monochromator and before and after the analyzer, and the reactor beam PG filter was once again employed. The vertical resolution was measured using a graphite crystal and found to be 0.16 Å<sup>−</sup>1. For lower temperatures (480 K to 610 K), PG(004) monochromator crystals

with vertical focusing and PG(004) analzyer crystals were used, and constant-*Q* scans were taken with a fixed incident energy of 35 meV. A velocity selector in the reactor beam was employed to suppress higher and lower order wavelengths. The same collimations as the *Ei* = 13.7 meV experiment were used and the vertical resolution at this higher energy was found to be 0.24 Å−1. The same scans taken at high temperatures were also taken at much lower temperatures (300 K for the *Ei* = 13.7 meV experiment and 250 K for the *Ei* = 35 meV experiment), and these data were used for background subtraction.

In the small-angle inelastic neutron scattering configuration, spin waves are probed in the long-wavelength (i.e., small-*q*) limit, and the dispersion for a ferromagnet is

$$
\hbar\omega(q) = \Delta + D(T)q^2,\tag{1}
$$

where Δ is any anisotropy gap and *D*(*T*) is the spin wave stiffness which in mean field theory is proportional to the magnetization. The kinematic constraints for the scattering severely restrict the range of energy transfers accessible, so that the spin waves can only be observed if there is little to no anisotropy gap. The point in reciprocal space where the spin waves are being probed can be viewed by the schematic in Figure 2a. Here, a parabolic dispersion about *Q* = 0 and energy transfer, *E* = 0, is shown. About this point, the dispersion of an isotropic ferromagnet powder sample is identical to that of a single crystal. Similar experiments on amorphous alloys [30] and powder samples of manganites [31] have been widely used to establish their isotropic nature.

**Figure 2.** (**a**) A three-dimensional schematic of an isotropic parabolic spin wave dispersion near *Q* = 0 and energy transfer, *E* = 0. The dispersion is shown as the orange surface, and dashed blue lines show the direction of constant-*Q* scans cutting through the dispersion surface along *E*. (**b**) A two-dimensional schematic demonstrating how the experiment captures the intensity from the spin wave excitations. The orange solid line represents the dispersion, *h*¯ *ω*(*q*) = Δ + *D*(*T*)*q*2, and the surrounding blue surface represents the full-width-at-half-maximum of the spread in energy of the dispersion, Γ(*q*), due to thermally induced magnon-magnon interactions. Dashed blue lines are examples of constant-*Q* scans made in the experiment. Overlayed on these lines are the instrumental resolution ellipses, *R*(*Q*, *E*), along *Q* = 0.07 Å−<sup>1</sup> and 0.11 Å<sup>−</sup>1. The left panel uses the refined values for the dispersion from the actual data at *T* = 580 K in the *Ei* = 35 meV experiment. The right panel uses the same parameters, but increased the gap to be 0.5 meV in order to demonstrate the sensitivity of the technique to the size of the gap. Here, the scans performed during the experiment wouldn't be able to reach the signal of the spin waves. (**c**) The actual data at *T* = 580 K in the *Ei* = 35 meV experiment. The solid orange lines are the refined fits to the data, shown as blue circles. The three constant-*Q* scans are vertically offset from one another for clarity.

One advantage of studying the spin waves in the long-wavelength limit (i.e., about *Q* = 0) is that the instrumental resolution is focused on both the energy gain and energy loss side, unlike with a Bragg point where there is a focused and de-focused side. More details on the resolution function in the small-angle limit can be found in Ref. [32] The intensity detected in a neutron scattering experiment represents a convolution of the instrumental resolution, *R*(*Q*, *E*), and the scattering function, *S*(*q*, *h*¯ *ω*), making it necessary to include the convolution when analyzing the data. Using the Cooper-Nathans approximation for the resolution, we used the program ResLib [33] to fit each set of data (where a set of data consists of all the constant-*Q* scans at a single temperature) to the scattering function,

$$S(q, \hbar \omega) \propto F(q, \hbar \omega) \frac{\hbar \omega}{1 - e^{-\hbar \omega/k\_B T}}\tag{2}$$

where *F*(*q*, *h*¯ *ω*) is the spectral weight function, *kB* is the Boltzmann constant, and *T* is the temperature. At finite temperatures, magnon-magnon interactions lead to damping effects in energy for the spin waves. For Heisenberg ferromagnets below *TC*, and in the energy regime ¯*hω kBT*, the excitation width in energy has been calculated [34] as

$$\Gamma(q) \propto q^4 T^2 \left\{ \frac{1}{6} \ln^2 \left( \frac{k\_B T}{\hbar \omega} \right) + \frac{5}{9} \ln \left( \frac{k\_B T}{\hbar \omega} \right) - 0.05 \right\}.\tag{3}$$

The shape of the broadening in energy is approximated using a Lorentzian function as the spectral weight function, *F*(*q*, *h*¯ *ω*), centered about *h*¯ *ω*(*q*), with Γ(*q*) as the full-widthat-half maximum.

#### **3. Results**

#### *3.1. Characterization*

Figure 1d shows the magnetization both above (770 K) and below (600 K) the ferromagnetic transition temperature. No coercivity was observed for either temperature, meaning Fe3Sn2 is a soft ferromagnet. Figure 1e shows the magnetic susceptibility as a function of temperature at an applied magnetic field of 0.1 T. The derivative of the susceptibility with respect to temperature shows the ferromagnetic transition to be ≈665 K, which is consistent with previous reports that show the Curie temperature to vary anywhere between 640 K and 660 K [5,16,24].

The observed NPD profile and Rietveld refined fit are shown in Figure 1f. The data were taken at 680 K and confirm the structure to be Fe3Sn2 with refined lattice parameters of *a* = *b* = 5.3787 ± 0.0004 Å and *c* = 19.863 ± 0.002 Å. The refined atomic positions for Fe are *x*/*a* = 0.4940 ± 0.0004, *y*/*b* = 0.5060 ± 0.0004, and *z*/*c* = 0.1132 ± 0002. The Sn positions are *z*/*c* = 0.1039 ± 0.0007 and *z*/*c* = 0.3318 ± 0.0007 for the Sn1 and Sn2 sites, respectively. The isotropic thermal parameters (*B*) were refined to 0.8 ± 0.1 Å2, 4.1 ± 0.5 Å2, and 2.6 ± 0.4 Å2 for the Fe, Sn1, and Sn2 sites, respectively. There were three small impurity peaks in the pattern that were unable to be identified. They are marked with an **\*** in Figure 1f.

#### *3.2. Inelastic Neutron Scattering*

We first demonstrate the sensitivity of the small-angle inelastic scattering configuration to the size of the gap in order to discern between isotropic and anisotropic ferromagnets. A schematic of a dispersion following Equation (1) near **Q** = 0 is shown in Figure 2a. The neutron scattering plane is defined by two arbitrary orthogonal vectors, **q***<sup>x</sup>* and **q***y*, and constant-*Q* cuts are shown as blue dashed lines to show how the experimental scans can cut through the dispersion along energy, *E*.

Each temperature set of constant-*Q* scans was fit globally to obtain the spin wave parameters, and the parameters for *T* = 580 K were used to create Figure 2b. The spin stiffness parameter was found to be *<sup>D</sup>*(*<sup>T</sup>* = 580 K) = <sup>135</sup> ± 3 meV Å<sup>2</sup> and the gap, Δ = 0.09 ± 0.02 meV, where the uncertainties throughout represent one standard deviation due to statistical counting. These parameters were used to create the solid orange line

representing *h*¯ *ω*(*q*). We note that even for an ideal isotropic spin system a small dipolar gap is expected due to ferromagnetic magnetization. The full-width-at-half-maximum of the spin waves in energy, Γ(*q*), is shown as the shaded blue region following Equation (3), and the instrumental resolution, *R*(*q*, *E*), is shown as black ellipses. The maxima and minima of the ellipses along energy represent the allowed scan region which satisfies the required conservation of momentum and energy represented by the scattering triangle (i.e., scanning farther in energy is not possible). Two representative constant-*Q* scans are shown as dashed blue lines at 0.07 Å−<sup>1</sup> and 0.11 Å<sup>−</sup>1, and with the small gap of 0.09 meV, the center of the instrumental resolution ellipses for both scans is able to pass over the peak of the dispersion on the energy gain and loss sides. This is not possible if the gap is increased to 0.5 meV, as shown in the right panel of Figure 2b (all other parameters from the *T* = 580 K fit were fixed). The actual data from *T* = 580 K are shown in Figure 2c as blue circles and the fits are shown as solid orange lines.

The spin stiffness parameter, *D*(*T*), was extracted from the fits for each temperature and is shown in Figure 3. The dashed line is a power law fit to the data: *D*(*T*) = *D*<sup>0</sup> - *TC*−*T TC ν*−*<sup>β</sup>* , where *<sup>D</sup>*<sup>0</sup> = <sup>271</sup> ± 9 meV Å2, *TC* = 662.4 ± 0.8 K, and *<sup>ν</sup>* − *<sup>β</sup>* = 0.34 ± 0.02. The solid line is a fit to the Dyson formalism of two spin-wave interactions in a Heisenberg ferromagnet, where *D*(*T*) = *D*<sup>0</sup> 1 − *A* - *kBT* 4*πD*<sup>0</sup> 5/2 *ζ* 5 2 , *ζ* 5 2 is the Riemann integral and *A* is a constant proportional to the interaction range [35]. The *T*5/2 temperature dependence is not valid near the critical regime, and only the lowest four temperatures were used in the fit, resulting in *<sup>D</sup>*<sup>0</sup> = <sup>231</sup> ± 7 meV Å2.

The gap was not found to have any meaningful temperature dependence, ranging between 0.06 meV and 0.09 meV, and was the same within plus or minus one standard deviation. It should also be noted that the instrumental resolution in energy for the scans taken is on the order of these values (see Figure 2b), meaning the exact fitted value for the gap is not well-defined. For example, in the *Ei* = 13.7 meV experiment, the resolution in energy at *Q* = 0.07 Å−<sup>1</sup> and *E* = 0 meV is 0.29 meV, and at *E* = 0.6 meV the resolution is 0.07 meV.

**Figure 3.** The temperature dependence of the spin wave stiffness parameter, *D*(*T*). Data from both the *Ei* = 13.7 meV and *Ei* = 35 meV experiments were included in the power law fit, *D*(*T*) = *D*<sup>0</sup> - *TC*−*T TC ν*−*<sup>β</sup>* , shown as the dashed line. Only the lowest four temperatures were used in the Dyson fit, *D*(*T*) = *D*<sup>0</sup> 1 − *A* - *kBT* 4*πD*<sup>0</sup> 5/2 *ζ* -5 2 , shown as the solid line.

#### **4. Discussion**

The temperature renormalization of the spin stiffness for Heisenberg ferromagnets is expected to follow a power law on approach to *TC* with the critical exponents *ν* − *β* = 0.34 [36], which is the exponent found in this study, further showing that magnetically, Fe3Sn2 is a

typical exchange ferromagnet at elevated temperatures. In fact, the *ν* − *β* exponent from the power law fit is strikingly similar to those for elemental Fe [37] (*ν* − *β* = 0.37) and Ni [38] (*ν* − *β* = 0.39), in addition to the amorphous iron magnets already mentioned in Ref. [30], using the same technique. We remark that the *TC* of 631 K for Ni is comparable to Fe3Sn2, but the spin stiffness of *D* = 550 meV Å<sup>2</sup> in the ground state of this itinerant magnet is much larger [39]. However, the extrapolation of *D*(*T*) to *T* = 0 K for Fe3Sn2 using the power law fit is likely over estimated due to its validity only near *TC*. The value of *<sup>D</sup>*(*<sup>T</sup>* = 0 K) = <sup>231</sup> ± 7 meV Å<sup>2</sup> using the Dyson formalism is a better estimate of the ground state spin stiffness, although the limited temperature range accessible in this study still results in a large extrapolation window.

All ferromagnets have a gap due to the magnetic dipole-dipole interaction between different atoms [40]. This gap is typically small and often out of the range of resolution for inelastic neutron scattering experiments. The dipole-dipole interaction or resolution/instrumental alignment effects are both probable reasons for the observation of a small gap in this study (on the order of 0.06 meV to 0.09 meV), which is quite small compared to the exchange energy rendering Fe3Sn2 an isotropic ferromagnet to an excellent approximation. However, the magnetic anisotropy energy due to magnetocrystalline anisotropy was recently calculated to be close to our gap value, at 0.037 meV per Fe atom for the ground state when the easy-axis and spins are oriented within the kagome plane [41]. The ground state for which this value was calculated, though, is in a different temperature regime and spin configuration than that of the present experiment, so it is unclear whether the gap observed is solely due to magnetic anisotropy energy coming from dipolar interactions and/or magnetocrystalline effects.

We now discuss the meaning surrounding the term "anisotropy" in our discussion. Previous studies have cited the uniaxial anisotropy in Fe3Sn2 as one of the necessary ingredients for the formation of the topologically protected skyrmionic bubbles [10], and many of the centrosymmetric skyrmion systems discovered thus far are well-known to be a result of competition between frustrated magnetic exchange and spin anisotropy [21,22]. Unsurprisingly, measurements of the anisotropy energy density, *Ku*, have therefore been published [11,14,15] and show that the onset of a magnetic anisotropy precedes the temperatures at which the skyrmion bubbles are found. However, anisotropy can range from preferred orientation of a spin—which all ordered crystalline magnets have—to an appreciable energy required to pull spins away from a preferred direction. As a soft ferromagnet, Fe3Sn2 falls into the former category and can be considered an isotropic ferromagnet in accordance with our results. This is in contrast to the large anisotropies required for permanent magnet devices for magnetostatic energy storage [42]. In fact, one of the appealing properties of skyrmion-based devices may be that the weak anisotropy requirements open up the field for potential skyrmion candidate materials, especially when considering inducing small anisotropies into materials via doping is quite common.

An example of a system that internally tunes its anisotropy is Nd2Fe14B, a hard uniaxial ferromagnet used in permanent magnet applications but also exhibits a spinreorientation transition such as that in Fe3Sn2. It was found that the rotating spins act to tune the overall anisotropy in the system [43], although in contrast with Fe3Sn2, the anisotropy is due to the lanthanide crystal field effect. It is also instructive to recall that spin anisotropy is not required for skyrmion formation in inversion symmetric systems [19,20], although these theories have not been specifically applied yet to Fe3Sn2. Another metallic breathing kagome lattice to host a skyrmion spin texture is Gd3Ru4Al12 [44]. In contrast to Fe3Sn2, the ordered magnetic state is antiferromagnetic and a weak anisotropy is of the easy-plane type. Future work on either of these kagome materials to directly probe the anisotropy gap in proximity to the skyrmion phases would be of interest to explore the role of the anisotropy versus magnetic exchange frustration. Inelastic neutron scattering can be used to achieve this below the temperatures accessible in the work presented here but would require a comparable mass of co-aligned single crystals and sub-meV instrumental resolution in a wide-angle scattering experiment.

**Author Contributions:** Material synthesis, N.J.G. and N.B.; Magnetization measurements, D.P.; neutron experiment design, J.W.L. and R.L.D.; neutron experiment analysis, J.W.L. and R.L.D.; writing—original draft preparation, R.L.D.; writing—review and editing, J.W.L., R.L.D., D.P., and N.J.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** Synthesis and characterization work (N.J.G.) were supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Science and Engineering Division. Work in the Materials Science Division at Argonne National Laboratory was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Science and Engineering Division.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data is available upon request to the corresponding author.

**Conflicts of Interest:** The identification of any commercial product or trade name does not imply endorsement or recommendation by the National Institute of Standards and Technology. The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


## *Article* **Fractional Power-Law Intraband Optical Conductivity in the Low-Dimensional Dirac Material CaMnBi2**

**M. B. Schilling 1, C. X. Wang 2, Y. G. Shi 2, R. K. Kremer 3, M. Dressel <sup>1</sup> and A. V. Pronin 1,\***


**Abstract:** We studied the broadband optical conductivity of CaMnBi2, a material with two-dimensional Dirac electronic bands, and found that both components of the intraband conductivity follow a universal power law as a function of frequency at low temperatures. This conductivity scaling differs from the Drude(-like) behavior, generally expected for free carriers, but matches the predictions for the intraband response of an electronic system in a quantum critical region. Since no other indications of quantum criticality are reported for CaMnBi2 so far, the cause of the observed unusual scaling remains an open question.

**Keywords:** Dirac materials; optical-conductivity scaling; topological semimetals

#### **1. Introduction**

CaMnBi2 and its sister compound SrMnBi2 are some of the first materials in which bulk electronic bands with Dirac-like dispersion were experimentally confirmed [1,2]. Both materials are arranged in layers with square nets of Bi atoms (space group P4/*nmm*). The materials are believed to possess two-dimensional Dirac bands that are anisotropic and slightly gapped due to spin-orbit coupling [1–5]. The materials have antiferromagnetic in-plane ordering of Mn ions with Néel temperatures between 270 and 290 K [1,3,6,7]. In CaMnBi2, another transition at Ts ≈ 50 K was detected by various experimental techniques including transport [2,3,8,9], magnetoresistance [2], susceptibility [3], thermopower [8], and optical [10,11] measurements. The signatures of Ts are often tiny and not always resolved in DC transport [6]. No indications of a phase transition were detected in specific-heat [3] and neutron measurements [7]. The anomaly at Ts was first tentatively attributed to either weak ferromagnetic order [2] or spin canting [3]. Based on optical and magnetic torque measurements in combination with band-structure calculations, Yang et al. [11] recently concluded on a spin-canting-induced band reconstruction at Ts, therefore clarifying the nature of this transition. In this paper, we report on the optical conductivity measurements in CaMnBi2. Below Ts, we found an unusual scaling of its intraband conductivity. This scaling was previously attributed to manifestations of quantum criticality. Hence, it might be an indication of quantum criticality in CaMnBi2, although other explanations cannot be excluded.

#### **2. Materials and Methods**

**Sample growth and characterization:** single crystals of CaMnBi2 were grown using a self-flux method similar to that described previously [6]. Elementary Ca (99.99%), Mn (99.9%), and Bi (99.99%) were mixed in the molar ratio Ca:Mn:Bi = 1:1:8 and put into an alumina tube before sealing it in a quartz tube. The mixture was heated up to 800 ◦C during 10 h, kept at this temperature for 5 h, then slowly cooled down to 450 ◦C at a

**Citation:** Schilling, M.B.; Wang, C.X.; Shi, Y.G.; Kremer, R.K.; Dressel, M.; Pronin, A.V. Fractional Power-Law Intraband Optical Conductivity in the Low-Dimensional Dirac Material CaMnBi2. *Crystals* **2021**, *11*, 428. https://doi.org/10.3390/ cryst11040428

Academic Editor: Andreas Hermann

Received: 30 March 2021 Accepted: 13 April 2021 Published: 16 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

rate of 3 ◦C/h. The excess Bi flux was decanted at this temperature in a centrifuge. As CaMnBi2 is somewhat air-sensitive, its handling was carried out in an inert gas atmosphere. The obtained samples were carefully characterized by X-ray, transport, magnetic, and specific-heat measurements as described in the Appendix A.

**Optical measurements:** the near-normal-incidence optical reflectivity R(ν) was measured from a large (roughly 2 by 3 mm) (001) surface of a CaMnBi2 crystal at a number of temperatures between 10 to 300 K over a broad frequency range from ν = 50 to 22,000 cm−<sup>1</sup> (≈6 meV–2.75 eV) using two Fourier-transform spectrometers (Bruker IFS 113v and Bruker Vertex 80v equipped with a Hyperion IR microscope (all three devices are from Bruker Corporation, Billerica, MA, USA)). At low frequencies, an in situ gold evaporation technique was utilized for reference measurements. For frequencies above 1000 cm<sup>−</sup>1, gold and protected silver mirrors served as references. The complex optical conductivity, σ(ν) = σ1(ν)+iσ2(ν), was obtained using Kramers–Kronig transformations. The highfrequency range was extended by involving the X-ray atomic scattering functions for high-frequency extrapolations [12]. The results of the four-point DC resistivity measurements were used for the low-frequency extrapolations. To avoid possible surface oxidation, all measurements were performed on freshly cleaved surfaces.

#### **3. Results**

The results of our optical experiments are shown in Figures 1 and 2. All measurements were obtained on (001) planes of CaMnBi2 (in-plane response). Let us note here that, although the in-plane Dirac bands of CaMnBi2 are known to be highly anisotropic [4,5], this band anisotropy is not expected to be seen in the linear optical response because of the fourfold in-plane symmetry. Indeed, our polarization-dependent reflectivity measurements do not reveal any optical anisotropy. Hence, we discuss the measurements performed with unpolarized light throughout the paper.

**Figure 1.** Frequency-dependent in-plane reflectivity R(ν) (panels (**a**,**b**)) and the real part of the optical conductivity σ1(ν) (panel (**c**)) of CaMnBi2 for select temperatures between 10 and 300 K. The development of a dip in R(ν) and a bump in <sup>σ</sup>1(ν) at around 1500 cm−<sup>1</sup> (≈200 meV) is clearly seen at T < 50 K and marked with the arrow in panel (**c**). Note the change in frequency scale at 8000 cm<sup>−</sup>1.

**Figure 2.** Scaling of the intraband complex optical conductivity in CaMnBi2. The real (panel (**a**)) and imaginary (panel (**b**)) parts of the optical conductivity follow the Drude behavior for T > Ts ~ 50 K and scale as 1/ω0.5 at low temperatures. The conductivity angle ϕ (panel (**c**)) is almost constant at these temperatures, while it increases quasi-linearly in the Drude case.

The obtained broadband optical spectra as functions of frequency are shown in Figure 1. The interband part of the optical conductivity has been analyzed in Reference [11]; our experimental findings are in full agreement with these results. The goal of our paper is to analyze the intraband low-energy response.

At low frequencies and for T > 50 K, the reflectivity R(ν) [ν = ω/(2π)] and both parts of the complex conductivity are dominated by intraband electronic transitions and are typically metallic: R(ν) approaches unity as ν diminishes, while σ1(ν) and σ2(ν) reveal typical Drude behavior, as can be seen best from Figure 2a,b. The conductivity angle, ϕ = arctan(σ2/σ1), is frequency dependent and follows the Drude model; see panel (c).

The low-energy interband transitions within the Dirac bands, which are known to provide a power-law contribution to low-frequency σ1(ν) [13–17], are not seen in our measurements. The low-frequency response of CaMnBi2 is completely dominated by free carriers. This situation is, in fact, rather typical for different Dirac systems, in which the Fermi level is situated far from the band crossings and/or the free-carrier contributions from non-Dirac bands are significant [18–23].

At the spin-canting temperature Ts ≈ 50 K, dramatic changes occur in the optical spectra. Apart from the formation of a low-frequency mode at approximately 200 meV (Figure 2a,b) and a corresponding dip in reflectivity (Figure 1a) that were reported previously [10,11], the intraband absorption also drastically changes its shape. The single Drude term is unable to describe the low-energy spectra. (Let us note that the two-Drude approach used in Reference [11] is able to provide only a very rough description of the experimental σ1(ω); see Figure 1b of Reference [11].) Instead, one can see that both components of the

optical conductivity follow a power law. This power-law behavior is most apparent at low temperatures: below 30 K. As evidenced from our fit, σ1(ω) and σ2(ω) both depend on frequency as 1/ω0.5; this power law is shown as black solid lines in Figure 2a,b. This behavior of conductivity strongly differs from the conventional Drude response. In particular, changes in the reactive part (σ2(ω)) are significant: instead of showing a broad maximum (which corresponds to the scattering rate in the Drude case), σ2(ω) is now monotonic in frequency. Furthermore, the pre-factors in the frequency dependencies of σ1(ω) and σ2(ω) are identical, i.e., σ<sup>1</sup> and σ<sup>2</sup> are equal at a given frequency. This provides that ϕ = π/4.

#### **4. Discussion**

The strong non-Drude intraband response of CaMnBi2 is rather surprising. Free electrons are generally expected to follow a Drude(-like) conductivity ansatz [24]. One can notice that the conductivity scaling observed here was widely discussed in the past in relation to quantum phase transitions (QPTs). The presence of a QPT leads to universal power-law scaling behaviors of the response functions [25]. In particular, the frequencydependent complex conductivity should follow such a behavior. For the frequency region where kBT < *h*¯ ω, the conductivity can be a universal function of frequency [25]. In our case, this inequality is fulfilled. Van der Marel et al. [26] argued that scale invariance, causality, and time reversal symmetry require that, for a quantum critical system, the complex conductivity in this frequency region follows:

$$\sigma(\omega) = |\sigma(\omega)| \mathbf{e}^{i\rho(\omega)} = \mathbb{C}\omega^{\gamma - 2} \mathbf{e}^{i\pi(1 - \gamma/2)},\tag{1}$$

where γ is a critical exponent and C is a constant. This ansatz implies that both σ1(ω) and σ2(ω) depend on frequency as ωγ−<sup>2</sup> and the phase ϕ is frequency independent and set by the same exponent γ. If we apply Equation (1) to the recorded spectra of σ<sup>1</sup> and σ2, we find the critical exponent γ to be 3/2. According to the scaling analysis, this value of the critical exponent should provide a frequency-independent value for the conductivity angle, ϕ = π/4. As noticed above, this result indeed follows from our data.

Interestingly, the same optical conductivity scaling (with γ = 3/2) was theoretically elaborated by Ioffe and Millis [27] in relation to a possible QPT in the superconducting cuprates. Van der Marel [28] suggested a generalized form of this relation (with the critical exponent not fixed at 3/2) that merges with a proposition of Anderson [29] for a one-dimensional Luttinger liquid in the collision-less limit.

It should be noted that the Dirac bands in CaMnBi2 are two-dimensional (due to the planar net of Bi atoms) and that they possess a very high anisotropy within this plane [4,5]. In fact, the electronic band structure can be viewed as a gapped dispersive nodal line in two dimensions, with the Fermi velocity in one direction being much smaller than the Fermi velocity in another one. It is known that the presence of a nodal line effectively reduces the dimensionality of electronic transport [16,17,30]. This reduction can possibly occur in CaMnBi2, leading to a quasi-one-dimensional situation and, eventually, to the realization of a quantum critical state.

Certainly, the nature of the observed scaling has to be clarified in further theoretical and experimental studies; the observed scaling is not necessarily related to quantum criticality. The goal of this paper is to report this unusual conductivity behavior and to suggest a possible explanation. Still, one can note that the possible quantum criticality in CaMnBi2 (if it is confirmed) should likely be related to the magnetism in this system and, particularly, to the spin-canting transition at Ts. This conclusion follows from the fact that the observed conductivity scaling appears only below this temperature.

#### **5. Conclusions**

In summary, we performed broadband optical conductivity measurements of CaMnBi2 —a highly anisotropic material with two-dimensional nets of Bi atoms and anisotropic Dirac bands. We detected the formation of a finite-frequency absorption mode at T<Ts = 50 K, which is in agreement with previous studies. Most importantly, the optical response of

itinerant electrons at these temperatures is not of the Drude type. Instead, it follows a fractional power-law behavior, σ1(ω)~σ2(ω)~ω−0.5, that is similar to the behavior proposed for quantum critical systems with the critical exponent γ *=* 3/2. These findings might indicate that CaMnBi2 is in the vicinity of a quantum phase transition. More input from the theory side and further experiments are necessary to confirm this proposition.

**Author Contributions:** Conceptualization, A.V.P.; sample growth, C.X.W. and Y.G.S.; measurements, M.B.S. and R.K.K.; data analysis, M.B.S. and A.V.P.; writing—original draft preparation, M.B.S.; writing—review and editing, M.D. and A.V.P.; funding acquisition, Y.G.S., M.D. and A.V.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partly funded by the Deutsche Forschungsgemeinschaft (DFG) via grant No. DR228/51-3 and by K. C. Wong Education Foundation (GJTD-2018-01).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data sets generated and analyzed during the current study are available from the corresponding author upon request.

**Acknowledgments:** We thank Vladimir Juriˇci´c and Run Yang for useful discussions and Gabriele Untereiner for technical support.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **Appendix A**

**Electronic transport and magnetic susceptibility:** temperature-dependent DC resistivity, ρ(T), was measured in a custom-made setup at temperatures down to 2 K. To highlight the discussed feature in resistivity near Ts, we plot ρ(T) on an enlarged scale in Figure A1a. To further confirm the Ts feature, we performed temperature-dependent measurements of magnetic susceptibility, χ(T), at 2 T down to 10 K with a commercial setup (MPMS, Quantum Design Inc., San Diego, CA, USA) based on a superconducting quantum interference device. The results of these measurements are shown in Figure A1b. To clearly observe the anomaly at Ts, we performed a polynomial fit of the χ(T) curve at T>Ts and subtracted the fit from the experimental data. The value of Δχ(T) obtained in this way is shown in Figure A1c. For T << Ts, Δχ(T) remains almost constant.

**Heat-capacity measurements:** heat capacity was measured as a function of temperature, employing the relaxation method (PPMS, Quantum Design Inc., San Diego, CA, USA). The sample was attached with Apiezon N vacuum grease to the sapphire platform, and the heat capacity of this platform (including the vacuum grease) was measured in advance and then subtracted from the total heat capacity. Figure A1d displays our results of the specific-heat measurements. The data reveal a small λ-type anomaly at 290 K, which is associated with the antiferromagnetic ordering of Mn atoms. The Néel temperature is in good agreement with that reported by Guo et al. [6] from neutron powder diffraction data and is slightly higher than the findings in other reports [3,7]. The low-temperature Cp/T data can be well fitted with a power law, Cp/T = γ + βT2 + δT<sup>4</sup> + εT6, with the Sommerfeld term γ = 6.78 mJ/(molK2). The higher powers of this polynomial represent the lattice and magnon contributions. They are very small and amount to β = 0.00126(1) J/(molK4), <sup>δ</sup> = 2.04(3) × <sup>10</sup>−<sup>5</sup> J/(molK6), and <sup>ε</sup> <sup>=</sup> −9.3(2) × <sup>10</sup>−<sup>8</sup> J/(molK8). Most importantly, Cp/T is featureless at Ts and in agreement with a previous report [3].

**Figure A1.** Characterization measurements of CaMnBi2. DC resistivity ρ(T) (**a**), magnetic susceptibility χ(T) (**b**), and the deviation of χ(T) from its high-temperature behavior (**c**). The data in panels (**a**–**c**) are shown at the temperatures near Ts to highlight the presence of a transition. Molar specific heat (**d**). The main frame shows a Sommerfeld plot (blue circles) with a fit (red line) to a polynomial in T<sup>2</sup> (see text). The insets provide an overview on a linear temperature scale and a magnification around the 290 K anomaly. The unit cell volume is shown as a function of temperature in panel (**e**). The red solid line represents a fit of the experimental data to Equation (A1) with the parameters given in the text. The inset displays a section of the diffraction pattern (λ = 0.709319 Å) versus Bragg angle 2Θ and temperature. The Bragg reflections shown are indexed as 200 at 18.15◦, 114 at 19.58◦, 105 at 20.60◦, 211 at 20.65◦, 203 at 21.30◦, and 212 at 21.63◦.

**X-ray measurements:** in order to check the structural aspect, we performed temperaturedependent X-ray powder diffraction measurements at 295, 100, 50, 20, and 5 K. The X-ray patterns were collected on a CaMnBi2 sample contained in a 0.3-mm diameter quartz glass capillary under He exchange gas using Mo *K*α<sup>1</sup> radiation. Temperatures between 295 and 5 K were adjusted in a home-built cryostat. As revealed in Figure A1e, there are no visible splittings or broadenings of the Bragg reflections, which would be indicative of a structural phase transition. The tetragonal lattice parameters were obtained from Rietveld profile refinements of the diffraction patterns assuming the space group P4/*nmm* (No. 129) and the atom and lattice parameters reported by Brechtel et al. [31] as starting parameters. They perfectly follow a simple Debye law:

$$\mathbf{V(T)} = \mathbf{V}\_0 + \mathbf{I}\_\mathbf{V} \mathbf{T} \frac{\mathbf{T}}{\Theta\_\mathbf{D}^3} \int\_0^{\Theta\_\mathbf{D}/\mathbf{T}} \frac{\mathbf{x}^3}{\mathbf{e}^\mathbf{x} - 1} d\mathbf{x},\tag{A1}$$

with V0 = 223.58(3) Å, Θ<sup>D</sup> = 279(4) K, and Iv = 0.0130(7) Å (here, V0 is the unit-cell volume at 0 K, Θ<sup>D</sup> is the Debye temperature, and the pre-factor Iv is a linear function of the Grüneisen parameter in the Debye approximation [32]). As seen from Figure A1e, no anomalies in the thermal expansion, which could indicate a structural phase transition, were detected. This is consistent with the absence of broadenings or splittings of the Bragg reflections. Thus, any detectable structural transition at Ts is excluded.

#### **References**


### *Article* **Topological Properties in a Λ/***V***-Type Dice Model**

**Shujie Cheng and Xianlong Gao \***

Department of Physics, Zhejiang Normal University, Jinhua 321004, China **\*** Correspondence: gaoxl@zjnu.edu.cn

**Abstract:** We studied a non-interacting Λ/*V*-type dice model composed of three triangular sublattices. By considering the isotropic nearest-neighbor hoppings and the next-nearest-neighbor hoppings with the phase, as well as the quasi-staggered on-site potential, we acquired the full phase diagrams under the different fillings of the energy bands. There are abundant topological non-trivial phases with different Chern numbers *C* = ±1, as well as higher ones ±2, ±3 and a metal phase in several regimes. In addition, we also checked the bulk–edge correspondence of the system by analyzing the edge-state energy spectrum.

**Keywords:** band structures; high Chern numbers; bulk-edge correspondence

#### **1. Introduction**

By means of an environment with a low temperature and strong magnetic field, Klitzing et al. discovered the quantum Hall effect in 1980 [1]. This discovery kicked off a wave of research on the quantum Hall effect, and a series of research on the topological features of condensed matter was inspired in the following decades [2–4]. These types of condensed matter are well classified according to their symmetry [5]. With these efforts, people have found some new topological quantum matter [6–8] with a quantum anomalous Hall effect (QAHE) [9], which releases the harsh condition of realizing a strong magnetic field.

A Chern insulator is a kind of insulator with QAHE that breaks the time-reversal symmetry. Its topological features can be directly reflected by the topological invariant, i.e., the Chern number (*C*). Thouless–Kohmoto–Nightingale–den Nijs (TKNN) first used the Chern number to describe the topological properties of a two-dimensional system [10]. In a gapped system, when the Fermi energy lies in a bulk band gap, the Chern number is always an integer and is equal to the quantized Hall conductance in units of *e*2/*h* [11–13]. Chern numbers can be used to distinguish whether or not a system has topological properties. In other words, *C* = 0 (*C* = 0) corresponds to the topological nontrivial (trivial) phase. Afterwards, by developing TKNN's theory, Berry provided an alternative way to calculate the topological invariant with the Berry gauge field in the Brillouin zone [14]. Interestingly, the Chern numbers are linked to gapless edge states, forming the so-called bulk-edge correspondence [15]. The magnitude of the Chern number indicates the number of edge states.

In 1988, Haldane first theoretically pioneered the idea of breaking the time-reversal symmetry by applying a zero net magnetic flux through each unit cell in a hexagonal lattice and engineered a topological nontrivial model with *C* = ±1, which is known as the Haldane model [16]. This model opens a gate for people to study QAHE and has more or less influenced other two-dimensional systems that appeared later, such as the Checkerboard lattice [17], Kagomé lattice [18–21] and Lieb lattice [22–25]. Moreover, by considering the long-range tunneling [26–28] or more a complex magnetic flux [29–31], one can obtain the topological phase with a higher Chern number.

Recently, people have realized the Haldane model by trapping ultracold atoms in an optical lattice formed by three standing-wave laser beams [32–34] or in a periodically modulated optical honeycomb lattice [35]. In addition to these experiments, there have

**Citation:** Cheng, S.; Gao, X. Topological Properties in a Λ/*V*-Type Dice Model. *Crystals* **2021**, *11*, 467. https://doi.org/10.3390/cryst11050467

Academic Editors: Artem Pronin and Kwang-Yong Choi

Received: 30 March 2021 Accepted: 19 April 2021 Published: 22 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

been several works about the dice model [30,31] with a nonzero phase in the process of nearest-neighbor (NN) hoppings in the dice lattice [36–43]. In this paper, we will focus on a non-interacting dice model with a nonzero phase in the next-nearest-neighbor (NNN) term, but not in the NN term. We will consider a Λ/*V*-type on-site potential, which will give us a more fruitful phase diagram and exciting topological phenomena.

By analyzing the dispersions of the bands, in this paper, we uncover that the system possesses two different phases—a metal phase and bulk insulating phase, which are both in the 1/3 filling and 2/3 filling cases. In the bulk insulating phase, the topological properties are further investigated. After calculating the Chern number, we find that there is a topological trivial phase with *C* = 0 and topological nontrivial phases with *C* = ±1, *C* = ±2, and *C* = ±3, which are separated by the band-crossing lines. Then, we obtain the full phase diagram and find that the metal phase and topological phases are symmetrically distributed within the phases. The relative positions of these phases in the phase diagram will be inverted when we change the filling from 2/3 to 1/3, or vice verse. Finally, we also solve the edge-state spectrum and check the rule of bulk-edge correspondence.

This paper is organized as follows. In Section 2, we present the Hamiltonian of the dice model in both real and momentum space. In Section 3, we first analyze the band structures of the system. Next, we calculate the Chern numbers and numerically obtain the full phase diagram. Further, we check the bulk-edge correspondence through the edge-state spectrum. We summarize our work in the final section.

#### **2. Model and Hamiltonian**

Here, a non-interacting Λ/*V*-type dice model is studied, and it is shown in Figure 1. There are three interpenetrating triangular sublattices, denoted by *R* (red dots), *B* (blue dots), and *G* (green dots). The lattice constant *a* is set as *a* = 1. The Hamiltonian of our model has three parts:

$$
\hat{H} = \hat{H}\_1 + \hat{H}\_2 + \hat{H}\_3.\tag{1}
$$

The first part *H*ˆ <sup>1</sup> shows the isotropic hoppings with the same hopping amplitudes between nearest-neighbor sites, which pertain to various sublattices, and this part is presented as 

$$\begin{split} \hat{H}\_{1} &= \sum\_{\langle \mathbf{R}\_{i}, \mathbf{B}\_{j} \rangle} t \Big( \boldsymbol{\varepsilon}\_{\mathbf{R}\_{i}}^{\dagger} \mathbf{\hat{c}}\_{\mathbf{B}\_{j}} + H.c. \Big) \\ &+ \sum\_{\langle \mathbf{R}\_{i}, \mathbf{G}\_{\ell} \rangle} t \Big( \boldsymbol{\varepsilon}\_{\mathbf{R}\_{i}}^{\dagger} \mathbf{\hat{c}}\_{\mathbf{G}\_{\ell}} + H.c. \Big) \\ &+ \sum\_{\langle \mathbf{B}\_{j}, \mathbf{G}\_{\ell} \rangle} t \Big( \boldsymbol{\varepsilon}\_{\mathbf{B}\_{j}}^{\dagger} \mathbf{\hat{c}}\_{\mathbf{G}\_{\ell}} + H.c. \Big), \end{split} \tag{2}$$

in which *t* denotes the hopping amplitude, and it is regarded as the unit of energy, *c*ˆ**R***<sup>i</sup>* , *c*ˆ**B***<sup>j</sup>* . *c*ˆ**G** are the fermionic annihilation operators, and **R***i*, **B***j*, and **G** denote the coordinates of relevant sublattice sites *R*, *B*, and *G*, respectively. · · · is the nearest-neighbor relation.

In the same *R* and *B* sublattices, we will consider the next-nearest-neighbor hoppings accompanied by a phase. Then, *H*ˆ <sup>2</sup> can be expressed as

$$\begin{split} \hat{H}\_{2} &= \sum\_{\langle \mathbf{R}\_{i}, \mathbf{R}\_{j} \rangle} \Big( t\_{2} \varepsilon^{i\phi} \mathcal{E}\_{\mathbf{R}\_{i}}^{\dagger} \mathcal{E}\_{\mathbf{R}\_{j}} + H.c. \Big) \\ &+ \sum\_{\langle \mathbf{B}\_{i}, \mathbf{B}\_{j} \rangle} \Big( t\_{2} \varepsilon^{i\phi} \mathcal{E}\_{\mathbf{B}\_{i}}^{\dagger} \mathcal{E}\_{\mathbf{B}\_{j}} + h.c. \Big), \end{split} \tag{3}$$

in which *<sup>t</sup>*2*e*±*i<sup>φ</sup>* is the hopping amplitude, where *<sup>φ</sup>* is the phase, and <sup>±</sup> represents the direction of the hoppings (+ is the clockwise direction and − is the counterclockwise direction). References [30–34,44,45] tell us that the hopping terms in our toy models can be realized through laser-assisted tunneling. Moreover, the phase accompanied by the tunneling can be modulated with the momentum recoil. We think that this method is helpful in selectively engineering such a phase in specific hopping terms in experiments. At the very least, these are the most natural imperfections that may arise in the experimental implementation of such a model.

**Figure 1.** (**a**) Geometry of the dice lattice. *R*, *B*, and *G* are sublattice sites, and are marked with red, blue, and green dots respectively. The lattice constant *a* is taken as *a* = 1. The vectors *an* (*n* = 1, 2, 3) link the nearest-neighbor sites that pertain to different sublattices. The vectors *bn* (*n* = 1, 2, 3) link the nearest-neighbor *R* sites or *B* sites. The circle arrows represent the next-nearest-neighbor tunnelings accompanied by a phase *eiφ*. (**b**) The first Brillouin zone of the model. Γ, K, and M are high-symmetry points in the high-symmetry path, and are connected by three red dashed lines.

The final part describes the on-site potentials with a special configuration

$$
\hat{H}\_3 = \Delta \sum\_{\mathbf{R}\_i} \mathcal{E}\_{\mathbf{R}\_i}^\dagger \mathcal{E}\_{\mathbf{R}\_i} + \Delta \sum\_{\mathbf{B}\_i} \mathcal{E}\_{\mathbf{B}\_i}^\dagger \mathcal{E}\_{\mathbf{B}\_i} - 2\Delta \sum\_{\mathbf{G}\_i} \mathcal{E}\_{\mathbf{G}\_i}^\dagger \mathcal{E}\_{\mathbf{G}\_i} \tag{4}
$$

in which Δ denotes the potential at the R and B sublattice sites, and −2Δ is the potential at the G sublattice sites. The configuration of the potential can be viewed as Λ/*V*-type (Δ < 0/Δ > 0) of three levels in a super atom, which is composed of three sites and forms a three-band model in the lattice case. This potential configuration differs from that in other dice models [30,31] and can also be realized by tuning single-beam lattice depths [34,35].

We consider a system that has discrete translational symmetry; thus, the single-particle Hamiltonian can be written in momentum space [46–48] as

$$
\hat{\mathcal{H}}(\mathbf{k}) = I(\mathbf{k}) + \mathbf{d}(\mathbf{k}) \cdot \vec{\lambda}, \tag{5}
$$

where *I*(**k**) is a scalar, **d**(**k**) is a real vector with eight components, and- *λ* denotes a vector consisting of Gell–Mann matrices [49]. As a matter of fact, the Chern number will not be affected by the scalar *I*(**k**), and is only determined by **d**(**k**). In order to obtain the <sup>H</sup>ˆ(**k**), we need to perform a discrete Fourier transformation on the three-component basis (*c*ˆ**k**,*R*, *c*ˆ**k**,*B*, *c*ˆ**k**,*G*) *T*:

$$
\begin{split}
\hat{\mathbf{c}}\_{\mathbf{k},R} &= \frac{1}{\sqrt{N}} \sum\_{\mathbf{R}\_{\hat{j}}} e^{-i\mathbf{k}\cdot\mathbf{R}\_{\hat{j}}} \hat{\mathbf{c}}\_{\mathbf{R}\_{\hat{j}'}} \\
\hat{\mathbf{c}}\_{\mathbf{k},B} &= \frac{1}{\sqrt{N}} \sum\_{\mathbf{B}\_{\hat{j}}} e^{-i\mathbf{k}\cdot\mathbf{B}\_{\hat{j}}} \hat{\mathbf{c}}\_{\mathbf{B}\_{\hat{j}'}} \\
\hat{\mathbf{c}}\_{\mathbf{k},G} &= \frac{1}{\sqrt{N}} \sum\_{\mathbf{G}\_{\hat{j}}} e^{-i\mathbf{k}\cdot\mathbf{G}\_{\hat{j}}} \hat{\mathbf{c}}\_{\mathbf{G}\_{\hat{j}}}.
\end{split}
\tag{6}
$$

After the derivation, we acquire the components of the vector **d**(**k**), which are presented as

$$\begin{aligned} d\_1 &= d\_6 = d\_4 = \sum\_n t \cos(\mathbf{k} \cdot a\_n), \\ d\_2 &= d\_7 = -d\_5 = \sum\_n t \sin(\mathbf{k} \cdot a\_n), \\ d\_3 &= -2t\_2 \sin \phi \sum\_n \sin(\mathbf{k} \cdot b\_n), \\ d\_8 &= \sqrt{3}\Delta + \frac{2t\_2}{\sqrt{3}} \cos \phi \sum\_n \cos(\mathbf{k} \cdot b\_n), \end{aligned} \tag{7}$$

in which the six-unit vectors *an* and *bn* (*n* = 1, 2, 3), which are shown in Figure 1a, are listed as

$$\begin{aligned} a\_1 &= \begin{pmatrix} 0 \\ -1 \end{pmatrix}, & a\_2 &= \frac{1}{2} \begin{pmatrix} \sqrt{3} \\ 1 \end{pmatrix}, & a\_3 &= \frac{1}{2} \begin{pmatrix} -\sqrt{3} \\ 1 \end{pmatrix}, \\\ b\_1 &= \begin{pmatrix} \sqrt{3} \\ 0 \end{pmatrix}, & b\_2 &= \frac{1}{2} \begin{pmatrix} -\sqrt{3} \\ 3 \end{pmatrix}, & b\_3 &= -\frac{1}{2} \begin{pmatrix} \sqrt{3} \\ 3 \end{pmatrix}. \end{aligned} \tag{8}$$

#### **3. Results and Discussion**

*3.1. Band Structures*

To begin, we study the band structures of the model described by Equation (1). The first Brillouin zone of this model is shown in Figure 1b, with Γ, K, and M being the highsymmetry points [50]. Without loss of generality, in the following numerical calculation, we set *t* = *t*<sup>2</sup> = 1. With the known components of **d**(**k**) in Equation (7), we can diagonalize <sup>H</sup>ˆ(**k**) and obtain its eigenvalues at each momentum **<sup>k</sup>**. According to this strategy, we finally acquire the energy dispersions along the high-symmetry path Γ-K-M-Γ with different Δ and phase *φ*.

As a matter of fact, the energy at each band has a maximal value or a minimal value at a high-symmetry point. According to this obvious feature, we numerically analyze the energy gap between two adjacent energy bands, and finally uncover that the system consists of a bulk insulating phase and a metal phase. Figure 2 shows the metal–insulator phase diagram at 1/3 filling, and Figure 3 shows that for 2/3 filling. The green dots, which are marked as *a*, *b*, *c*, and *d*, are the four chosen typical parameter points. The terms 1/3 filling and 2/3 filling mean that the Fermi energies are selected to ensure that these three bands are filled just enough by 1/3 and 2/3, respectively. There is no doubt that the system is in the metal phase (surrounded by the blue solid line and marked by *M*) when the gap is closed. The bulk insulating phase (marked by *I*) requires that the gap is be open. For instance, when the parameter is tuned to the *a* point, the gap is closed at 1/3 filling and is open at 2/3 filling. The consequence at the *b* point is the opposite of that at the *a* point. The system remains gapless at the *c* point and remains gapped at the *d* point. It is worth noting that the metal phase and bulk insulating phase are symmetrically distributed.

**Figure 2.** The metal-insulator phase diagrams in the cases of 1/3 filling (**a**) and 2/3 filling (**b**). *M* refers to the metal phase and *I* refers to the bulk insulating phase. In each case, there are four chosen typical points, *a*, *b*, *c*, and *d*, which are marked by green dots and are discussed in the main text.

Intuitively, the dispersion curves at these four chosen parameter points, (Δ*a*, *φa*), (Δ*b*, *φb*), (Δ*c*, *φc*), and (Δ*d*, *φd*), are plotted in Figure 3a–d, respectively. In these figures, we choose two Fermi energies to ensure the 1/3 filling (black solid line) and the 2/3 filling (black dashed line), respectively. In Figure 3a, because the top of the lowest band is higher than the bottom of the middle band, the Fermi energy crosses the middle band and forms a fully occupied lowest band and a partially occupied middle band. Accordingly, the system presents a metallic property in the 1/3 filling case. If the Fermi energy is located at the 2/3 filling line, the middle band is fully occupied and the highest band is empty. Then, in this case, the system will be in the bulk insulating phase. For the same reasons, one can easily comprehend that the system is a band insulator at 1/3 filling and a metal at 2/3 filling in the case shown in Figure 3b. Moreover, in Figure 3c,d, the system is stable in the metal phase and bulk insulating phase, respectively, no matter what the filling is. All of the results agree with our metal–insulator phase diagrams in Figure 2.

**Figure 3.** Dispersions of the Λ/*V*-type dice model along the high-symmetry path Γ-K-M-Γ. (**a**) Δ*<sup>a</sup>* = −2, *φ<sup>a</sup>* = −2.5, (**b**) Δ*<sup>b</sup>* = 2, *φ<sup>b</sup>* = 0.5, (**c**) Δ*<sup>c</sup>* = −0.5, *φ<sup>c</sup>* = 0.1, (**d**) Δ*<sup>d</sup>* = −2, *φ<sup>d</sup>* = 1.5. The red, blue, and green solid lines show the dispersions. The lower Fermi energy (black dashed line) and higher Fermi energy (black solid line) show the cases of 1/3 and 2/3 filling, respectively.

#### *3.2. Chern Numbers and the Edge-State Spectrum*

The nontrivial topological features of the Haldane model [16] and other dice models [30,31,43] motivate us to make it clear what topological phases exist in the bulk insulating region. According to the energy band theory [10,11,14], the n-th occupied band's Chern number is defined as a contour integral along the boundary of the first Brillouin zone,

$$\mathbb{C}\_n = \frac{1}{2\pi} \oint \mathbf{A}\_n(\mathbf{k}) \cdot d\mathbf{k},\tag{9}$$

in which *n* ∈ {1, 2, 3} denotes the band index and **A***<sup>n</sup>* = −*iψn*(**k**)|∇**k**|*ψn*(**k**), with <sup>|</sup>*ψn*(**k**) being the associated eigenvector of <sup>H</sup>ˆ(**k**). The ascending order of *<sup>n</sup>* indicates the band from the bottom to the top. The topological properties of this model can be reflected by two quantities, i.e., *C*<sup>1</sup> 3 and *C*<sup>2</sup> 3 , which satisfy *C*<sup>1</sup> 3 = *C*<sup>1</sup> and *C*<sup>2</sup> 3 = *C*<sup>1</sup> + *C*2. We calculate the Chern number and obtain the full phase diagrams of the system, which are shown in Figure 4a,b for 1/3 and 2/3 filling of the system, respectively. There are several phase boundary lines in the bulk insulating phase, which are captured by the closing of the energy band. In Fact, the energy-crossing-lines also appear in the metal phase, but do not change the intrinsic properties of the metal phase. In each diagram, the black solid lines separate the topological nontrivial phase from the topological trivial phase, and the red solid lines distinguish the topological nontrivial phase with various *C*<sup>1</sup> 3 and *C*<sup>2</sup> 3 .

**Figure 4.** Two full phase diagrams with (**a**) 1/3 filling and (**b**) 2/3 filling, respectively. The blue solid lines and axes surround the metallic region (*M*), and the bulk insulating region consists of several topological regions. The black solid lines separate the topological nontrivial phase from the topological trivial phase, and the red solid lines distinguish the topological nontrivial phase with various *C*<sup>1</sup> 3 and *C*<sup>2</sup> 3 .

As can be seen from the phase diagrams, abundant quantum phases exist in our model. In Figure 4a, intuitively, we know that, except for the metal phase, there are topological nontrivial phases with *C*<sup>1</sup> 3 = *C*<sup>1</sup> = ±1 and *C*<sup>1</sup> 3 = *C*<sup>1</sup> = ±2, as well as the topological trivial phase with *C*<sup>1</sup> 3 = *C*<sup>1</sup> = 0. A similar circumstance also appears in Figure 4. Furthermore, when the parameters are tuned continuously, the system will undergo abundant phases. We take Δ = −0.7 as an example. As *φ* increases, at 2/3 filling, the system will undergo a closed circle with six different phases:

$$\begin{array}{ccccc} \mathbb{C}\_{\frac{2}{3}} = -1 & \Rightarrow & \mathbb{C}\_{\frac{2}{3}} = +2 \\ & \uparrow & & \Downarrow \\ & \mathcal{M} & & \mathcal{M} \\ & \uparrow & & \Downarrow \\ \mathbb{C}\_{\frac{2}{3}} = +1 & \Leftarrow & \mathbb{C}\_{\frac{2}{3}} = -2 \end{array}$$

where *M* stands for the metal phase. The distribution of the phase diagram in Figure 4b (2/3 filling) can be loosely regarded as the inversion of that with the 1/3 filling. Similarly, the system can also undergo rich phases when we tune the Fermi energy.

In the following, we check the rule of the bulk–edge correspondence. We keep the *x* direction of the system periodic and make the system possess a zigzag edge along the *y* direction. Therefore, *kx* is a good quantum number. Finally, we solved the edge-state spectrum [4,15], which is plotted in Figure 5. We find that the system also obeys the rule of bulk–edge correspondence [15] for the reason that the magnitudes of *C*<sup>1</sup> 3 and *C*<sup>2</sup> 3 can be reflected from the intersections of the edge-state spectra. When we choose the parameter (Δ, *φ*)=(−2, 1.5), the edge-state spectrum intersects at *kx* = ±*π*/ <sup>√</sup>3 (see Figure 5a), which means that there is only a pair of edge modes, corresponding to *C*<sup>1</sup> 3 = 1. The edge-state

spectrum is gapped at 2/3 filling, which means that there are no edge modes, corresponding to *C*<sup>2</sup> 3 = 0. Meanwhile, at this parameter point, the Chern number of the lowest band is *C*<sup>1</sup> = −1 and that of the middle band is *C*<sup>2</sup> = *C*<sup>2</sup> 3 − *C*<sup>1</sup> = 0 − (−1) = 1. In addition, it is known from the symmetry of the phase diagram that the Chern numbers become *C*<sup>1</sup> = 1 and *C*<sup>2</sup> = −1 at (Δ, *φ*)=(−2, −1.5), and will also obtain the same edge-state spectrum as that in Figure 5a. When we select the parameter at (Δ, *φ*) = (1, −1.3), there are two different intersections at 1/3 filling (see Figure 5b), which means that there are two pairs of edge modes, corresponding to *C*<sup>1</sup> 3 = 2. The spectrum intersects at *kx* = ±*π*/ <sup>√</sup>3, which means that there is a pair of edge modes, corresponding to *C*<sup>2</sup> 3 = −1. By the definitions of *C*<sup>1</sup> 3 and *C*<sup>2</sup> 3 , we get *C*<sup>1</sup> = 2 and *C*<sup>2</sup> = −3 in this case. Similarly, when we take (Δ, *φ*) = (1, −1.3), we will know that *C*<sup>1</sup> = −2 and *C*<sup>2</sup> = 3 and we will obtain the same correspondence as that shown in Figure 5b.

**Figure 5.** Two edge-state spectra of a cylindrical geometry with a zigzag edge. (**a**) Δ = −2, *φ* = 1.5. There is a pair of edge modes in the case of 1/3 filling, corresponding to *C*<sup>1</sup> 3 = −1, while there are no edge modes in the case of 2/3 filling, which corresponds to *C*<sup>2</sup> 3 = 0; (**b**) Δ = 1, *φ* = −1.3. There are two pairs of edge modes in the case of 1/3 filling, which corresponds to *C*<sup>1</sup> 3 = 2, and there is only a pair of edge modes in the case of 2/3 filling, which corresponds to *C*<sup>2</sup> 3 = −1.

#### **4. Conclusions**

To conclude, we have studied the band structures and topological properties of a Λ/V-type dice model. First, we investigated the energy spectrum characteristics under 1/3 filling and 2/3 filling. In the Δ-*φ* parameter space, the system can be divided into two parts: the metal phase and the bulk insulating phase. Furthermore, we calculated the Chern numbers according to the energy band theory in the bulk insulating phase and obtained a fruitful phase diagram. Interestingly, there are many topological nontrivial phases that are separated by some energy-level-crossing lines with different Chern numbers, such as *C* = ±1, *C* = ±2, and *C* = ±3. Finally, by solving the edge-state spectrum, we found that the system obeys the rule of bulk-edge correspondence.

In spite of the existence of dice lattice structures in several electronic materials, such as SrTiO3/SrIrO3/SrTiO3 [51], Ba2CoRe2O12 [52], and Gd2CCl2 [53–55], the high free modulation of the parameters in cold-atom experiments will be convenient for us to study the topological phases by manipulating the neutral atoms, which never occurred in the aforementioned research. For this reason, we hope that our system can be realized in cold-atom experiments and that phases with higher Chern numbers will be observed.

**Author Contributions:** Conceptualization, X.G. and S.C.; methodology, X.G. and S.C.; software, X.G. and S.C.; validation, S.C.; formal analysis, S.C.; investigation, S.C.; resources, X.G.; data curation, S.C.; writing—original draft preparation, S.C.; writing—review and editing, X.G. and S.C.; visualization, X.G. and S.C.; supervision, X.G.; project administration, X.G.; funding acquisition, X.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by NSFC under Grants No. 11835011 and No. 11774316.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the authors. The data are not publicly available due to intellectual property protection.

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

#### **References**


## *Article* **Low-Energy Optical Conductivity of TaP: Comparison of Theory and Experiment**

**Alexander Yaresko 1,\* and Artem V. Pronin <sup>2</sup>**


**\*** Correspondence: a.yaresko@fkf.mpg.de

**Abstract:** The ab-plane optical conductivity of the Weyl semimetal TaP is calculated from the band structure and compared to the experimental data. The overall agreement between theory and experiment is found to be best when the Fermi level is slightly (20 to 60 meV) shifted upwards in the calculations. This confirms a small unintentional doping of TaP, reported earlier, and allows a natural explanation of the strong low-energy (50 meV) peak seen in the experimental ab-plane optical conductivity: this peak originates from transitions between the almost parallel non-degenerate electronic bands split by spin-orbit coupling. The temperature evolution of the peak can be reasonably well reproduce by calculations using an analog of the Mott formula.

**Keywords:** Weyl semimetals; band-structure calculations; optical response

#### **1. Introduction**

Weyl fermions [1] are known to be observed as elementary excitations in condensedmatter systems—the Weyl semimetals (WSMs) [2–8]. In WSMs, the valence and conduction bands touch each other at selected points of the Brillouin zone (BZ), the Weyl nodes.

TaP belongs to the currently most studied family of WSMs, which also includes NbP, TaAs, and NbAs. These materials are nonmagnetic non-centrosymmetric WSMs with 24 Weyl nodes of two different types, usually dubbed as W1 (8 nodes) and W2 (16 nodes) [5,9–13]. The available band-structure calculations predict that in TaP the W1 nodes are situated some 40 to 55 meV below the Fermi level *EF*, while the W2 nodes are at 12 to 20 meV above it [11–13], see Figure 1.

The low-energy band structure of TaP and other WSMs determines their peculiar physical properties [14]. One way to experimentally probe the band structure at low energies is optical spectroscopy in the infrared (particularly, in the far-infrared) region [15]. The frequency-dependent conductivity, *σ*(*ω*) = *σ*1(*ω*) + i*σ*2(*ω*), of three-dimensional linear bands has been well studied theoretically using model Hamiltonians [16–23]. It has been shown that the interband portion of *σ*1(*ω*) for a single isotropic Weyl band has to follow a linear frequency dependence [16–18]:

$$
\sigma\_1(\omega) = \frac{\varepsilon^2}{12h} \frac{\omega}{v\_F} \tag{1}
$$

where *vF* is the band Fermi velocity (this formula is obtained assuming the electron-hole symmetry). Such linear behavior of *σ*1(*ω*) at low energies has indeed been observed in a number of established three-dimensional Weyl and Dirac semimetals, as well as in candidate materials [24–30].

In many real materials of this type, however, the linear interband conductivity at low energies is (partly) masked by other features, such as intraband (Drude) conductivity or resonance-like interband contributions [30–35]. Particularly strong low-energy peak was

**Citation:** Yaresko, A.; Pronin, A.V. Low-Energy Optical Conductivity of TaP: Comparison of Theory and Experiment. *Crystals* **2021**, *11*, 567. https://doi.org/10.3390/cryst 11050567

Academic Editor: Kwang-Yong Choi

Received: 25 March 2021 Accepted: 17 May 2021 Published: 19 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

observed in TaP [30,35]. In Ref. [30], this peak was assigned to electron-hole pair excitations near the saddle points of the crossing bands, which form the Weyl nodes (Figure 1). To our knowledge, this assignment doesn't have a direct support from optical-conductivity calculations based on band structure. Also, the total number of states near the saddle points is relatively low. Hence, only relatively small kinks in the optical conductivity, rather than strong peaks, are expected in this situation [20,34]. On the other hand, our earlier study of the sister compound NbP [34] has demonstrated that the low-energy peaks, similar to the one in TaP, appear in NbP and are due to multiple transitions between almost parallel bands split by spin-orbit coupling (SOC). Based on our band structure calculations, we argue in this paper that the same explanation of its low-energy peak holds also for TaP.

**Figure 1.** (**a**) Crystallographic structure of TaP and (**b**) a schematic diagram of its Weyl bands. Possible transitions between the saddle points of the merging Weyl bands and between the SOC-split bands are indicated as arrows.

#### **2. Results and Discussion**

In Figure 2, we plot our experimental optical spectra presented earlier in Ref. [35]. The measurements have been done on the isotropic *ab* plane of TaP (cf. Figure 1). The prominent low-energy peak is clearly seen in the real part of the optical conductivity at 50 meV (it corresponds to a deep in the optical reflectivity).

To gain insight into the origin of this peak, we carried out band structure calculations within the local density approximation (LDA) based on the experimental crystal structure of TaP [36]. The calculations were performed using the linear muffin-tin orbital (LMTO) method [37] with the Perdew-Wang exchange-correlation potential [38]. We used the relativistic PY LMTO computer code [39] with SOC added to the LMTO Hamiltonian in the variational step. BZ integrations were done using the improved tetrahedron method [40]. Additional empty spheres (E) were inserted at the 8b Wyckoff positions in order to minimize the effect of atomic sphere overlap. The Ta, P, and E states up to the maximal orbital quantum number *lmax* = 3, 2, and, 1, respectively, were included into the LMTO basis set which is essential for calculation of the dipole matrix elements for the interband transitions involving the Ta *d*- and the P *p*-derived bands. When calculating the real part of the optical conductivity, we used the tetrahedron method on a dense 128 × 128 × 128 *k*-mesh in order to resolve interband transitions between the SOC-split bands close to Weyl points [31,34]. No broadening has been applied to the computed spectra.

In Figure 3 we show the results of our band-structure calculations as well as the BZ of TaP. Four non-spin-degenerate bands, numbered 19 to 22, can be resolved in the vicinity of *EF* (at every given **k** point, the bands are numbered with increasing energy staring from the lowest calculated band). Note that the bands in each of the two pairs, (19, 20) and (21, 22), are split by SOC because of the lack of inversion symmetry. Our results reproduce well the published band structures of TaP calculated using the full-potential codes [11,13].

**Figure 2.** (**a**) Experimental in-plane reflectivity and (**b**) the corresponding real part of the optical conductivity of TaP at selected temperatures [35]. The arrows mark the feature discussed in this paper. The increased *σ*1*xx* at low energies is due to intraband (Drude-like) absorption.

**Figure 3.** (**a**) Brillouin zone of TaP. (**b**) Fermi surface cross sections calculated for Δ*EF* = 50 meV. Black circles show approximate positions of projections of Weyl points onto *ky* = 0 plane. (**c**) Band structure of TaP. Four bands closest to *EF* (marked 19 to 22) are shown in different colors. Black and red horizontal dashed lines show the as-calculated position of *EF* and the Fermi level shifted upwards by Δ*EF* = 50 meV, respectively.

Before we discuss the calculated optical conductivity spectra, we would like to note that in WSMs the match between the measured and the calculated optical conductivity at low frequencies is typically only qualitative: the calculations catch the major features observed in the experimental spectra, but are unable to reproduce the exact frequency positions of the features and their spectral shapes [13,15,30,31,34].

Another important point to be mentioned here is the unintentional (self-)doping, which is inherent to many real materials, where impurities, crystallographic defects, and vacancies may slightly change the position of *EF*. Such unintentional doping, varying from sample to sample, has been shown to be relevant to TaP [41]. On the other hand, band structure calculations themselves have finite accuracy (cf. the spread in the calculated energy positions of the Weyl nodes, mentioned above). These considerations justify small variation of the position of *EF* to get a better match between theory and experiment. Hereafter, we vary the Fermi-level position within Δ*EF* = ±100 meV.

Figure 4 presents the results of our interband-conductivity calculations starting from the self-consistent band structure shown in Figure 3. The black solid line in Figure 4 shows *σ*1*xx*(*h*¯ *ω*) obtained with as-calculated *EF*. The overall run of the experimental interband conductivity is well reproduced by this curve: *σ*1*xx*(*h*¯ *ω*) increases with frequency and reaches a maximum at 1 eV (cf. Figure 2 and note that the intraband (Drude) contribution has not been taken into account in the band-structure computations). Nevertheless, no peak is visible in these computations at around 50 meV. A slight variation of *EF* provides such a peak, but only if the Fermi level is shifted upwards (red and blue curves). Shifting *EF* downwards does not change the *σ*1*xx* spectra in the desirable way (magenta and cyan curves). The hight of the 50-meV peak reaches the experimental value of 2.5 × 103 <sup>Ω</sup>−1cm−<sup>1</sup> at Δ*EF* = 50 meV.

**Figure 4.** Low-energy optical conductivity of TaP calculated from its band structure. Lines of different colors correspond to different positions of the Fermi level, as indicated. The conductivity calculated for smaller positive Δ*EF* is plotted in the inset. The contributions of 21 → 22 transitions are shown by dashed lines.

In the inset in Figure 4, we present an expanded view of the low-frequency optical conductivity calculated for small positive Δ*EF*. It is obvious, that already a very small *EF* shift of 20 meV is sufficient to produce the 50-meV peak. Note, that for all three curves yet another experimental feature—a broad shoulder at 0.3–0.5 eV—is also evident in the calculated spectra. Thus, we can conclude that a tiny shift of the Fermi level allows one to obtain a very reasonable overall description of the experimental *σ*1*xx*(*h*¯ *ω*), including the strong peak at 50 meV.

To understand what interband optical transitions are responsible for this peak, one can take a look at Figure 3, where the original and shifted by 50 meV Fermi level positions are shown by black and red dashed lines, respectively. In the vicinity of Weyl points, i.e., near the S point and along the N–M line, band 21 is above the as-calculated *EF*. The lowfrequency interband conductivity is dominated by the transitions between the initial band 20 and the final band 21. The shift of *EF* to higher energy leads to partial occupation of band 21. This suppresses the 20 → 21 transitions at low energy and, at the same time, allows transitions from band 21 to band 22, which remains mostly empty. As these SOC-split bands are almost parallel, the energies of such transitions are expected to be roughly the same for different momenta. Thus, a sharp peak may occur in *σ*1*xx*(*h*¯ *ω*).

To confirm this observation, we performed band-resolved optical-conductivity calculations for the transitions between bands 21 and 22. The results of these calculations for three Δ*EF* are plotted by dashed lines in the inset of Figure 4. It is apparent, that the 21 → 22 transitions provide the major contribution to the 50-meV peak, confirming our proposition. A 21 → 22 contribution coming from the **k** volume near the middle of the Γ–X line appears also for the as-calculated *EF*, but it is too weak to be responsible for the 50-meV peak.

In order to model the experimentally observed temperature evolution of the 50 meV peak, we introduced a temperature dependence of the calculated interband optical conductivity by multiplying the interband transition probabilities with the factor *f*(*εi***k**)[1 − *f*(*ε<sup>f</sup>* **<sup>k</sup>**)], where *f*(*ε*) is the Fermi-Dirac distribution function and *εi***<sup>k</sup>** and *ε<sup>f</sup>* **<sup>k</sup>** are the energies of initial and final states, respectively. Figure 5a shows that even this simple approach allows one to reproduce the experimentally observed reduction of the 50-meV peak with increasing temperature. A better agreement between theory and experiment is obtained, if the optical conductivity is calculated using an analog of the Mott formula [42], which is widely used to study the thermoelectric properties of metals. In this approximation,

$$
\sigma(\omega) = \int \sigma(E, \omega) \left( -\frac{\partial f(E)}{\partial E} \right) dE\_\prime \tag{2}
$$

where *<sup>∂</sup> <sup>f</sup>*(*E*) *<sup>∂</sup><sup>E</sup>* is the energy derivative of the Fermi-Dirac function and *σ*(*E*, *ω*) is calculated with *E* being the energy which discriminates between the initial and final bands, so that *E* = *EF* at *T* = 0. The results of these computations are shown in Figure 5b. We note that we compute the temperature dependence of the interband contribution only. In order to reproduce the upturn of the measured conductivity at low photon energies, one needs to consider the temperature dependence of the intraband Drude term, which is beyond the scope of this work.

**Figure 5.** Temperature dependence of the optical conductivity calculated (**a**) by multiplying the interband transition probabilities with the Fermi-Dirac function and (**b**) using the Mott formula.

Finally, we would like to emphasize the importance of the transitions between the SOC-split bands. Such transitions can be considered forbidden in some models [20], while in the real WSMs they play an important role, as we have shown earlier for NbP [34].

These transitions are allowed, because the electronic bands can be characterized by their well-defined spin polarization, *s*±1/2, only for **k**-vectors faraway from the Weyl nodes; closer to the nodes, SOC is strong and spin polarization is much less perfect. Thus, transitions between *any* pair of bands are allowed there.

#### **3. Conclusions**

Summarizing, we have calculated the low-energy optical conductivity of the Weyl semimetal TaP (in the ab plane) and compared it to the experimental results. The best match between theory and experiment is found for a slightly shifted Fermi level (+20 to 60 meV). This shift confirms a small unintentional doping of TaP, discussed earlier [35,41], and offers a natural explanation of the strong low-energy (50 meV) peak reported in the experimental data [30,35]: the peak is due to transitions between the almost parallel non-degenerate electronic bands split by spin-orbit coupling.

**Author Contributions:** Conceptualization, A.Y. and A.V.P.; calculations, A.Y.; providing experimental data, A.V.P.; writing, A.Y. and A.V.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded in part by the Deutsche Forschungsgemeinschaft via grant number DR228/51-3.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are available upon request to the corresponding author.

**Acknowledgments:** We thank Sascha Polatkan for technical assistance in manuscript preparation.

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

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Crystals* Editorial Office E-mail: crystals@mdpi.com www.mdpi.com/journal/crystals

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-4236-2