# Observation of the Galactic Center PeVatron Beyond 100 TeV with HAWC

A. Albert  Physics Division, Los Alamos National Laboratory, Los Alamos, NM,
USA  R. Alfaro  Instituto de Física, Universidad Nacional Autónoma de México,
Ciudad de México, México  C. Alvarez  Universidad Autónoma de Chiapas, Tuxtla
Gutiérrez, Chiapas, México  A. Andrés  Instituto de Astronomía, Universidad
Nacional Autónoma de México, Ciudad de México, México  J.C. Arteaga-Velázquez
Universidad Michoacana de San Nicolás de Hidalgo, Morelia, México  D. Avila
Rojas  Instituto de Física, Universidad Nacional Autónoma de México, Ciudad de
México, México  H.A. Ayala Solares  Department of Physics, Pennsylvania State
University, University Park, PA, USA  R. Babu Department of Physics and
Astronomy, Michigan State University, East Lansing, MI, USA  E. Belmont-Moreno
Instituto de Física, Universidad Nacional Autónoma de México, Ciudad de
México, México  A. Bernal  Instituto de Astronomía, Universidad Nacional
Autónoma de México, Ciudad de México, México  K.S. Caballero-Mora  Universidad
Autónoma de Chiapas, Tuxtla Gutiérrez, Chiapas, México  T. Capistrán
Instituto de Astronomía, Universidad Nacional Autónoma de México, Ciudad de
México, México  A. Carramiñana  Instituto Nacional de Astrofísica, óptica y
Electrónica, Puebla, México  S. Casanova  Instytut Fizyki Jadrowej im Henryka
Niewodniczanskiego Polskiej Akademii Nauk, IFJ-PAN, Krakow, Poland  U. Cotti
Universidad Michoacana de San Nicolás de Hidalgo, Morelia, México  J. Cotzomi
Facultad de Ciencias Físico Matemáticas, Benemérita Universidad Autónoma de
Puebla, Puebla, México  S. Coutiño de León  Department of Physics, University
of Wisconsin-Madison, Madison, WI, USA  E. De la Fuente  Departamento de
Física, Centro Universitario de Ciencias Exactase Ingenierias, Universidad de
Guadalajara, Guadalajara, México  C. de León  Universidad Michoacana de San
Nicolás de Hidalgo, Morelia, México  D. Depaoli  Max-Planck Institute for
Nuclear Physics, 69117 Heidelberg, Germany  N. Di Lalla  Department of
Physics, Stanford University: Stanford, CA 94305–4060, USA  R. Diaz Hernandez
Instituto Nacional de Astrofísica, óptica y Electrónica, Puebla, México  B.L.
Dingus  Physics Division, Los Alamos National Laboratory, Los Alamos, NM, USA
M.A. DuVernois  Department of Physics, University of Wisconsin-Madison,
Madison, WI, USA  J.C. Díaz-Vélez  Department of Physics, University of
Wisconsin-Madison, Madison, WI, USA  K. Engel  Department of Physics,
University of Maryland, College Park, MD, USA  T. Ergin  Department of Physics
and Astronomy, Michigan State University, East Lansing, MI, USA  C. Espinoza
Instituto de Física, Universidad Nacional Autónoma de México, Ciudad de
México, México  K.L. Fan Department of Physics, University of Maryland,
College Park, MD, USA  K. Fang  Department of Physics, University of
Wisconsin-Madison, Madison, WI, USA  N. Fraija  Instituto de Astronomía,
Universidad Nacional Autónoma de México, Ciudad de México, México  S. Fraija
Instituto de Astronomía, Universidad Nacional Autónoma de México, Ciudad de
México, México  J.A. García-González  Tecnologico de Monterrey, Escuela de
Ingenieria y Ciencias, Monterrey, N. L., México, 64849  F. Garfias  Instituto
de Astronomía, Universidad Nacional Autónoma de México, Ciudad de México,
México  H. Goksu  Max-Planck Institute for Nuclear Physics, 69117 Heidelberg,
Germany  M.M. González  Instituto de Astronomía, Universidad Nacional Autónoma
de México, Ciudad de México, México  J.A. Goodman  Department of Physics,
University of Maryland, College Park, MD, USA  S. Groetsch  Department of
Physics, Michigan Technological University, Houghton, MI, USA  J.P. Harding
Physics Division, Los Alamos National Laboratory, Los Alamos, NM, USA  S.
Hernández-Cadena  Tsung-Dao Lee Institute & School of Physics and Astronomy,
Shanghai Jiao Tong University  I. Herzog  Department of Physics and Astronomy,
Michigan State University, East Lansing, MI, USA  J. Hinton  Max-Planck
Institute for Nuclear Physics, 69117 Heidelberg, Germany  D. Huang Department
of Physics, University of Maryland, College Park, MD, USA  F. Hueyotl-
Zahuantitla  Universidad Autónoma de Chiapas, Tuxtla Gutiérrez, Chiapas,
México  T.B. Humensky  NASA Goddard Space Flight Center, Greenbelt, MD 20771,
USA  P. Hüntemeyer  Department of Physics, Michigan Technological University,
Houghton, MI, USA  A. Iriarte  Instituto de Astronomía, Universidad Nacional
Autónoma de México, Ciudad de México, México  S. Kaufmann  Universidad
Politecnica de Pachuca, Pachuca, Hgo, México  D. Kieda  Department of Physics
and Astronomy, University of Utah, Salt Lake City, UT, USA  A. Lara  Instituto
de Geofísica, Universidad Nacional Autónoma de México, Ciudad de México,
México  W.H. Lee  Instituto de Astronomía, Universidad Nacional Autónoma de
México, Ciudad de México, México  J. Lee  University of Seoul, Seoul, Republic
of Korea  H. León Vargas  Instituto de Física, Universidad Nacional Autónoma
de México, Ciudad de México, México  J.T. Linnemann  Department of Physics and
Astronomy, Michigan State University, East Lansing, MI, USA  A.L. Longinotti
Instituto de Astronomía, Universidad Nacional Autónoma de México, Ciudad de
México, México  G. Luis-Raya  Universidad Politecnica de Pachuca, Pachuca,
Hgo, México  K. Malone  Physics Division, Los Alamos National Laboratory, Los
Alamos, NM, USA  O. Martinez  Facultad de Ciencias Físico Matemáticas,
Benemérita Universidad Autónoma de Puebla, Puebla, México  J. Martínez-Castro
Centro de Investigación en Computación, Instituto Politécnico Nacional, México
City, México.  J.A. Matthews  Dept of Physics and Astronomy, University of New
México, Albuquerque, NM, USA  P. Miranda-Romagnoli  Universidad Autónoma del
Estado de Hidalgo, Pachuca, México  J.A. Montes  Instituto de Astronomía,
Universidad Nacional Autónoma de México, Ciudad de México, México  J.A.
Morales-Soto  Universidad Michoacana de San Nicolás de Hidalgo, Morelia,
México  E. Moreno  Facultad de Ciencias Físico Matemáticas, Benemérita
Universidad Autónoma de Puebla, Puebla, México  M. Mostafá  Department of
Physics, Temple University, Philadelphia, Pennsylvania, USA  M. Najafi
Department of Physics, Michigan Technological University, Houghton, MI, USA
L. Nellen  Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de
México, Ciudad de México, México  M. Newbold  Department of Physics and
Astronomy, University of Utah, Salt Lake City, UT, USA  M.U. Nisa  Department
of Physics and Astronomy, Michigan State University, East Lansing, MI, USA  R.
Noriega-Papaqui  Universidad Autónoma del Estado de Hidalgo, Pachuca, México
L. Olivera-Nieto  Max-Planck Institute for Nuclear Physics, 69117 Heidelberg,
Germany  N. Omodei  Department of Physics, Stanford University: Stanford, CA
94305–4060, USA  M. Osorio-Archila  Instituto de Física, Universidad Nacional
Autónoma de México, Ciudad de México, México  Y. Pérez Araujo  Instituto de
Física, Universidad Nacional Autónoma de México, Ciudad de México, México
E.G. Pérez-Pérez  Universidad Politecnica de Pachuca, Pachuca, Hgo, México
C.D. Rho  Department of Physics, Sungkyunkwan University, Suwon 16419, South
Korea  D. Rosa-González  Instituto Nacional de Astrofísica, óptica y
Electrónica, Puebla, México  E. Ruiz-Velasco  Max-Planck Institute for Nuclear
Physics, 69117 Heidelberg, Germany  H. Salazar  Facultad de Ciencias Físico
Matemáticas, Benemérita Universidad Autónoma de Puebla, Puebla, México  D.
Salazar-Gallegos  Department of Physics and Astronomy, Michigan State
University, East Lansing, MI, USA  A. Sandoval  Instituto de Física,
Universidad Nacional Autónoma de México, Ciudad de México, México  M.
Schneider  Department of Physics, University of Maryland, College Park, MD,
USA  G. Schwefer  Max-Planck Institute for Nuclear Physics, 69117 Heidelberg,
Germany  J. Serna-Franco  Instituto de Física, Universidad Nacional Autónoma
de México, Ciudad de México, México  A.J. Smith  Department of Physics,
University of Maryland, College Park, MD, USA  Y. Son  University of Seoul,
Seoul, Republic of Korea  R.W. Springer  Department of Physics and Astronomy,
University of Utah, Salt Lake City, UT, USA  O. Tibolla  Universidad
Politecnica de Pachuca, Pachuca, Hgo, México  K. Tollefson  Department of
Physics and Astronomy, Michigan State University, East Lansing, MI, USA  I.
Torres  Instituto Nacional de Astrofísica, óptica y Electrónica, Puebla,
México  R. Torres-Escobedo  Tsung-Dao Lee Institute & School of Physics and
Astronomy, Shanghai Jiao Tong University  R. Turner  Department of Physics,
Michigan Technological University, Houghton, MI, USA  F. Ureña-Mena  Instituto
Nacional de Astrofísica, óptica y Electrónica, Puebla, México  E. Varela
Facultad de Ciencias Físico Matemáticas, Benemérita Universidad Autónoma de
Puebla, Puebla, México  X. Wang  Department of Physics, Michigan Technological
University, Houghton, MI, USA  Z. Wang  Department of Physics, University of
Maryland, College Park, MD, USA  I.J. Watson  University of Seoul, Seoul,
Republic of Korea  E. Willox  Department of Physics, University of Maryland,
College Park, MD, USA  H. Wu  Department of Physics, University of Wisconsin-
Madison, Madison, WI, USA  S. Yu  Department of Physics, Pennsylvania State
University, University Park, PA, USA  S. Yun-Cárcamo Department of Physics,
University of Maryland, College Park, MD, USA  H. Zhou  Tsung-Dao Lee
Institute & School of Physics and Astronomy, Shanghai Jiao Tong University

###### Abstract

We report an observation of ultra-high energy (UHE) gamma rays from the
Galactic Center region, using seven years of data collected by the High-
Altitude Water Cherenkov (HAWC) Observatory. The HAWC data are best described
as a point-like source (HAWC J1746-2856) with a power-law spectrum
($\mathrm{d}N/\mathrm{d}E=\phi(E/26\,\text{TeV})^{\gamma}$), where
$\gamma=-2.88\pm 0.15_{\text{stat}}-0.1_{\text{sys}}$ and $\phi=1.5\times
10^{-15}$ (TeV cm2s)-1
$\pm\,0.3_{\text{stat}}\,^{+0.08_{\text{sys}}}{}_{-0.13_{\text{sys}}}$
extending from 6 to 114 TeV. We find no evidence of a spectral cutoff up to
$100$ TeV using HAWC data. Two known point-like gamma-ray sources are
spatially coincident with the HAWC gamma-ray excess: Sgr A∗ (HESS J1745-290)
and the Arc (HESS J1746-285). We subtract the known flux contribution of these
point sources from the measured flux of HAWC J1746-2856 to exclude their
contamination and show that the excess observed by HAWC remains significant
($>$5$\sigma$) with the spectrum extending to $>$100 TeV. Our result supports
that these detected UHE gamma rays can originate via hadronic interaction of
PeV cosmic-ray protons with the dense ambient gas and confirms the presence of
a proton PeVatron at the Galactic Center.

##  1 Introduction

The Galactic sources of cosmic-ray acceleration to petaelectronvolt (PeV)
energies—known as PeVatrons—remain unidentified and are still subject to
discussion (Aharonian et al., 2019; Cristofari, 2021; Sudoh & Beacom, 2023; de
Oña Wilhelmi et al., 2024; Fang & Halzen, 2024; Blasi, 2013; Amato, 2014;
Gabici et al., 2019). Previous studies suggest that cosmic rays are actively
accelerated in the Galactic Center (GC) region (H.E.S.S. Collaboration et al.,
2016). The arc-minute angular resolution of the Imaging Atmospheric Cherenkov
telescopes (IACTs), e.g., the High Energy Stereoscopic System (H.E.S.S.;
Aharonian et al., 2006a; H.E.S.S. Collaboration et al., 2016; Abdalla et al.,
2018), the Major Atmospheric Gamma-Ray Imaging Cherenkov (MAGIC) telescopes
(Acciari et al., 2020), and the Very-Energetic Radiation Imaging Telescope
Array System (VERITAS; Adams et al., 2021), has allowed them to measure gamma-
ray emission up to $\sim$20 TeV from the two point sources of interest in the
region: Sgr A∗ (HESS J1745-290), the supermassive black hole at the center of
the Galaxy, and the unidentified source HESS J1746-285, which is spatially
coincident with the Galactic radio arc (Yusef-Zadeh & Morris, 1987; Yusef-
Zadeh et al., 2004). The observation of the point-like supernova remnant (SNR)
G0.9+0.1 (Abdalla et al., 2018; Acciari et al., 2020; Adams et al., 2021) and
the unidentified extended source HESS J1745-303 (Aharonian et al., 2006b) were
reported as well about $1^{\circ}$ away from the GC. These IACTs have also
observed very-high-energy (VHE) gamma rays from the GC ridge (H.E.S.S.
Collaboration et al., 2016; Abdalla et al., 2018; Adams et al., 2021; Acciari
et al., 2020). This diffuse emission spatially correlates to the Central
Molecular Zone (CMZ) morphology (Aharonian et al., 2006b), which is derived
from dense gas tracers (Tsuboi et al., 1999). This correlation suggests a
hadronic origin for the observed gamma-ray emission given the severe energy
losses via synchrotron emission in the leptonic scenario (Aharonian et al.,
2006a; H.E.S.S. Collaboration et al., 2016; Abdalla et al., 2018).

In this work, we use seven years of data from the HAWC Gamma-Ray Observatory
to study the gamma-ray emission from the GC region. Our analysis extends the
previous observations to energies $>$100 TeV, which allows the PeV cosmic-ray
interaction to be directly probed. We show that the UHE emission observed by
HAWC is most likely from the Galactic ridge emission by subtracting the flux
contribution from HESS J1745-290, as reported in H.E.S.S. Collaboration et al.
(2016), with good agreement in location and spectrum to other observations
(Adams et al., 2021; Acciari et al., 2020; Abe et al., 2023), and from HESS
J1746-285, as reported in Abdalla et al. (2018). The latter has been observed
by several IACTs also with good agreement, but only in Abdalla et al. (2018)
is the contribution of underlying diffuse emission additionally taken into
account for the source HESS J1746-285. The flux contribution of these two
sources at energies $>$100 TeV is extremely small. This indicates that these
sources do not contribute solely to the origin of the observed VHE gamma rays.
Our result provides evidence of a PeVatron at the center of our Galaxy with
the first measurement of nearly 100 gamma-ray events with energies $>$100 TeV.

This Letter is organized as follows: Section 2 briefly describes the HAWC data
set used in this analysis, Sections 3 and 4 present and discuss the results of
the analysis, and in Section 5 we present our conclusions.

##  2 HAWC data

The HAWC Observatory—located on the side of the Sierra Negra volcano in
Puebla, Mexico at 4100 m asl—is made up of 300 water Cherenkov detectors
(Abeysekara et al., 2023). We apply signal-topology-based cuts to reduce the
cosmic-ray background (99.9% of events detected). We recently updated HAWC’s
reconstruction algorithms (“Pass 5”), improving its effective area, angular
resolution, and gamma/hadron separation at the highest energies and zenith
angles. With these improvements, HAWC is able to observe the GC, which
culminates at 48∘ zenith (Albert et al., 2024). As a further check, we
verified that the results obtained when reconstructing data from the Crab
Nebula when it reaches zenith angles greater than $45^{\circ}$ are in
agreement with those reported in the study by Albert et al. (2024).

Using 2546 days of HAWC data, we detect gamma-ray emission from the GC region,
with a maximum significance of 6.5$\sigma$ above the background. We analyzed
the data with the Pass 5 version of the neural network energy estimator
(Abeysekara et al., 2019; Albert et al., 2024) and included off-array events,
which are showers whose cores fall off the main array up to 1.5 times its
physical area and improve the sensitivity of HAWC to high zenith angles and
high energies (Albert et al., 2024).

To model the gamma-ray flux from the GC region we employed the HAWC
Accelerated Likelihood (HAL) plugin with the Multi-Mission Maximum Likelihood
(3ML)111https://github.com/threeML/threeML framework (Younk et al., 2015;
Abeysekara et al., 2022)—a forward-folded maximum-likelihood approach
(Vianello et al., 2016)—within a rectangular region of interest (ROI) $\pm
3^{\circ}$ in latitude and $\pm 2.5^{\circ}$ in longitude.

We define our test statistics (TS) as:

$\text{TS}=2\ln{\left(\frac{L_{\text{model}}}{L_{\text{bkg}}}\right)}\,,$ (1)

where $L_{\text{model}}$ denotes the maximum likelihood from the source model
and $L_{\text{bkg}}$ is background only. According to Wilks’ theorem (Wilks,
1938), which applies to HAWC data (Abeysekara et al., 2017b), the TS is
asymptotically $\chi^{2}$ distributed, with the degrees of freedom equal to
the difference in the number of free parameters of the nested models. Thus,
under the case of one free parameter, $\sqrt{\text{TS}}$ can be used as a
measure of significance, $\sigma$.

The extended-source assumption was tested and no strong preference was found
($\Delta\text{TS}=6.62$ compared to the point-source assumption), hence the
simplest model was chosen. From the extended-source fit, we estimated an upper
limit (UL) on the source extension (radius) at the 68% confidence level (CL)
for the Gaussian width of the source ($\sim$0.46∘, see dashed circle in Figure
1(a)). Adding curvature to the spectrum did not significantly improve the test
statistic either ($\Delta\text{TS}=0.44$).

##  3 Main analysis results

The best fit to the data is a point source with a simple-power-law spectrum
(TS$=49$ for four free parameters—position and spectral parameters):

$\frac{\text{d}N}{\text{d}E}=\phi\left(\frac{E}{26\text{TeV}}\right)^{\gamma}\,,$
(2)

where $\phi$ is the flux normalization at the pivot energy and $\gamma$ is the
power-law index. The pivot energy of 26 TeV is calculated such that it
minimizes the correlation between the flux normalization and spectral index.
We summarize the best-fit parameters of HAWC J1746-2856 in Table 1 and include
both statistical and systematic uncertainties. The latter account for the
contribution of four non-negligible independent systematic uncertainties that
were identified in the previous energy-dependent study of the Crab (Abeysekara
et al., 2019) and are estimated by producing instrument response functions
with different detector configurations to investigate any potential mis-
modeling of the detector. The results were then compared with the standard
HAWC analysis and the uncertainties were added in quadrature. We do not model
the diffuse background emission because the flux detected by HAWC from the GC
is much higher than the expected Galactic diffuse emission (Cao et al., 2023;
Alfaro et al., 2024).

We calculated an UL on the minimum energy at 6 TeV and a lower limit on the
maximum energy at 114 TeV, both at the 68% CL. Above 100 TeV, 3474 events
passed trigger conditions for reconstruction, from which 98 events passed HAWC
gamma/hadron separation cuts. At 100 TeV, the energy resolution is 10% in
$\log_{10}(E/\text{TeV})$ (Abeysekara et al., 2019) and the hadron retention
after gamma/hadron separation cuts is $<1$ %. To count the events, we used a
circular region centered in the best-fit position (see Table 1) with the
radius set at the UL on the source extension.

Parameter estimated | Best fit | Statistical | Systematic
---|---|---|---
|  | uncertainties | uncertainties
RA (∘) | 266.28 |  $\pm$ 0.05 | $+0.09,-0.06$
Dec (∘) | $-28.94$ |  $\pm$ 0.04 | $+0.03,-0.02$
Flux norm. ($\phi$) $\times 10^{-15}$ (TeV cm2s)-1 | 1.5 |  $\pm$ 0.30 | $+0.08,-0.13$
Index ($\gamma$) | $-2.88$ |  $\pm$ 0.15 | $-0.1$
Table 1: Best-fit results for HAWC J1746-2856 with statistical and systematic
uncertainties. The spectrum is best described by a simple power law
$\mathrm{d}N/\mathrm{d}E=\phi(E/26\,\text{TeV})^{\gamma}$. See Section 3 for
details. In Galactic coordinates, the best-fit position of HAWC J1746-2856 is
$(l,b)=(0.06^{\circ}$, $0.09^{\circ})$.

In Figure 1(a), we show a significance map from the GC region obtained with
HAWC data by calculating the TS of every pixel as the ratio of the logarithm
of the likelihoods of the signal measured over the expected background (Younk
et al., 2015; Abeysekara et al., 2017a). We also include the location of HESS
J1745-290 (Sgr A∗;Aharonian et al., 2006b) and HESS J1746-285 (the Arc;
Abdalla et al., 2018), which are relevant to this study as they are inside of
the HAWC J1746-2856 extension UL radius and excluded from the diffuse emission
region used in H.E.S.S. Collaboration et al. (2016). In addition, we show the
positions of SNR G0.9+0.1 and HESS J1745-303. The reported gamma-ray flux
level of the SNR falls below the sensitivity of HAWC at this declination
(Albert et al., 2024). No significant excess is observed by HAWC at the
reported location of the SNR or HESS J1745-303. The 4.5$\sigma$ hot spot above
the SNR location is not coincident with any known gamma-ray sources, but it
aligns with a candidate open stellar cluster (Dutra et al., 2003). Although
gamma rays are observed in the vicinity of stellar clusters (Abramowski et
al., 2012; Aharonian et al., 2022; Abeysekara et al., 2021), the analysis
cannot rule out contributions from other unresolved sources.

Figure 1(b) shows the best-fit spectrum for HAWC J1746-2856 (see Table 1 for
systematic uncertainties and best-fit position) compared to the H.E.S.S.
measured spectra of Sgr A∗ and the radio Arc. Since HAWC cannot resolve these
point sources, we conservatively assume that their spectra extends and covers
the entire HAWC energy range, which is represented with dashed lines in Figure
1(b). In Figure 1(c), we show the significance map obtained after subtracting
the estimated excess of the H.E.S.S. point sources from the HAWC data. The
predicted event count is calculated by convolving a model consisting of the
reported best-fit parameters for the H.E.S.S. sources with the HAWC instrument
response function. We also include contours of carbon monosulfide (CS) line
emission—integrated from $-200$ km/s to 200 km/s—to show the spatial
correlation of the HAWC central excess with the density distribution of the
ambient dense gas (Tsuboi et al., 1999). Thus, the residual shown in Figure
1(c) is likely emission from the GC ridge diffusion and, in smaller
contribution, unresolved sources.

In Figure 1(d), we subtract the flux from the two H.E.S.S. point sources from
the HAWC best-fit spectrum (shown separately in Figure 1(b)). The error band
illustrates the combination of HAWC and H.E.S.S. uncertainties in quadrature.
We also compare our measurement with the diffuse emission flux points
estimated in H.E.S.S. Collaboration et al. (2016), where the diffuse emission
was derived within an annulus of inner radius 0.15∘—to exclude HESS
J1745-290—and outer radius 0.45∘. In that study, a sector ($\sim$66∘) of the
annulus is excluded to avoid HESS J1746-285. These excluded regions may
explain the slightly larger flux detected by HAWC, although the uncertainties
from both results are still compatible (see Figure 1(d)). Other IACTs have
also measured the diffuse emission, however, they use regions with
significantly different morphologies (Acciari et al., 2020; Adams et al.,
2021).

In summary, we have shown that the measured flux of HAWC J1746-2856 is
significantly higher than that of HESS J1745-290 and HESS J1746-285.
Therefore, even after excluding their contributions, the spectrum extends
beyond 100 TeV.

|
---|---
(a) | (b)
|
(c) | (d)

Figure 1: Galactic Center analysis results. (a) Significance map obtained
using the HAWC neural network energy estimator (on- and off-array events)
(Abeysekara et al., 2019) and the position of the three main point sources and
one extended source in the GC region as measured by H.E.S.S.. The dashed
circle outlines the extension UL at 68% CL. We also include the diffuse region
used in the H.E.S.S. analysis (H.E.S.S. Collaboration et al., 2016). (b)
Spectra of the two H.E.S.S. sources, along with the best-fit spectrum of HAWC
J1746-2856. The dashed lines for the H.E.S.S. sources show the extrapolation
of their best-fit to the HAWC energy range. The flux points are calculated for
each energy bin (Albert et al., 2024) by fixing all the fit parameters except
for the flux normalization. (c) HAWC emission after subtracting the two
H.E.S.S. point sources. We also show the density distribution contours of the
ambient gas as traced by CS (J1-0) line emission (Tsuboi et al., 1999). (d)
Original best-fit HAWC spectral energy distribution (SED) and the result after
subtracting the two H.E.S.S. point-source spectra. As a reference, we include
the diffuse emission measured by H.E.S.S. Collaboration et al. (2016), as
their region is almost spatially coincident with our model. See Section 3 for
details.

##  4 Discussion

The HAWC detection of photons with energies exceeding 100 TeV further
strengthens the hadronic-origin interpretation suggested by H.E.S.S.
Collaboration et al. (2016), where relativistic protons ($\gtrsim$1 PeV)
collide with the surrounding dense ambient gas.

In the leptonic scenario, the gamma-ray emission comes from the inverse
Compton scattering of electrons with energies $E_{e}>100$ TeV. In the GC
region, these electrons have a short lifetime, mostly due to synchrotron
radiation. Assuming a magnetic field strength of 100 $\mu$G (Crocker et al.,
2010), the cooling time is:

$t_{\text{cool}}\approx
13\left(\frac{E_{\text{e}}}{100\,\text{TeV}}\right)^{-1}\left(\frac{B}{100\,\mu\text{G}}\right)^{-2}\,\text{yr}\,,$
(3)

corresponding to a maximum distance that the electrons may travel
$c\,t_{\text{cool}}=4\,\text{pc}$, even assuming the extreme case of ballistic
movement. Such a distance is significantly smaller than the size of the CMZ,
which is hundreds of parsec. Therefore, the HAWC observation strongly
disfavors the leptonic scenario. The only way to make such a scenario work
would be to have tens of unresolved electron accelerators co-existing in the
region.

In the hadronic scenario, although $\pi_{0}$ decay is the dominant cool-down
channel (Aharonian et al., 2009; Longair, 2010), the cooling time is so much
larger than the escape time (by several orders of magnitude) that the proton-
cooling effect is negligible (Scherer et al., 2023). The escape time of
$E_{\text{p}}=1$ PeV protons can be roughly estimated as:

$t_{\text{escape}}\approx\frac{r^{2}}{2D}\approx
100\left(\frac{r}{40\,\text{pc}}\right)^{2}\left(\frac{E_{\text{p}}}{1\,\text{PeV}}\right)^{-0.3}\,\text{yr}\,,$
(4)

where $D\sim 1.2\times
10^{30}(E_{\text{p}}/100\text{TeV})^{0.3}\,\text{cm}^{2}/\text{s}$ (Strong et
al., 2007) is the diffusion coefficient in the interstellar medium (ISM) and
$r\sim 40$ pc is the radius of the diffuse emission region used in H.E.S.S.
Collaboration et al. (2016). As the magnetic field at the GC is much higher
than that of the average ISM (Crocker et al., 2010), protons are likely
confined therein for a longer time. Nonetheless, $t_{\text{escape}}$ is much
shorter than the age of the Galaxy, implying that the proton source(s) are
either very young or injecting protons into the CMZ in a recent burst.
Therefore, the only plausible explanation is that one or more sources quasi-
continuously accelerate and inject high-energy protons into the CMZ at rates
that exceed the escape time.

Finally, we estimated the gamma-ray luminosity ($L_{\gamma}(E_{\gamma}\geq
10\text{TeV})=2.24\times 10^{34}$ erg/s) by integrating the differential flux
of the HAWC central source between 10 and 114 TeV, subtracting the
contribution of H.E.S.S. point sources and assuming an 8.5 kpc distance to the
GC region. With this result, we calculated the energy density of cosmic-ray
protons using our measurement of the gamma-ray flux above 10 TeV to be:

$w_{\text{p}}(\geq 10E_{\gamma})=1.8\times
10^{-2}\left(\frac{\eta_{N}}{1.5}\right)^{-1}\left(\frac{L_{\gamma}(E_{\gamma}\geq
10\text{TeV})}{10^{34}\,\text{erg/s}}\right)\left(\frac{M}{10^{6}M_{\odot}}\right)^{-1}\,\text{eV/cm}^{3}\approx
8.1\times 10^{-3}\,\text{eV/cm}^{3}\,,$ (5)

where the CS total mass of the gas ($5\times 10^{6}M_{\odot}$) is the sum of
CS mass in the three H.E.S.S. annuli that are roughly coincident with the HAWC
region (H.E.S.S. Collaboration et al., 2016) and $\eta_{N}=1.5$ considers the
existence of nuclei heavier than hydrogen in cosmic rays and the interstellar
matter. This energy density obtained for $>$100 TeV protons is larger than the
$1\times 10^{-3}$ eV/cm3 local measurement by the Alpha Magnetic Spectrometer
(AMS; Aguilar et al., 2015; Abeysekara et al., 2021). Additionally, we
calculate the total energy budget of protons with energies $>$100 TeV:

$W_{p}\approx L_{\gamma}(E_{\gamma}\geq 10\text{TeV})\,t_{\text{pp}}\approx
3.53\times 10^{49}n^{-1}\,\text{erg}\,,$ (6)

where $t_{\text{pp}}\approx 5\times 10^{7}n^{-1}$ yr is the cooling time for
proton-proton (pp) interactions assuming the relative velocity of the
interacting protons to be equivalent to the speed of light ($c$) and an
ambient gas density of $n$, in units of cm-3. We estimated the cosmic-ray
energy density from H.E.S.S. measurements using the diffuse region shown in
Figure 1. By integrating the protons with energies between 100 TeV and 1140
TeV, we found the integral cosmic-ray density to be $\approx 2.1\times
10^{49}n^{-1}\,\text{erg}$, which is compatible with HAWC’s results. Our
interpretation is consistent with the steady proton source scenario suggested
by H.E.S.S. Collaboration et al. (2016). Therefore, we attribute the UHE gamma
rays to the freshly accelerated proton cosmic rays from the local accelerators
within the GC region, which continuously inject protons with PeV energies.

##  5 Conclusions

We report the first detection of $>$100 TeV gamma rays from the GC region with
a number of nearly 100 events. This HAWC result extends the highest energy
reported from the GC by the IACTs by more than a factor of two. The best-fit
model for seven years of HAWC data from the GC is a point source with a
simple-power-law spectrum
($\mathrm{d}N/\mathrm{d}E=\phi(E/26\,\text{TeV})^{\gamma}$), where
$\gamma=-2.88\pm 0.15_{\text{stat}}-0.1_{\text{sys}}$ and $\phi=1.5\times
10^{-15}$ (TeV cm2s)-1
$\pm\,0.3_{\text{stat}}\,^{+0.08_{\text{sys}}}{}_{-0.13_{\text{sys}}}$, with
no signs of a cutoff. After subtracting the small contribution of HESS
J1745-290 and HESS J1746-285 from the HAWC best-fit spectrum, the remaining
flux—likely from the Galactic ridge diffuse emission—maintains the power-law
shape, extending to at least 114 TeV. Extending the power-law spectrum to
these energies reveals a PeVatron at the GC, as first suggested by H.E.S.S.
Collaboration et al. (2016), with photons up to $\sim$30 TeV. Although our
analysis does not resolve the object accelerating protons to PeV energies, we
can confirm the existence of a PeVatron at the GC. Additionally, we discuss
the possible origin of such high-energy gamma rays—using model-independent
arguments—and conclude that the hadronic mechanism and quasi-continuous
injection scenarios are preferred. Moreover, we calculate the gamma-ray
luminosity of the PeVatron and find that the cosmic-ray energy density is
above the average, which clearly suggests the presence of freshly accelerated
0.1–1 PeV protons in the GC region. Finally, we show that the total energy
budget of protons with energies $>$100 TeV calculated with HAWC data is
compatible with H.E.S.S. measurements.

Several specific sites of proton acceleration have been proposed within the
HAWC J1746-2856 emission region, in particular near the vicinity of Sgr A∗
(H.E.S.S. Collaboration et al., 2016) and within the compact star clusters,
the Arches and Quintuplet clusters (Aharonian et al., 2019), which we did not
resolve in this analysis. Recently, there has been progress in modeling the
CMZ with more realistic cosmic-ray dynamics in agreement with existing data
(Scherer et al., 2023). The next generation of experiments, such as the
Cherenkov Telescope Array (CTA; Consortium et al., 2019) and the Southern
Wide-field Gamma-ray Observatory (SWGO; Albert et al., 2019), could better
differentiate and constrain these models with improved gamma-ray observations.

We acknowledge the support from: the US National Science Foundation (NSF); the
US Department of Energy Office of High-Energy Physics; the Laboratory Directed
Research and Development (LDRD) program of Los Alamos National Laboratory;
Consejo Nacional de Ciencia y Tecnología (CONACyT), México, grants 271051,
232656, 260378, 179588, 254964, 258865, 243290, 132197, A1-S-46288,
A1-S-22784, CF-2023-I-645, cátedras 873, 1563, 341, 323, Red HAWC, México;
DGAPA-UNAM grants IG101323, IN111716-3, IN111419, IA102019, IN106521,
IN114924, IN110521 , IN102223; VIEP-BUAP; PIFI 2012, 2013, PROFOCIE 2014,
2015; the University of Wisconsin Alumni Research Foundation; the Institute of
Geophysics, Planetary Physics, and Signatures at Los Alamos National
Laboratory; Polish Science Centre grant, DEC-2017/27/B/ST9/02272; Coordinación
de la Investigación Científica de la Universidad Michoacana; Royal Society -
Newton Advanced Fellowship 180385; Generalitat Valenciana, grant
CIDEGENT/2018/034; The Program Management Unit for Human Resources &
Institutional Development, Research and Innovation, NXPO (grant number
B16F630069); Coordinación General Académica e Innovación (CGAI-UdeG), PRODEP-
SEP UDG-CA-499; Institute of Cosmic Ray Research (ICRR), University of Tokyo.
H.F. acknowledges support by NASA under award number 80GSFC21M0002; National
Research Foundation of Korea RS-2023-00280210. We also acknowledge the
significant contributions over many years of Stefan Westerhoff, Gaurang Yodh
and Arnulfo Zepeda Domínguez, all deceased members of the HAWC collaboration.
Thanks to Scott Delay, Luciano Díaz and Eduardo Murrieta for technical
support. S. Yun-Cárcamo analyzed the data, performed the maximum likelihood
analysis and prepared the original manuscript. S. Yun-Cárcamo and D. Huang
carried out the discussion section calculations. R. Babu, K. L. Fan, and D.
Huang helped with analysis tools and interpretation of results. The full HAWC
Collaboration has contributed through the construction, calibration, and
operation of the detector, the development and maintenance of reconstruction
and analysis software, and vetting of the analysis presented in this
manuscript. All authors have reviewed, discussed, and commented on the results
and the manuscript.

## References

  * Abdalla et al. (2018) Abdalla, H., Abramowski, A., Aharonian, F., et al. 2018, Astronomy & Astrophysics, 612, A9
  * Abe et al. (2023) Abe, K., Abe, S., Aguasca-Cabot, A., et al. 2023, PoS, ICRC2023, 574, doi: 10.22323/1.444.0574
  * Abeysekara et al. (2017a) Abeysekara, A., Albert, A., Alfaro, R., et al. 2017a, The Astrophysical Journal, 843, 39
  * Abeysekara et al. (2019) —. 2019, The Astrophysical Journal, 881, 134
  * Abeysekara et al. (2021) —. 2021, Nature Astronomy, 5, 465
  * Abeysekara et al. (2023) —. 2023, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 1052, 168253, doi: https://doi.org/10.1016/j.nima.2023.168253
  * Abeysekara et al. (2017b) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017b, The Astrophysical Journal, 843, 40, doi: 10.3847/1538-4357/aa7556
  * Abeysekara et al. (2022) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2022, in 37th International Cosmic Ray Conference.
  * Abramowski et al. (2012) Abramowski, A., Acero, F., Aharonian, F., et al. 2012, Astronomy & Astrophysics, 537, A114
  * Acciari et al. (2020) Acciari, V. A., Ansoldi, S., Antonelli, L., et al. 2020, Astronomy & astrophysics, 642, A190
  * Adams et al. (2021) Adams, C., Benbow, W., Brill, A., et al. 2021, The Astrophysical Journal, 913, 115
  * Aguilar et al. (2015) Aguilar, M., Aisa, D., Alpat, B., et al. 2015, Physical review letters, 114, 171103
  * Aharonian et al. (2019) Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nature astronomy, 3, 561
  * Aharonian et al. (2006a) Aharonian, F., Akhperjanian, A., Bazer-Bachi, A., et al. 2006a, The Astrophysical Journal, 636, 777
  * Aharonian et al. (2006b) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006b, Nature, 439, 695–698, doi: 10.1038/nature04467
  * Aharonian et al. (2009) Aharonian, F., Akhperjanian, A., Anton, G., et al. 2009, Astronomy & Astrophysics, 503, 817
  * Aharonian et al. (2022) Aharonian, F., Ashkar, H., Backes, M., et al. 2022, Astronomy & Astrophysics, 666, A124
  * Albert et al. (2024) Albert, A. ., Alfaro, R., Alvarez, C., et al. 2024, Performance of the HAWC Observatory and TeV Gamma-Ray Measurements of the Crab Nebula with Improved Extensive Air Shower Reconstruction Algorithms.  https://arxiv.org/abs/2405.06050
  * Albert et al. (2019) Albert, A., Alfaro, R., Ashkar, H., et al. 2019, arXiv e-prints, arXiv:1902.08429, doi: 10.48550/arXiv.1902.08429
  * Alfaro et al. (2024) Alfaro, R., Alvarez, C., Arteaga-Velázquez, J. C., et al. 2024, The Astrophysical Journal, 961, 104, doi: 10.3847/1538-4357/ad00b6
  * Amato (2014) Amato, E. 2014, International Journal of Modern Physics D, 23, 1430013
  * Blasi (2013) Blasi, P. 2013, The Astronomy and Astrophysics Review, 21, 1
  * Cao et al. (2023) Cao, Z., Aharonian, F., An, Q., et al. 2023, Phys. Rev. Lett., 131, 151001, doi: 10.1103/PhysRevLett.131.151001
  * Consortium et al. (2019) Consortium, T. C., et al. 2019, arXiv preprint arXiv:1904.12196
  * Cristofari (2021) Cristofari, P. 2021, Universe, 7, 324, doi: 10.3390/universe7090324
  * Crocker et al. (2010) Crocker, R. M., Jones, D. I., Melia, F., Ott, J., & Protheroe, R. J. 2010, Nature, 463, 65
  * de Oña Wilhelmi et al. (2024) de Oña Wilhelmi, E., López-Coto, R., Aharonian, F., et al. 2024, Nature Astronomy, 8, 425, doi: 10.1038/s41550-024-02224-9
  * Dutra et al. (2003) Dutra, C. M., Ortolani, S., Bica, E., et al. 2003, Astronomy & Astrophysics, 408, 127
  * Fang & Halzen (2024) Fang, K., & Halzen, F. 2024, arXiv e-prints, arXiv:2404.15944, doi: 10.48550/arXiv.2404.15944
  * Gabici et al. (2019) Gabici, S., Evoli, C., Gaggero, D., et al. 2019, International Journal of Modern Physics D, 28, 1930022
  * H.E.S.S. Collaboration et al. (2016) H.E.S.S. Collaboration, Abramowski, A., Aharonian, F., et al. 2016, Nature, 531, 476, doi: 10.1038/nature17147
  * Longair (2010) Longair, M. S. 2010, High energy astrophysics (Cambridge university press)
  * Scherer et al. (2023) Scherer, A., Cuadra, J., & Bauer, F. E. 2023, Astronomy & Astrophysics, 679, A114
  * Strong et al. (2007) Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Annu. Rev. Nucl. Part. Sci., 57, 285
  * Sudoh & Beacom (2023) Sudoh, T., & Beacom, J. F. 2023, Phys. Rev. D, 107, 043002, doi: 10.1103/PhysRevD.107.043002
  * Tsuboi et al. (1999) Tsuboi, M., Handa, T., & Ukita, N. 1999, The Astrophysical Journal Supplement Series, 120, 1
  * Vianello et al. (2016) Vianello, G., Lauer, R., Younk, P., et al. 2016, PoS, ICRC2015, 1042, doi: 10.22323/1.236.1042
  * Wilks (1938) Wilks, S. S. 1938, The annals of mathematical statistics, 9, 60
  * Younk et al. (2015) Younk, P. W., Lauer, R. J., Vianello, G., et al. 2015, arXiv preprint arXiv:1508.07479
  * Yusef-Zadeh et al. (2004) Yusef-Zadeh, F., Hewitt, J., & Cotton, W. 2004, The Astrophysical Journal Supplement Series, 155, 421
  * Yusef-Zadeh & Morris (1987) Yusef-Zadeh, F., & Morris, M. 1987, Astrophysical Journal, Part 1 (ISSN 0004-637X), vol. 322, Nov. 15, 1987, p. 721-728., 322, 721


††thanks: Corresponding author<EMAIL_ADDRESS>
# CaTe: a new topological node-line and Dirac semimetal

Yongping Du1, Feng Tang1, Di Wang1, Li Sheng1,2, Er-jun Kan3, Chun-Gang Duan4,
Sergey Y. Savrasov5 and Xiangang Wan1,2 1National Laboratory of Solid State
Microstructures and Department of Physics, Nanjing University, Nanjing 210093,
China
2Collaborative Innovation Center of Advanced Microstructures, Nanjing
University, Nanjing 210093, China
3Key Laboratory of Soft Chemistry and Functional Materials (Ministry of
Education), and Department of Applied Physics, Nanjing University of Science
and Technology, Nanjing, Jiangsu 210094, P. R. China.
4Key Laboratory of Polar Materials and Devices, Ministry of Education, East
China Normal University, Shanghai 200062, China
5Department of Physics, University of California, Davis, One Shields Avenue,
Davis, California 95616, USA

###### Abstract

Topological semimetals recently stimulate intense research activities.
Combining first-principles calculations and effective model analysis, we
predict that CaTe is topological node-line semimetal when spin-orbit coupling
(SOC) is ignored. We also obtain the nearly flat surface state which has the
drumhead characteristic. When SOC is included, three node lines evolve into a
pair of Dirac points along the $M-R$ line. These Dirac points are robust and
protected by $C_{4}$ rotation symmetry. Once this crystal symmetry is broken,
the Dirac points will be eliminated, and the system becomes a strong
topological insulator.

###### pacs:

71.20.-b, 73.20.-r, 71.55.Ak, 73.43.-f

##  I Introduction

Topological insulator (TI) has attracted broad interest in recent yearsRE-1 ;
RE-2 . The unique property of TI is that the bulk state has a energy gap while
the surface state is gapless. The topological property also have been proposed
for three dimensional (3D) semimetalWeyl-1 ; Weyl-2 ; RE-3 ; Weyl-response ;
DSM-rev . Up to now, three kinds of topological semimetal have been
discovered, i.e., 3D Dirac semimetal (DSM)DSM-rev ; DSM-1 ; Na3Bi ; Cd3As2 ;
BaAgBi ; Cava , Weyl semimetal (WSM)Weyl-1 ; Weyl-2 ; HgCr2Se4 ; balets  and
node-line semimetal (NLS)NSL-1 ; NSL-2 ; NSL-3 ; NSL-4 . The 3D DSM has four-
fold degeneracy point formed by two double degeneracy linear band crossing.
Combining the crystal symmetry and time reversal symmetry, 3D DSM can be
robust against external perturbations. Based on band structural calculation,
several materials have been proposed to be 3D DSMDSM-rev ; DSM-1 ; Na3Bi ;
Cd3As2 ; BaAgBi ; Cava  and some of them already been confirmed by
experimentsDSMExp-1 ; DSMExp-3 ; DSMExp-4 ; DSMExp-5 . If one breaks time
reversal symmetryWeyl-1 ; HgCr2Se4 ; balets ; axion  or inversion
symmetryWeyl-inver-1 ; Weyl-inver-2 ; TaAs ; TaAs-hasan , the double
degeneracy bands will split, consequently the 3D DSM evolves into WSM. Very
recently, the predictions about WSM in TaAs family TaAs ; TaAs-hasan  had been
confirmed experimentally TaAs-exp-1 ; TaAs-exp-2 ; TaAs-exp-3 ; TaAs-exp-4 .

Unlike DSM and WSM whose band crossing points distribute at separate $k$
points in the Brillouin zone (BZ), for the NLS, the crossing points around the
Fermi level form a closed loop. Several compounds had been proposed as NLS
included MTCNSL-2 , Bernal graphiteNSL-mat-2 , hyperhoneycomb latticesNSL-
mat-3  and antiperovskite Cu3PdNNSL-3 ; NSL-4  and Cu3NZnNSL-4 . When SOC is
neglected, for the system with band inversion, time reversal symmetry together
with inversion symmetry or mirror symmetry will guarantee node line in 3D
BZNSL-2 ; NSL-3 ; NSL-4 ; TaAs ; NSL-mat-4 ; NSL-mat-5 . Same with TI and WSM,
NLS also has an characteristic surface state, namely, _drumhead_ like
stateNSL-1 ; NSL-2 ; NSL-3 ; NSL-4 . Such 2D flat band surface state may
become a route to achieve high temperature superconductivityNSL-app-1 ; NSL-
app-2 .

In this article, based on first-principles calculations and effective model
analysis, we propose that CaTe in CsCl-type structure is a NLS with drumhead
like surface flat bands when SOC is ignored. As shown in Fig. 1(b), around the
$M$ point, there are three node-line rings, which is perpendicular to each
other. When SOC is included, these three node-line rings evolve into two Dirac
points along the $M-R$ line. The Dirac points are robust and protected by the
$C_{4}$ rotational symmetry. If the $C_{4}$ symmetry is broken, the system
becomes a strong topological insulator with $Z_{2}$ indices (1;000).

Figure 1: (color online). (a) Crystal structure of CaTe in CsCl type phase.
(b)The 3D BZ and projected (001) two dimensional (2D) BZ of CaTe. Three dash
circles are the scheme of the three line nodes around the $M$ point. The blue
circle is parallel to $k_{z}=0$ plane, red circle is parallel to $k_{x}=0$
plane and the green circle is parallel to $k_{y}=0$ plane.

##  II Crystal structure and method

As one member of the alkaline-earth chalcogenides, CaTe have attracted
tremendous interests because of its technological applications ranging from
catalysis to luminescenceCrystal-1 ; Crystal-2 ; Crystal-3 ; Crystal-4 ;
Crystal-5 . CaTe undergoes a phase transition from NaCl-type structure at
ambient conditions to CsCl-type structure at hydrostatic pressure about 33
GPaCrystal-1 ; Crystal-2 . The structure of CaTe in CsCl-type is shown in
Fig.1(a). The space group of this phase is $Pm\overline{3}m$ (NO. 221). The
electronic band structure calculations have been carried out using the full
potential linearized augmented plane wave method as implemented in WIEN2K
package WIEN2K . To obtain accurate band inversion strength and band order,
the modified Becke-Johnson exchange potential together with local-density
approximation for the correlation potential (MBJLDA) has been applied mBJ .
The plane-wave cutoff parameter RMTKmax is set to 7 and a $16\times 16\times
16$ mesh was used for the BZ integral. The SOC interaction is included by
using the second-order variational procedure.

##  III Electronic structure

Firstly, we calculate the band structure of CaTe and show the result without
SOC in Fig. 2(a). By checking the wave functions, we find that the valence
bands and conduction bands are mainly contributed by 5$p_{z}$ (blue) state of
Te and 3$d_{z^{2}}$ (red) state of Ca, respectively, as shown in Fig. 2(a).
The band inversion happened at $M$ point where the energy of Te-5$p_{z}$ state
is higher than the energy of Ca-3$d_{z^{2}}$ state by about 0.75 eV.
Interestingly, this kind of band inversion is not caused by the SOC, which is
different from most topological materialsRE-1 ; RE-2 . We calculate the
electronic structure of CaTe by applying tensile strain to check the origin of
band inversion at $M$ point. The energy difference between Te-5$p_{z}$ state
and Ca-3$d_{z^{2}}$ state decreases as increasing the tensile strain. We find
that when a$\geq$1.05a0, the band inversion at M point disappear. With the
time reversal symmetry and inversion symmetry, the band inversion results in
node lines as proved in Ref. NSL-2 .

Figure 2: (color online). (a) Electronic structure of CaTe without SOC. The
weights of Te-5$p_{z}$ (Ca-3$d_{z^{2}}$)state is proportional to the width of
blue (red) curves. (b) Electronic structure of CaTe with SOC. Along
$\Gamma-M$, $M-X$ and $M-M^{\prime}$, a small gap is opened. The Dirac point
at the $M-R$ line is protected by $C_{4}$ rotational symmetry. The inset is
the detailed structure around the M point. (see main text for detailed
description).

The effective Hamiltonian of node line around $M$ point can be established by
using the $\overset{\rightharpoonup}{k}\cdot$ $\overset{\rightharpoonup}{p}$
method. Considering the crystal symmetry at $M$ point and time reversal
symmetry, the effective Hamiltonian can be written as following form:

$H(\overset{\rightharpoonup}{k})=g_{0}(\overset{\rightharpoonup}{k})\sigma_{0}+g_{x}(\overset{\rightharpoonup}{k})\sigma_{x}+g_{z}(\overset{\rightharpoonup}{k})\sigma_{z},$

where the $\sigma_{x}\ $and $\sigma_{z}$ are Pauli matrices, $\sigma_{0}$ is
unit matrix.
$g_{0}(\overset{\rightharpoonup}{k})=M_{0}-B_{0}(k_{x}^{2}+k_{y}^{2})-C_{0}k_{z}^{2}$,
$g_{x}(\overset{\rightharpoonup}{k})=Ak_{x}k_{y}k_{z}$,
$g_{z}(\overset{\rightharpoonup}{k})=M_{z}-B_{z}(k_{x}^{2}+k_{y}^{2})-C_{z}k_{z}^{2}$.
This system has both of time reversal symmetry and inversion symmetry, thus,
the component of $\sigma_{y}$ is zeroNSL-2 . The eigenvalues of this $2\times
2$ effective Hamiltonian are
$E(\overset{\rightharpoonup}{k})=g_{0}(\overset{\rightharpoonup}{k})\pm\sqrt{g_{x}^{2}(\overset{\rightharpoonup}{k})+g_{z}^{2}(\overset{\rightharpoonup}{k})}$.
When $g_{x}(\overset{\rightharpoonup}{k})=0$ and
$g_{z}(\overset{\rightharpoonup}{k})=0$, the nodal line will emergent. It can
be easily checked that the equation $g_{z}(\overset{\rightharpoonup}{k})=0$
has solution only when $M_{z}B_{z}>0$ and $M_{z}C_{z}>0$. And $M_{z}B_{z}>0$
and $M_{z}C_{z}>0$ are also the condition of band inversion. On the other
hand, $g_{x}(\overset{\rightharpoonup}{k})=0$ confine the node lines in three
mutually perpendicular planes (namely, $k_{x}$=0 plane, $k_{y}$=0 plane and
$k_{z}$=0 plane ) as illustrated in Fig. 1(b). Due to the fact that
$g_{0}(\overset{\rightharpoonup}{k})$ does not equal to zero which breaks the
electron-hole symmetry, consequently, the node lines have finite energy
dispersion.

When SOC is considered, three node lines evolve into two Dirac points at $M-R$
line as shown in Fig. 2(b). At $M$ point, the two states near Fermi level
belong to irreducible representation $\Gamma_{7}^{-}$ and $\Gamma_{6}^{+}$,
respectively. While along the $M-X$ line, two bands have the same irreducible
representation $\Gamma_{5}\ $as shown in Fig. 2(b), thus they can hybridize
with each other and open a small gap (about 50meV). For both $\Gamma-M$ line
and $M-M^{\prime}$ line, the two bands around Fermi level are also belong to
the same irreducible representation, thus there is not band crossing along
$\Gamma-M$ line and $M-M^{\prime}$ line. Since the band splitting is
determined by the SOC, thus one can achieve the NLS by doping the lighter
atoms such as Se, S.

Along the $M-R$ line, which reserve the $C_{4}$ rotation symmetry, two states
with $\Gamma_{7}^{-}$ and $\Gamma_{6}^{+}$ at $M$ point evolve into
$\Gamma_{7}$ and $\Gamma_{6}$, thus the hybridization between these two bands
is forbidden, there is a Dirac point as shown in Fig.2(b). When the $C_{4}$
rotational symmetry is broken, like by strain effect, the band crossing point
will disappear, and this 3D DSM will become a strong TI with topological
indices $Z_{2}$ to be(1;0,0,0).

To understand the band inversion at $M$ point and the topological property of
this system, we derive a low energy effective Hamiltonian at $M$ point based
on the projection-operator methodBaAgBi . $M$ point has $D_{4h}$ symmetry and
also time reversal symmetry. As discussed above, at $M$ point,
$\Gamma_{7}^{-}$ symmetry state has angular momentum $j_{z}=\pm 3/2$ and
$\Gamma_{6}^{+}$ symmetry state has angular momentum $j_{z}=\pm 1/2$.
Therefore using the basis of
($|j_{z}=-\frac{1}{2}\rangle_{d},|j_{z}=+\frac{1}{2}\rangle_{d},|j_{z}=-\frac{3}{2}\rangle_{p},|j_{z}=+\frac{3}{2}\rangle_{p}$,
the effective Hamiltonian around $M$ point can be written as (see APPENDIX for
detail.):

$H_{eff}=\left(\begin{matrix}M_{1}(\overset{\rightharpoonup}{k})&0&Ak_{+}+B(\overset{\rightharpoonup}{k})&D(\overset{\rightharpoonup}{k})\\\
0&M_{1}(\overset{\rightharpoonup}{k})&D^{\ast}(\overset{\rightharpoonup}{k})&-Ak_{-}-B^{\ast}(\overset{\rightharpoonup}{k})\\\
Ak_{-}+B^{\ast}(\overset{\rightharpoonup}{k})&D(\overset{\rightharpoonup}{k})&M_{2}(\overset{\rightharpoonup}{k})&0\\\
D^{\ast}(\overset{\rightharpoonup}{k})&-Ak_{+}-B(\overset{\rightharpoonup}{k})&0&M_{2}(\overset{\rightharpoonup}{k}),\end{matrix}\right)$

where
$M_{1}(\overset{\rightharpoonup}{k})=M_{10}+M_{11}(k_{x}^{2}+k_{y}^{2})+M_{12}k_{z}^{2}$,
$M_{2}(\overset{\rightharpoonup}{k})=M_{20}+M_{21}(k_{x}^{2}+k_{y}^{2})+M_{22}k_{z}^{2}$,
$B(\overset{\rightharpoonup}{k})=B_{1}k_{+}k_{z}^{2}+B_{2}(k_{x}^{3}+ik_{y}^{3})+iB_{3}k_{-}k_{x}k_{y}$,
$D(\overset{\rightharpoonup}{k})=D_{1}(k_{x}^{2}-k_{y}^{2})k_{z}+iD_{2}k_{x}k_{y}k_{z}$,
$k_{\pm}=k_{x}\pm ik_{y}$. Along the $k_{z}$ axis (where $k_{x}$=0, $k_{y}$=0)
the effective Hamiltonian is diagonal, and the eigenvalues are
$E(\overset{\rightharpoonup}{k})=M_{1}(\overset{\rightharpoonup}{k})$ and
$E(\overset{\rightharpoonup}{k})=M_{2}(\overset{\rightharpoonup}{k})$. As
mentioned above, the Dirac point is on the $M-R$ line, thus it is interesting
to discuss the effective model along this line. Since there is the band
inversion between $\left|j_{z}=\pm\frac{1}{2}\right\rangle_{d}$ and
$\left|j_{z}=\pm\frac{3}{2}\right\rangle_{p}$ at $M$ point, it is easy to
obtain that $M_{10}<M_{20}$, $M_{22}<0<M_{12}$, and the Dirac points locate at
$\overset{\rightharpoonup}{k_{c}}=(\frac{\pi}{a},\frac{\pi}{a},k_{zc}$=$\pm\sqrt{\frac{M_{20}-M_{10}}{M_{12}-M_{22}}})$.
Neglecting the high-order terms,
$E(\overset{\rightharpoonup}{k_{c}}+\overset{\rightharpoonup}{\delta k})$ can
be expressed as $(M_{12}+M_{22})k_{zc}\delta
k_{z}\pm\sqrt{(M_{12}-M_{22})^{2}k_{zc}^{2}\delta k_{z}^{2}+A^{2}(\delta
k_{x}^{2}+\delta k_{y}^{2})}$, where $\delta k_{x,y,z}$ are small displacement
from $\overset{\rightharpoonup}{k_{c}}$. In the vicinity of
$\overset{\rightharpoonup}{k_{c}}$, the band dispersion is a linear, thus our
effective Hamiltonian is nothing but 3D massless Dirac fermions.

Figure 3: (color online). The surface state of (001) surface of CaTe. (a) The
surface state is calculated without SOC. The flat surface state around Fermi
level is denoted by red color. The inset is the detailed band structure around
the $\overline{M}$ point. (b) The surface state is calculated with SOC. The
red dots are the projected bulk Dirac nodes. The red lines between the bulk
gap at $\overline{X}$ point are surface states and the blue circle denotes the
surface Dirac cones. The inset is the detailed band structure around the
$\overline{X}$ point.

The band inversion at $M$ point and the Dirac nodes in CaTe suggest the
existence of topological nontrivial surface state. To study the surface states
in CaTe we use a 200-unit-cells-thick (001) slab with top (bottom) surface
terminated by Ca (Te) atoms. The surface state is then calculated by using the
tight-binding method. The hopping parameters are determined from a maximally
localized Wannier functions (MLWFs)MLWFs , which are projected from the Bloch
state derived from first-principles calculations.

Fig. 3(a)/(b) shows the surface state of CaTe (001) surface without/with SOC,
respectively. When SOC is ignored, the system is a NLS, and possess nearly
flat surface band around the Fermi energy. As shown in Fig. 3(a), our
numerical results find that the nearly flat surface ”drumhead” state appears
in the interiors of the projected nodal line rings on the (001) surface around
the $\overline{M}$ point. Since the slab we used has two surfaces, there are
two surface states as shown in the red lines in the Fig. 3(a). The particle-
hole symmetry is broken by nonzero term $g_{0}(\overset{\rightharpoonup}{k})$,
thus these two surface bands are not perfect flat with about 70 meV bandwidth.
This type of 2D flat bands are proposed as a novel route to achieve high
temperature superconductivityNSL-app-1 ; NSL-app-2 .

When the SOC is included, three node lines are gapped out and become a pair of
Dirac points along the $M-R$ line, thus the NLS become a 3D DSM. There is bulk
Dirac node projected on $\overline{M}$ point (the red dot) as shown in Fig.
3(b). Along the $\overline{M}-\overline{X}$, there is also a projected bulk
Dirac node, which locate near the $\overline{X}$ point denoted by red dots.
Fig. 3(b) is clearly shows that the gapped bulk states along the
$\overline{\Gamma}-\overline{X}$ direction and the existence of surface Dirac
cones (in the blue circle) due to the topologically nontrivial $Z_{2}$
indices, like the same case in Na3BiNa3Bi  and Cd3As2Cd3As2 .

##  IV Conclusion

In summary, based on first-principles calculations and effective model
analysis, we suggest that CaTe in CsCl-type structure is a NLS when SOC is
ignored. There are three node-line rings which are perpendicular to each other
around the $M$ point. With band inversion at $M$ point, this NLS is protected
by the time reversal symmetry and inversion symmetry. When the SOC is
included, three node-line rings become a pair of Dirac points. These Dirac
nodes are robust and protected by the $C_{4}$ crystal symmetry and the system
become a DSM.

##  V Acknowledgement

The work was supported by the National Key Project for Basic Research of China
(Grant No. 2014CB921104), NSFC under Grants No. 11525417 and 11374147. The
project is also funded by Priority Academic Program Development of Jiangsu
Higher Education Institutions. S.S. were supported by NSF DMR (Grant No.
1411336). Y.D. is supported by the program B for Outstanding PhD candidate of
Nanjing University.

##  VI APPENDIX

The conduction and valence bands of CaTe at $M$ point are mainly contributed
by four states: $|j_{z}=-\frac{1}{2}\rangle_{d}$,
$|j_{z}=+\frac{1}{2}\rangle_{d}$, $|j_{z}=-\frac{3}{2}\rangle_{p}$ and
$|j_{z}=+\frac{3}{2}\rangle_{p}$, we thus use these states as the basis to
build the effective model Hamiltonian at the $M$ point of BZ. As a $4\times 4$
Hermitian matrix, the effective Hamiltonian can be expended as
$H=\underset{i=1}{\overset{16}{\sum}}f_{i}(\overset{\rightharpoonup}{k})\Gamma_{i}$,
where $f_{i}(\overset{\rightharpoonup}{k})$ are function of momentum k.
$\Gamma_{i}$ are Dirac matrices, which can be written as the direct product of
$\sigma_{i}$ and $\tau_{j}$ ($\sigma_{i=1,2,3,4}$, $\tau_{j=1,2,3,4}$ are unit
matrix $\sigma_{0}$ and three Pauli matrices $\sigma_{x}$, $\sigma_{y}$ and
$\sigma_{z}$).

Under the operation of crystal symmetry and time reversal symmetry, the
Hamiltonian should be invariant. This requires the function
$f_{i}(\overset{\rightharpoonup}{k})$ and the associated $\Gamma_{i}$ matrices
belong to the same irreducible representation. Thus the key problem is to
determine the irreducible representation for
$f_{i}(\overset{\rightharpoonup}{k})$ and $\Gamma_{i}$ matrices, which can be
done by the projection-operator method.

$\Gamma$ matrices | representation | T
---|---|---
$\Gamma_{1}$,$\Gamma_{13}$ | $R_{1}$ | +
$\Gamma_{4},\Gamma_{16}$ | $R_{2}$ | -
$\\{\Gamma_{2},\Gamma_{3}\\},\\{\Gamma_{14},\Gamma_{15}\\}$ | $R_{5}$ | -
$\Gamma_{7}$ | $R_{10}$ | -
$\Gamma_{11}$ | $R_{10}$ | +
$\Gamma_{6}$ | $R_{11}$ | -
$\Gamma_{10}$ | $R_{11}$ | +
$\\{\Gamma_{8},\Gamma_{9}\\}$ | $R_{12}$ | -
$\\{\Gamma_{5},\Gamma_{12}\\}$ | $R_{12}$ | +
$d(k)$ | representation | T
$C$, $k_{z}^{2}$, $k_{x}^{2}+k_{y}^{2}$ | $R_{1}$ | +
$k_{x}^{2}-k_{y}^{2}$ | $R_{3}$ | +
$k_{x}k_{y}$ | $R_{4}$ | +
$\\{k_{x}k_{z}$, $k_{y}k_{z}\\}$ | $R_{5}$ | +
$k_{z},k_{z}^{3},(k_{x}^{2}+k_{y}^{2})k_{z}$ | $R_{9}$ | -
$k_{x}k_{y}k_{z}$ | $R_{10}$ | -
$\\{k_{x}^{2}-k_{y}^{2}\\}k_{z}$ | $R_{11}$ | -
$(k_{x},k_{y}),(k_{x}^{3},k_{y}^{3}),(k_{x}^{2}k_{y},k_{y}^{2}k_{x}),(k_{x}k_{z}^{2},k_{y}k_{z}^{2})$ | $R_{12}$ | -
Table 1: The character table of Dirac $\Gamma$ matrices and the polynomials of
the momentum $k$ for CaTe at M point.

Because the SOC is included, we use the double space group. Under the
projection-operator method, we present the irreducible representation of Dirac
$\Gamma$ matrices and polynomials of $\overset{\rightharpoonup}{k}$, and their
transformation under time reversal in Table I. With the Table I, the effective
model Hamiltonian of CaTe at $M$ point can be easily expressed as:
$H=f_{1}(\overset{\rightharpoonup}{k})\Gamma_{1}+f_{13}(\overset{\rightharpoonup}{k})\Gamma_{13}+f_{8}(\overset{\rightharpoonup}{k})\Gamma_{8}-f_{9}(\overset{\rightharpoonup}{k})\Gamma_{9}+D_{1}(k_{x}^{2}-k_{y}^{2})k_{z}\Gamma_{6}-D_{2}k_{x}k_{y}k_{z}\Gamma_{7}$,
where
$f_{1}(\overset{\rightharpoonup}{k})=C_{1}+m_{1}(k_{x}^{2}+k_{y}^{2})+n_{1}k_{z}^{2}$,
$f_{13}(\overset{\rightharpoonup}{k})=C_{13}+m_{13}(k_{x}^{2}+k_{y}^{2})+n_{13}k_{z}^{2}$,
$f_{8}(\overset{\rightharpoonup}{k})=Ak_{x}+B_{1}k_{x}^{3}+B_{2}k_{x}k_{z}^{2}+B_{3}k_{y}^{2}k_{x}$,
$f_{9}(\overset{\rightharpoonup}{k})=Ak_{y}+B_{1}k_{y}^{3}+B_{2}k_{y}k_{z}^{2}+B_{3}k_{x}^{2}k_{y}$.
$\Gamma_{1}=\sigma_{0}\otimes\tau_{0}$,
$\Gamma_{13}=\sigma_{3}\otimes\tau_{0}$,
$\Gamma_{8}=\sigma_{1}\otimes\tau_{3}$,
$\Gamma_{9}=\sigma_{2}\otimes\tau_{0}$,
$\Gamma_{6}=\sigma_{1}\otimes\tau_{1}$,
$\Gamma_{7}=\sigma_{1}\otimes\tau_{2}$. Compare with the effective
Hamiltonian, we have $M_{10}(M_{20})=C_{1}\pm C_{13}$,
$M_{11}(M_{21})=m_{1}\pm m_{13}$, $M_{12}(M_{22})=n_{1}\pm n_{13}$.

## References

  * (1) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (2) X. L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (3) X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Phys. Rev. B 83, 205101 (2011).
  * (4) L. Balents, Physics 4, 36 (2011).
  * (5) Y. Ando, J. Phys. Soc. of Jpn. 82, 102001 (2013).
  * (6) K.-Y. Yang, Y.-M. Lu, and Y. Ran, Phys. Rev. B 84, 075129 (2011).
  * (7) T. O. Wehling, A. M. Black-Schaffer, and A. V. Balatsky, Adv. Phys. 63, 1 (2014).
  * (8) G. Xu, H. M. Weng, Z. Wang, X. Dai, and Z. Fang, Phys. Rev. Lett. 107, 186806 (2011).
  * (9) A. A. Burkov and Leon Balents, Phys. Rev. Lett. 107, 127205 (2011)
  * (10) S. M. Young, S. Zaheer, J. C. Y. Teo, C. L. Kane, E. J. Mele, and A. M. Rappe, Phys. Rrev. Lett. 108, 140405 (2012).
  * (11) Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu, H. M. Weng, X. Dai, and Z. Fang, Phys. Rev. B 85, 195320 (2012).
  * (12) Z. Wang, H. M. Weng, Q. Wu, X. Dai, and Z. Fang, Phys. Rev. B 88, 125427 (2013).
  * (13) Y. Du, B. Wan, D. Wang, L. Sheng, C. Duan and X. Wan, Sci. Rep. 5, 14423 (2015).
  * (14) Q. D. Gibson, L. M. Schoop, L. Muechler, L. S. Xie, M. Hirschberger, N. P. Ong, R. Car, and R. J. Cava Phys. Rev. B 91, 205128 (2015).
  * (15) A. A. Burkov, M.D. Hook, and L. Balents, Phys. Rev. B 84, 235126 (2011)
  * (16) H. Weng, Y. Liang, Q. Xu, R. Yu, Z. Fang, X. Dai, and Y. Kawazoe, Phys. Rev. B 92, 045108 (2015)
  * (17) R. Yu, H. Weng, Z. Fang, X. Dai and X. Hu, Phys. Rev. Lett. 115, 036807 (2015)
  * (18) Y. Kim, B. J. Wieder, C. L. Kane, and A. M. Rappe, Phys. Rev. Lett. 115, 036806 (2015)
  * (19) S. Jeon, B. B. Zhou, A. Gyenis, B. E. Feldman, I. Kimchi, A. C. Potter, Q. D. Gibson, R. J. Cava, A. Vishwanath, and A. Yazdani, Nat. Mater. 13, 851 (2014).
  * (20) M. Neupane, S.-Y. Xu, R. Sankar, N. Alidoust, G. Bian, C. Liu, I. Belopolski, T.-R. Chang, H.-T. Jeng, H. Lin, A. Bansil, F. Chou, and M. Z. Hasan, Nat. commun. 5, 3786 (2014).
  * (21) Z. K. Liu, B. Zhou, Y. Zhang, Z. J. Wang, H. M. Weng, D. Prabhakaran, S.-K. Mo, Z. X. Shen, Z. Fang, X. Dai, Z. Hussain, and Y. L. Chen, Science 343, 864 (2014).
  * (22) Z. K. Liu, J. Jiang, B. Zhou, Z. J. Wang, Y. Zhang, H. M. Weng, D. Prabhakaran, S. K. Mo, H. Peng, P. Dudin, T. Kim, M. Hoesch, Z. Fang, X. Dai, Z. X. Shen, D. L. Feng, Z. Hussain, and Y. L. Chen, Nature Materials 13, 677 (2014).
  * (23) X. Wan, A. Vishwanath, and S. Y. Savrasov, Phys. Rev. Lett. 108, 146601 (2012).
  * (24) G. B. Halasz and L. Balents, Phys. Rev. B 85, 035103 (2012).
  * (25) J. Liu and D. Vanderbilt, Phys. Rev. B 90, 155316 (2014).
  * (26) H. Weng, C. Fang, Z. Fang, B. A. Bernevig, and X. Dai, Phys. Rev. X 5, 011029 (2015).
  * (27) S. M. Huang, S. Y. Xu, I. Belopolski, C. C. Lee, G. Chang, B. Wang, N. Alidoust, G. Bian, M. Neupane, C. Zhang, S. Jia, A. Bansil, H. Lin, and M. Z. Hasan, Nat. commun. 6, 7373 (2015).
  * (28) B. Q. Lv, H. M. Weng, B. B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, Z. Fang, X. Dai, T. Qian, and H. Ding, Phys. Rev. X 5, 031013 (2015).
  * (29) B. Q. Lv, N. Xu, H. M. Weng, J. Z. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, C. E. Matt, F. Bisti, V. N. Strocov, J. Mesot, Z. Fang, X. Dai, T. Qian, M. Shi and H. Ding, Nat. Phys. 11, 724–727 (2015).
  * (30) S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, S.-M. Huang, H. Zheng, J. Ma, D. S. Sanchez, B. Wang, A. Bansil, F. Chou, P. P. Shibayev, H. Lin, S. Jia, M. Z. Hasan, Science 349, 613 (2015).
  * (31) S.-Y. Xu, N. Alidoust, I. Belopolski, Z. Yuan, G. Bian, T.-R. Chang, H. Zheng, V. N. Strocov, D. S. Sanchez, G. Chang, C. Zhang, D. Mou, Y. Wu, L. Huang, C.-C. Lee, S.-M. Huang, B. Wang, A. Bansil, H.-T. Jeng, T. Neupert, A. Kaminski, H. Lin, S. Jia and M. Z. Hasan, Nat. Phys. 11, 748–754 (2015).
  * (32) T. T. Heikkila and G. E. Volovik, arXiv: 1504.05824v1 (2015.)
  * (33) K. Mullen, B. Uchoa, and D. T. Glatzhofer, Phys. Rev. Lett. 115, 026403 (2015).
  * (34) L. S. Xie, L. M. Schoop, E. M. Seibel, Q. D. Gibson, W. Xie, and R. J. Cava, APL Materials 3, 083602 (2015).
  * (35) M. Zeng, C. Fang, G. Chang, Y. A. Chen, T. Hsieh, A. Bansil, H. Lin, L. Fu, arXiv:1504.03492 (2015)
  * (36) T. T. Heikkila and G. E. Volovik, arXiv: 1504.05824 (2015)
  * (37) N. B. Kopnin, T. T. Heikkila, and G. E. Volovik, Phys. Rev. B 83, 220503 (2011).
  * (38) H. G. Zimmer, H. Winzen, and K. Syassen, Phys. Rev. B 32, 4066 (1985).
  * (39) H. Luo, R. G. Greene, K. Ghandehari, T. Li, and A. L. Ruoff, Phys. Rev. B 50, 16232 (1994).
  * (40) M. Dadsetani and A. Pourghazi, Phys. Rev. B 73, 195102 (2006).
  * (41) S. Boucenna, Y. Medkour, L. Louail, M. Boucenna, A. Hachemi, A. Roumili, Comput. Mater. Sci. 68, 325 (2013).
  * (42) D. Varshneya, V. Rathore, R. Kinge, R. K. Singh, J. Alloys Compd. 484, 239 (2009).
  * (43) P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, and J. Luitz, WIEN2K, An Augmented Plane Wave+ Local Orbitals Program for Calculating Crystal Properties (Karlheinz Schwarz, Technische Universitat Wien, Austria, 2001).
  * (44) F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009).
  * (45) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012)