#  Hemispheric Tectonics on super-Earth LHS 3844b

Tobias G. Meier Center for Space and Habitability, University of Bern,
Gesellschaftsstrasse 6, 3012 Bern, Switzerland  Dan J. Bower Center for Space
and Habitability, University of Bern, Gesellschaftsstrasse 6, 3012 Bern,
Switzerland  Tim Lichtenberg Atmospheric, Oceanic and Planetary Physics,
Department of Physics, University of Oxford, Parks Road, Oxford OX1 3PU,
United Kingdom  Paul J. Tackley Institute of Geophysics, Department of Earth
Sciences, ETH Zurich, Sonneggstrasse 5, Zurich 8092, Switzerland  Brice-
Olivier Demory Center for Space and Habitability, University of Bern,
Gesellschaftsstrasse 6, 3012 Bern, Switzerland

(Received November 20, 2020; Revised January 22, 2021; Accepted February 4,
2021)

###### Abstract

The tectonic regime of rocky planets fundamentally influences their long-term
evolution and cycling of volatiles between interior and atmosphere. Earth is
the only known planet with active plate tectonics, but observations of
exoplanets may deliver insights into the diversity of tectonic regimes beyond
the solar system. Observations of the thermal phase curve of super-Earth LHS
3844b reveal a solid surface and lack of a substantial atmosphere, with a
temperature contrast between the substellar and antistellar point of around
1000 K. Here, we use these constraints on the planet’s surface to constrain
the interior dynamics and tectonic regimes of LHS 3844b using numerical models
of interior flow. We investigate the style of interior convection by assessing
how upwellings and downwellings are organized and how tectonic regimes
manifest. We discover three viable convective regimes with a mobile surface:
(1) spatially uniform distribution of upwellings and downwellings, (2)
prominent downwelling on the dayside and upwellings on the nightside, and (3)
prominent downwelling on the nightside and upwellings on the dayside.
Hemispheric tectonics is observed for regimes (2) and (3) as a direct
consequence of the day-to-night temperature contrast. Such a tectonic mode is
absent in the present-day solar system and has never been inferred from
astrophysical observations of exoplanets. Our models offer distinct
predictions for volcanism and outgassing linked to the tectonic regime, which
may explain secondary features in phase curves and allow future observations
to constrain the diversity of super-Earth interiors.

exoplanets, super-Earths, planetary interior, plate tectonics

††journal: ApJL

##  1 Introduction

Plate tectonics is the unifying theory of Earth Science that explains the
geological and surface evolution of Earth for at least the past 3 Gyr. Plate
tectonics is a fundamental component of long-term (geological) cycles that
enable exchange of volatiles between the interior and atmosphere. These cycles
regulate climate and provide the necessary ingredients to nurture and sustain
life on Earth and are thus essential to understanding the habitability of
distant worlds. The discovery of plate tectonics on Earth arose from seafloor
mapping and seismology, but these techniques cannot be used to discern
tectonic regimes on rocky exoplanets. Instead, observations of the thermal
phase curve are available for select planets. Here, we link phase curve
observations to numerical models of interior flow and constrain the possible
tectonic regimes of super-Earth LHS 3844b.

Since the first thermal map of a super-Earth was constructed for super-Earth
55 Cnc e (Demory et al., 2012, 2016), other super-Earths have been targeted
with observations to constrain their thermal emission and thus constrain their
day- and nightside temperatures. Recent efforts focused on LHS 3844b have
found coincidence between the substellar point and the observed hotspot, and
inferred a dayside temperature of $1040\pm 40$ K and a nightside temperature
around 0 K (Kreidberg et al., 2019). These suggest that heat redistribution is
inefficient; therefore the planet has neither substantial melt at its surface
nor an active atmosphere. The observational data suggests the emission
originates from bare rock with a low albedo (Kreidberg et al., 2019).
Therefore, LHS 3844b has established itself as a case study for understanding
the interior dynamics and tectonics of ultra-short period super-Earths.

Whether or not plate tectonics (mobile lid convection) operates on super-
Earths has been a long-standing debate (e.g., Valencia et al., 2006; O’Neill &
Lenardic, 2007), but the debate is notably hindered by the lack of
observational constraints to supplement theoretical and numerical modeling
efforts. Tectonic regimes manifest from the interior dynamics of a planet,
specifically the style and vigor of convection in the outermost silicate shell
(the mantle). Hence tectonic regimes are the surface expression of mantle flow
that extends deep in the planet. Therefore, introducing observational
constraints into models is essential to provide new insights into the viable
tectonic regimes operating on rocky exoplanets such as LHS 3844b.

It is unknown how the strong temperature contrast imposed by stellar
irradiation on LHS 3844b influences its interior flow and hence its tectonic
regimes. Thus we utilize advanced models of interior flow coupled with thermal
phase curve observations to determine the viable convective regimes operating
in the interior of an ultra-short period super-Earth. This enables us to probe
the coupling between the surface and interior of LHS 3844b and thereby infer
observational strategies for further geological characterization.

##  2 Constraining interior dynamics with observations

LHS 3844b is 1.3 Earth radii (Vanderspek et al., 2019) and the dayside and
nightside temperatures are 1040 K and $\approx 0$ K, respectively (Kreidberg
et al., 2019). We estimate the longitudinal temperature variation by assuming
a blackbody at equilibrium and that only the dayside reradiates energy
received from the star. Since the thermal phase curve only constrains the
average longitudinal dependence of surface temperature, we construct models of
interior flow within 2D spherical annulus geometry (Hernlund & Tackley, 2008).
Mantles behave as highly viscous fluids, so we solve for Stokes flow (mass,
momentum, and energy conservation) using the mantle convection code StagYY
(Tackley, 2008). Since the mass of LHS 3844b has not been measured, we assume
its mantle is silicate rock and has the same relative thickness as Earth. We
adopt a hydrostatic reference state to provide the mantle profiles of material
properties (Tackley et al., 2013). Utilizing an Arrhenius-type viscosity law,
the mantle is modeled with an upper mantle, a perovskite layer, and a post-
perovskite layer assuming a lower bound estimate of the viscosity (Tackley et
al., 2013). The perovskite-post-perovskite interface occurs at a depth around
$1680$ km (total mantle depth is $\approx$3500 km).

Lithospheric strength is modeled by a plastic yielding criteria ($\sigma_{\rm
duct}$) to obtain self-consistent plate-like behavior (Moresi & Solomatov,
1998; Tackley, 2000). Laboratory experiments estimate the yield stress to be a
few hundred MPa (Kohlstedt et al., 1995) and numerical simulations employ a
yield stress less than $\sim 150$ MPa to obtain plate-like behavior for Earth-
like planets (Tackley, 2000). We vary the ductile yield stress between 10 MPa
and 300 MPa to model a lithosphere that is weak and strong, respectively.
Planetary mantles can be heated from below from cooling of an iron-rich
(geophysical) core and heated from within from decay of radionuclides or tidal
heating. We therefore explore two heating modes: (1) basal heating, where the
mantle is heated exclusively by the geophysical core, and (2) mixed-mode
heating, where constant mantle heat production (Earth-like at $5.2\times
10^{-12}$ W kg-1) supplements heating by the core. Kane et al. (2020) estimate
the age of the host star as $7.8\pm 1.6$ Gyr, but an estimate of the stellar
radionuclide abundances is presently not available to provide a constraint on
the radiogenic heat budget of the planet (e.g., Unterborn et al., 2015).
Therefore, we assume an Earth-like internal heat budget, in accordance with
previous modeling efforts (e.g., Kane et al., 2020).

##  3 Tectonic regimes and interior flow

Our simulations discover three mobile lid tectonic regimes, each of which is
associated with a distinct interior temperature and flow field: (1) uniform
distribution of upwellings and downwellings (Figs. 1A and 2A), (2)
downwellings on the dayside and upwellings on the nightside (Figs. 1B, D and
2B, D), and (3) downwellings on the nightside and upwellings on the dayside
(Figs. 1C and 2C). Figure 1 shows the mantle temperature for three times of
each model where the substellar point is located at 0 degrees. Figure 2 shows
evolutionary tracks that summarize the distribution of upwellings (red tracks)
and downwellings (blue tracks) with longitude and time. The time of the
corresponding temperature fields in Fig. 1 is indicated by horizontal dashed
lines in Fig. 2. The internal heating ratio (e.g., Korenaga, 2017) is zero (by
definition) for the basally heated models and between 0.7 and 0.8 for the
internally heated cases. The time-averaged core–mantle boundary (CMB) heat
flux is $\approx 70$ mW m-2 and $\approx 50$ mW m-2 for the internally heated
weak and strong lithosphere models, respectively. For the basally heated
models, the CMB heat flux is $\approx 95$ mW m-2 and $\approx 65$ mW m-2 for
the weak and strong lithosphere models, respectively.

Figure 1: Mantle temperature for tectonic regimes of LHS 3844b for (A) weak
lithosphere ($\sigma_{\rm duct}$= 10 MPa) and basal heating, (B) strong
lithosphere ($\sigma_{\rm duct}$= 300 MPa) and basal heating, (C) weak
lithosphere ($\sigma_{\rm duct}$= 10 MPa) and mixed heating, and (D) strong
lithosphere ($\sigma_{\rm duct}$= 300 MPa) and mixed heating. Upper left label
in each cross-section (A3, B3, etc.) indicates the time of the cross-section
in Figure 2. The substellar point is at 0∘ longitude and the nightside (90∘ –
270∘) is denoted by a gray background. White dashed line in the mid-mantle in
B2, C2, and D2 shows the boundary between perovskite (low pressure) and post-
perovskite (high pressure).

###  3.1 Uniform distribution of upwellings and downwellings

A uniform distribution of upwellings and downwellings is predicted if the
planet is dominantly basal heated and has a weak lithosphere (Figs. 1A and
2A). However, downwellings are stronger on the nightside than the dayside
because the upper thermal boundary layer imposes a larger temperature contrast
between the surface and the interior (Fig. 1A). This results in strong
downwellings since the contrast determines both thermal buoyancy as well as
material strength through viscosity. High viscosity downwellings dictate the
long-wavelength pattern of flow and thus shepherd the upwellings (also known
as plumes, orange/pink structures in Fig. 1A) along the core–mantle boundary
(CMB) into position between the downwellings. The pinning of plumes by
downwellings is evident by the longitudinal stability of the upwelling and
downwelling tracks over several Gyrs (Fig. 2A).

Both the downwellings and the upwellings exhibit little lateral migration in
longitude, only oscillating back and forth by 45–60 degrees. Low viscosity
upwellings mirror the migration of high viscosity downwellings, maintaining a
constant separation at all times (Fig. 2A). This occurs for the dayside and
the nightside, even though the downwellings are weaker on the dayside. The
evolutionary tracks in Fig. 2A also reveal that upwellings and downwellings
are maintained by a constant draining of the upper and lower thermal boundary
layer, respectively. Hence once the arrangement of downwellings and upwellings
is established in the interior, they tend to persist and no new boundary layer
instabilities occur. Therefore both the day- and nightside only exhibit a
small amount of time-dependent flow.

###  3.2 Hemispheric tectonics: downwellings on dayside

Figure 2: Evolutionary tracks for tectonic regimes of LHS 3844b. Cross-
sections (A3, B3, C3, and D3) show arrangement of upwellings and downwellings
in the mantle. Below, the evolution charts show the longitude of upwellings
(red) and downwellings (blue) as a function of time. Horizontal dashed lines
denote the times of the cross-sections including those in Fig. 1. The
substellar point is at 0∘ longitude and the nightside (90∘ – 270∘) is denoted
by a gray background. Figure 3: Evolutionary tracks for the strong lithosphere
case with a higher nightside temperature (based on observational uncertainty)
and reduced internal heating rate compared to the reference cases (Fig. 2: B,
D). (A) basal heating and $T_{\rm Nightside}=355$ K, (B) basal heating and
$T_{\rm Nightside}=710$ K, (C) mixed heating (half the reference rate) and
$T_{\rm Nightside}=355$ K, and (D) mixed heating (half the reference rate) and
$T_{\rm Nightside}=710$ K. Hemispheric tectonics persist throughout this
parameter variation.

A dominant downwelling on the dayside and upwellings on the nightside are
predicted with a strong lithosphere, independent of the heating mode (Figs.
1B, D and 2B, D). A degree-1 convection pattern is established and hence
hemispheric tectonics operates at the surface. Figures 1B and 2B show the
temperature for the basally heated model with a strong lithosphere. A
prominent downwelling forms on the dayside and descends into the deep mantle
(Fig. 1: B1). It flushes hot material from the lower thermal boundary layer
around the CMB from the dayside to the nightside, thereby promoting plumes on
the nightside that rise to the surface (Fig. 1: B2, B3). This deep mantle flow
from dayside to nightside is accommodated in the post-perovskite layer,
whereas the return flow is established in the perovskite layer above and
delivers cold material from the nightside to the dayside (Fig. 1: B2). The
near-surface advection of cold material from the nightside to the dayside
thickens the upper thermal boundary layer on the dayside, which helps to
sustain the prominent downwelling.

Figure 2B shows that the downwelling (blue track) remains close to the
substellar point (0∘) rather than to the day–night terminator (90∘ or -90∘).
The downwelling displaces hot material laterally along the CMB, thereby
promoting the thickening of the lower thermal boundary layer and hence plume
formation. This leads to the initiation of upwellings on the dayside that are
pushed toward the nightside (Fig. 1: B3). This is evident in the evolutionary
tracks (Fig. 2B, red), where upwellings migrate up to 180∘ from near the
substellar point to rise near the antistellar point (180∘).

A similar degree-1 convection pattern is observed for a planet with a strong
lithosphere and a mixed heating mode (Fig. 1D). The single downwelling is weak
since a majority of the lithosphere is sufficiently strong to resist dynamic
instability (Fig. 1: D3). Furthermore, high mantle temperature ($\approx 3488$
K) due to internal heating causes the downwelling to dissipate quickly and
also reduces the temperature contrast (hence thermal buoyancy) of upwellings.
Although the downwelling is weak, the high mantle temperature ensures material
flows readily due to the temperature dependence of viscosity. Therefore, a
degree-1 flow regime is established with comparable interior velocities as for
a purely basally heated planet where the downwelling is stronger but the
mantle cooler (compare Fig. 1: B2 and Fig. 1: D2). The sweeping of hot
material to the nightside by the downwelling is evident in Fig. 1(D1) and Fig.
2(D).

We also investigated if hemispheric tectonics persist for higher nightside
temperatures based on the uncertainty from the observations (Kreidberg et al.,
2019). For basal heating with nightside temperatures of $355$ and $710$ K, we
find that hemispheric tectonics (downwellings on dayside) persist (Figs. 3A,
B). Similarly, hemispheric tectonics also persist when the internal heating
rate is a factor of $2$ less than the reference models (Figs. 3C, D).

Figure 4: Surface mobility, which is the ratio of the root mean square (rms)
surface velocity to the RMS of the convective velocity in the mantle (Earth
mobility is around 1.3). Red and blue lines show the dayside and nightside,
respectively. The shaded regions on either side of the line represent the
minimum and maximum values.

###  3.3 Hemispheric tectonics: downwellings on the nightside

We observe another hemispheric tectonic regime if the planet has a weak
lithosphere and mixed heating mode (Figs. 1C and 2C). A degree-1 pattern is
established and now strong downwellings occur on the nightside and upwellings
accumulate on the dayside (Fig. 1: C3). Unlike the previous hemispheric regime
(a single downwelling on the nightside), there are several downwellings on the
nightside as demonstrated by two prominent downwelling tracks either side of
the antistellar point (Fig. 2C, blue) and transient downwellings that occur
closer to the day–night terminator. Weak downwellings form on the dayside,
although these tend to remain in the perovskite layer and do not propagate to
depth in the post-perovskite layer (Fig. 1: C3, C2). A prominent upwelling
forms at the antistellar point (180∘) that remains stable and exhibits minimal
lateral movement (Fig. 1: C3), anchored by the two prominent downwellings on
either side (Fig. 2C). The degree-1 pattern of convection is once again driven
by downwellings that propagate to depth, in this case on the nightside. The
downwellings sweep hot material in the post-perovskite layer along the CMB
toward the dayside, where plumes rise and merge at the substellar point. A
return flow is established in the upper perovskite layer where warm material
from the surface of the dayside is advected to the cooler antistellar point
(Fig. 1: C2).

##  4 Implications

###  4.1 Hemispheric tectonics and surface mobility

Our numerical simulations that are constrained by phase curve observations
suggest that LHS 3844b may exhibit persistent hemispheric tectonics due to the
approximately 1000 K temperature contrast between its day- and nightside.
Hemispheric tectonics are not inferred for any solar system planet at present
day. Currently Earth has a dominant degree-2 pattern of convection, as
dictated by downwellings in the Pacific and broad hot anomalies (possibly
upwellings) beneath Africa and the Pacific Ocean. Transient hemispheric
tectonics on Earth may occur during supercontinent formation (landmasses
assemble on one hemisphere), but according to the geological record this is
not a stable and long-lived configuration. Venus does not have a mobile
surface driven by persistent downwellings, but rather its crustal thickness
and mean surface age suggests episodic overturn of the surface (Rolf et al.,
2018). Observations support a prominent role for upwellings to explain surface
features (coronae) and active upwellings may cluster in Venus’ southern
hemisphere (Gülcher et al., 2020). However, these upwellings are not
shepherded into the southern hemisphere by persistent downwellings in the
northern hemisphere.

For Mars, a transient degree-1 convection pattern may have initiated following
a giant impact that generated massive magmatism to source a Tharsis-like
volcanic province (Golabek et al., 2011). However, a mobile surface on Mars
lasted for only half a billion years or so, after which it transitioned to a
stagnant lid (no surface mobility; (Zhang & O’Neill, 2015)). In numerical
models of a generic tidally locked Earth-sized planet, van Summeren et al.
(2011) observe downwellings on the nightside and upwellings on the dayside; we
now find that a similar regime can manifest for a larger rocky planet (LHS
3844b) using phase curve observations and fluid simulations appropriate for a
super-Earth interior. We also discover a previously unrecognized regime of
hemispheric tectonics that is characterized by downwellings on the dayside and
upwellings on the nightside. Hence our data-constrained models emphasize the
importance of hemispheric tectonic regimes to interpret observations of rocky
exoplanets.

For a strong lithosphere, convective stresses in the interior are insufficient
to yield the nightside lithosphere and promote downwellings. Instead, the
nightside lithosphere establishes a surface return flow to the dayside that
accommodates the deep flow of material from the day- to the nightside
initiated by downwellings (Figs. 1B, D). As the cold lithosphere is
transported to the dayside, it becomes warmer and accommodates more
deformation. All models produce mobile surfaces, so even the strong
lithosphere models do not prevent surface motion (Fig. 4). Neither the dayside
nor nightside transition to a stagnant lid with zero surface velocity (Fig.
4). We recover a stagnant lid regime in models (not shown) without plastic
yielding (i.e. $\sigma_{\rm duct}\rightarrow\infty$ MPa) and find no interior
dichotomy due to the lack of downwellings. Instead, weak upwellings are
approximately uniformly distributed but become stifled by interior heat
production that raises the mantle temperature.

We consider lithospheric strength to depend on lithospheric composition, which
we assume to be uniform. Future work could investigate the effects of
temperature dependent yield stress. For our models with a strong lithosphere,
the high temperature on the dayside reduces the viscosity of the lithosphere
sufficiently to allow it to flow and sink. High temperature would likely
decrease the plastic yielding parameter, further facilitating downwellings on
the dayside (Hansen et al., 2019) and preventing downwellings on the nightside
(where the yielding parameter would increase). Hence the asymmetry between the
flow on the dayside and nightside would be further enhanced and increase the
likelihood of hemispheric tectonics. Finally, melting and crustal production
would likely increase the overall mobility of the lithosphere (Lourenço et
al., 2016).

###  4.2 Plumes, volcanism, and observations

Upwellings (plumes) form through two types of interactions with downwellings
(e.g., Tan et al., 2002). First, the leading edge of downwellings sweep hot
material along the CMB to form upwellings (e.g., Fig. 1(B1)). For hemispheric
tectonics, this causes the partitioning of upwellings in one hemisphere and
downwellings in the other hemisphere. Second, upwellings can form beneath
downwellings since cold downwellings trap hot material beneath, thus enabling
incipient plumes to accumulate buoyancy before breaking through and rising to
the surface (e.g., Fig. 1(B3), 30∘). A long-lived cluster of upwellings in a
particular hemisphere may promote extensive magmatism and volcanism. For
example, upwellings can explain the origin of some large igneous provinces
(LIPs) on Earth, which are large surface emplacements of basaltic magma (e.g.,
Coffin & Eldholm, 1994). LIPs drive extensive outgassing (e.g. CO2, Black &
Gibson, 2019) and climate modification. Our models reveal that plumes can be
evenly distributed in the interior, preferentially on the dayside, or
preferentially on the nightside. For a uniform distribution, where upwellings
are pinned in place by downwellings, upwellings only sample deep mantle
material from a small region at the CMB (around 60∘). For hemispheric
tectonics, plumes are instead flushed around the CMB to the opposite
hemisphere, thereby sampling a much larger fraction of the deep interior and
delivering this chemical signature to the near-surface.

We time-average the surface heat flux from the interior and determine its
contribution to the thermal phase curve. The thermal phase curve that is
applied as the surface boundary condition in our reference models is shown in
Fig. 5a. Figure 5b shows the nightside thermal phase curve including the
interior contribution. The dayside phase curve is not shown since the interior
heat flux is negligible compared to stellar insolation such that it has little
influence on the surface temperature. The nightside is a favorable target for
observing the heat contribution from the interior due to temperatures around
$0$ K for an atmosphere-less planet (no reflected light or reradiation of the
host star’s light). An interior heat contribution, as dictated by the interior
convective regime, could produce a variation of the geometric transit depth
(Kipping & Tinetti, 2010). However, the expected spectral signatures are
significantly below the ppm-level in the near-to-mid infrared, which is beyond
the capabilities of current or planned instrumentation (Fig. 5).

On the dayside, which is observed during the secondary transit, the
contribution of the interior heat flux to the eclipse depth would be on the
same order of magnitude as for the transit depth. A further complication for
the secondary transit is that deconvolving the interior flux from stellar
contributions requires knowledge of the stellar spectrum and planetary surface
properties (to determine the reflected, absorbed, and reradiated light). Thus,
it may be difficult to observe the contribution from the interior heat to the
thermal phase curve if heat is conducted through a solid lid. However, if
massive outpourings of melt occur on the surface associated with plumes (akin
to LIPs), this will generate regions of high heat flux and possibly provide a
means to test the existence of hemispheric tectonics with future astrophysical
observations. Furthermore, if upwellings facilitate or source outgassing of
volatiles on one hemisphere, this may perturb the atmospheric composition and
properties on this hemisphere. This may similarly produce a secondary signal
in multiwavelength phase curve observations, in addition to the dominant
signal from the day–night temperature contrast.

Figure 5: (a) Imposed surface temperature motivated by the thermal phase curve
of LHS 3844b. (b) Thermal phase curves of the nightside of LHS 3844b that
include the contribution of the heat flux from the interior.

Parameterized stagnant lid models of volcanic outgassing and atmospheric
erosion show that the lack of an atmosphere for LHS 3844b is consistent with a
volatile poor mantle (Kane et al., 2020). Our stagnant lid models produce a
uniform distribution of upwellings, which implies spatially uniform melting
and outgassing at the surface. Upwellings in the stagnant lid regime are
laterally mobile because they are not anchored in place by downwellings. Kane
et al. (2020) find that the heat-producing element budget, mantle viscosity,
and initial temperature, do not have a significant influence on the size of
the atmosphere of LHS 3844b. However, hemispheric tectonics could lead to
differences in melt production and outgassing between the dayside and
nightside, potentially modifying the size and chemical composition of an
atmosphere. If LHS 3844b is devoid of an atmosphere and has a volatile poor
mantle, we can preclude weakening of the lithosphere by surface water. This
suggests LHS 3844b has a strong lithosphere, lending support to hemispheric
tectonics characterized by a downwelling on the dayside and upwellings on the
nightside.

##  5 Conclusions

Our numerical experiments suggest that hemispheric tectonics may operate on
LHS 3844b, where one hemisphere is characterized by downwellings and the
opposite hemisphere by upwellings. For hemispheric tectonics, upwellings may
lead to preferential melt generation and outgassing on one hemisphere that
could manifest as a secondary signal in phase curve observations. However, the
contribution to the thermal phase curve from the interior flux is on the order
of 15–30 K, which would produce spectral signatures that are significantly
below the ppm level in the near-to-mid infrared, and will therefore be
challenging to detect by current and near-future observations. Nevertheless,
outpourings of melt (extrusive volcanism) with high temperature ($>1000$ K)
fueled by deep mantle upwellings could imprint a signature in the thermal
phase curve. If melting is more prevalent on one hemisphere, the associated
degassing of volatiles could preferentially perturb the composition and
properties of a thin and transient atmosphere on that side of the planet.

T.G.M. and D.J.B. acknowledge SNSF Ambizione Grant 173992. Calculations were
performed on UBELIX (http://www.id.unibe.ch/hpc), the HPC cluster at the
University of Bern. T.L. was supported by the Simons Foundation (SCOL award
#611576 to T.L.) and the SNSF (Early Postdoc.Mobility fellowship
#P2EZP2-178621 to T.L.). B.-O. D. acknowledges support from the SNSF
(PP00P2-190080). This research benefited from discussions and interactions
within the framework of the National Center for Competence in Research (NCCR)
PlanetS supported by the SNSF. Constructive comments from an anonymous
reviewer further sharpened the manuscript.

## References

  * Black & Gibson (2019) Black, B. A., & Gibson, S. A. 2019, Elements, 15, 319, doi: 10.2138/gselements.15.5.319
  * Coffin & Eldholm (1994) Coffin, M. F., & Eldholm, O. 1994, Rev. Geophys., 32, 1, doi: 10.1029/93RG02508
  * Demory et al. (2012) Demory, B.-O., Gillon, M., Seager, S., et al. 2012, Astrophys. J. Lett., 751, L28, doi: 10.1088/2041-8205/751/2/L28
  * Demory et al. (2016) Demory, B.-O., Gillon, M., de Wit, J., et al. 2016, Nature, 532, 207, doi: 10.1038/nature17169
  * Golabek et al. (2011) Golabek, G. J., Keller, T., Gerya, T. V., et al. 2011, Icarus, 215, 346 , doi: 10.1016/j.icarus.2011.06.012
  * Gülcher et al. (2020) Gülcher, A. J. P., Gerya, T. V., Montési, L. G. J., & Munch, J. 2020, Nat. Geosci., 13, 547, doi: 10.1038/s41561-020-0606-1
  * Hansen et al. (2019) Hansen, L. N., Kumamoto, K. M., Thom, C. A., et al. 2019, Journal of Geophysical Research: Solid Earth, 124, 5427, doi: https://doi.org/10.1029/2018JB016736
  * Hernlund & Tackley (2008) Hernlund, J. W., & Tackley, P. J. 2008, Phys. Earth Planet. Inter., 171, 48, doi: 10.1016/j.pepi.2008.07.037
  * Kane et al. (2020) Kane, S. R., Roettenbacher, R. M., Unterborn, C. T., Foley, B. J., & Hill, M. L. 2020, Planet. Sci. J., 1, 36, doi: 10.3847/psj/abaab5
  * Kipping & Tinetti (2010) Kipping, D. M., & Tinetti, G. 2010, Monthly Notices of the Royal Astronomical Society, 407, 2589, doi: 10.1111/j.1365-2966.2010.17094.x
  * Kohlstedt et al. (1995) Kohlstedt, D. L., Evans, B., & Mackwell, S. J. 1995, J. Geophys. Res.-Sol. Ea., 100, 17587, doi: 10.1029/95JB01460
  * Korenaga (2017) Korenaga, J. 2017, J. Geophys. Res.-Sol. Ea., 122, 4064, doi: 10.1002/2016JB013850
  * Kreidberg et al. (2019) Kreidberg, L., Koll, D. D. B., Morley, C., et al. 2019, Nature, 573, 87, doi: 10.1038/s41586-019-1497-4
  * Lourenço et al. (2016) Lourenço, D. L., Rozel, A., & Tackley, P. J. 2016, Earth Planet. Sci. Lett., 439, 18, doi: 10.1016/j.epsl.2016.01.024
  * Moresi & Solomatov (1998) Moresi, L., & Solomatov, V. 1998, Geophys. J. Int., 133, 669, doi: 10.1046/j.1365-246X.1998.00521.x
  * O’Neill & Lenardic (2007) O’Neill, C., & Lenardic, A. 2007, Geophys. Res. Lett., 34, doi: 10.1029/2007GL030598
  * Rolf et al. (2018) Rolf, T., Steinberger, B., Sruthi, U., & Werner, S. 2018, Icarus, 313, 107, doi: 10.1016/j.icarus.2018.05.014
  * Tackley (2000) Tackley, P. 2000, Geochem. Geophy. Geosys., 1, doi: 10.1029/2000GC000036
  * Tackley et al. (2013) Tackley, P., Ammann, M., Brodholt, J., Dobson, D., & Valencia, D. 2013, Icarus, 225, 50, doi: 10.1016/j.icarus.2013.03.013
  * Tackley (2008) Tackley, P. J. 2008, Phys. Earth Planet. Inter., 171, 7, doi: 10.1016/j.pepi.2008.08.005
  * Tan et al. (2002) Tan, E., Gurnis, M., & Han, L. 2002, Geochem. Geophy. Geosys., 3, 1, doi: 10.1029/2001GC000238
  * Unterborn et al. (2015) Unterborn, C. T., Johnson, J. A., & Panero, W. R. 2015, Astrophys. J., 806, 139, doi: 10.1088/0004-637X/806/1/139
  * Valencia et al. (2006) Valencia, D., O’Connell, R. J., & Sasselov, D. 2006, Icarus, 181, 545, doi: 10.1016/j.icarus.2005.11.021
  * van Summeren et al. (2011) van Summeren, J., Conrad, C. P., & Gaidos, E. 2011, Astrophys. J. Lett., 736, L15, doi: 10.1088/2041-8205/736/1/L15
  * Vanderspek et al. (2019) Vanderspek, R., Huang, C. X., Vanderburg, A., et al. 2019, Astrophys. J., 871, L24, doi: 10.3847/2041-8213/aafb7a
  * Zhang & O’Neill (2015) Zhang, S., & O’Neill, C. 2015, Icarus, 265, 187, doi: 10.1016/j.icarus.2015.10.019


# Catalogue of topological electrons and phonons in all allotropes of carbon

Qing-Bo Liu  Hubei Key Laboratory of Optical Information and Pattern
Recognition, School of Optical Information and Energy Engineering, School of
Mathematics and Physics, Wuhan Institute of Technology, Wuhan, 430073, China
Xiang-Feng Yang  Zhe-Qi Wang  School of Physics and Wuhan National High
Magnetic Field Center, Huazhong University of Science and Technology, Wuhan
430074, People’s Republic of China.  Institute for Quantum Science and
Engineering, Huazhong University of Science and Technology, Wuhan, Hubei
430074, China.  Ziyang Yu<EMAIL_ADDRESS>Lun Xiong<EMAIL_ADDRESS>Hubei Key Laboratory of Optical Information and Pattern Recognition, School of
Optical Information and Energy Engineering, School of Mathematics and Physics,
Wuhan Institute of Technology, Wuhan, 430073, China  Hua-Hua Fu
<EMAIL_ADDRESS>School of Physics and Wuhan National High Magnetic Field
Center, Huazhong University of Science and Technology, Wuhan 430074, People’s
Republic of China.  Institute for Quantum Science and Engineering, Huazhong
University of Science and Technology, Wuhan, Hubei 430074, China.

###### Abstract

Carbon, as one of the most common element in the earth, constructs hundreds of
allotropic phases to present rich physical nature. In this work, by combining
the ab inito calculations and symmetry analyses method, we systematically
study a large number of allotropes of carbon (703), and discovered 315 ideal
topological phononic materials and 32 topological electronic materials. The
ideal topological phononic nature includes single, charge-two, three, four
Weyl phonons, the Dirac or Weyl nodal lines phonons, and nodal surfaces
phonons. And the topological electron nature includes topological insulator,
(Type-II) Dirac points, triple nodal points, the Dirac (Weyl) nodal lines,
quadratic nodal lines and so on. For convenience, we take the $uni$ in SG 178
and $pbg$ in SG 230 as the examples to describe the topological features in
the main. We find that it is the coexistence of single pair Weyl phonons and
one-nodal surfaces phonons in the $uni$ in SG 178, which can form the single
surface arc in the (100) surface BZ and isolated double-helix surface states
(IDHSSs)in the (110) surface BZ. In topological semimetal $pbg$ in SG 230, we
find that the perfect triple degenerate nodal point can be found in the near
Fermi level, and it can form the clear surface states in the (001) and (110)
surface BZ. Our work not only greatly expands the topological features in all
allotropes of carbon, but also provide many ideal platforms to study the
topological electrons and phonons.

1\. Introduction.

Topological materials, including topological insulators (TIs), topological
semimetals (TSMs) [1, 2, 3, 4, 5, 6, 7, 8, 9], topological superconductors
(TSCs) [9, 10, 11, 12] and so on, have been received intense studies due to
their protected boundary states and prospects for the future applications in
quantum devices in the past a dozen years decades. And nodal point phonons,
nodal lines phonons and nodal surface phonons have all been attracting the
tremendous research interests, because the topological phonons can realize
particular phonon-based devices [13, 14], such as phonon diode effect in
hexagonal honeycomb lattice. Recently years, single Weyl phonons with Chern
number $\pm$1 [15, 16, 17, 18, 19, 20], charge-2 [21, 22], 3 [23, 24] and 4
[3, 25] Weyl phonons with charges of $\pm$2, $\pm$3 and $\pm$4 have been
observed in the realist materials and they can form the single, double, triple
and quadruple-helicoid surface states. At the same time, the spin-1 Weyl
phonons, charge-2 Dirac phonons and helicoid nodal lines phonons are all
observed in the experiment [26, 27]. So, to explore novel topological
electrons and phonons has been another central topic in topological physics
and materials.

Besides the well-known the allotropes of carbon of graphite [28], diamond
[29], carbon nanotubes [30], graphene [31] and fullerenes [32], more than 700
carbon allotropes have been theoretically predicted or experimentally
synthesized. In recent years, the magic angle twisted bilayer graphene (TBG)
has attracted the an intense interest of researchers, because it has many
interesting physical properties, such as superconductivity [33, 34, 35] and
topologically nontrivial electronic states. Beyond TBG, the topological
features of carbon allotropes have been widely studied in the electronic and
phononic systems [36, 37, 38], For example, topological semimetals phases have
been predicted in nanostructure carbon allotropes, body centered orthorhombic
C16 [39] monoclinic C16 [40] and C40 [41], and others [42, 43], including the
topological features of nodal points [42], nodal nets [41, 42], nodal rings
[36] and nodal surfaces [43]. And topological phonons are also predicted in
the allotropes of carbon, such as nodal rings phonons [44], straight nodal
lines phonons [45], single Weyl phonons [45], charge-2 Weyl phonons [45].
However, no reports have reported the all topological electrons and phonons of
all allotropes of carbon so far.

Figure 1: The schematic procedure for discovering the topological electrons
and phonons of all carbon allotropes.

In this work, by combining the ab initio calculation and symmetry analyses
method, we firstly systematically study a large number of carbon allotropes
(703), and discovered 315 ideal topological phononic materials and 32
topological electronic materials. The ideal topological phononic features
include charge-one, two, three, four Weyl phonons, Dirac or Weyl phonons,
Dirac or Weyl nodal lines phonons, and nodal surfaces phonons. And the
topological electron nature includes topological insulators, Dirac or Weyl
points, triple nodal points, the Dirac or Weyl nodal lines, quadratic nodal
lines in all allotropes of carbon. We mainly discuss the topological phononic
(electronic) nature of $uni$ ($pbg$) in SG 178 (230) in the main. We find that
it is the coexistence of single pair Weyl phonons with charge -1 and one-nodal
surfaces phonons in the $uni$ in SG 178. Moreover, it can form the single
surface arc in the (100) surface BZ and isolated double-helix surface states
in the (110) surface BZ. In addition, another example of topological semimetal
$pbg$ in SG 230 can form the ideal triple nodal point near the Fermi level.
And it can also form clear surface states in (001) and (110) surfaces states.
We discuss the topological electronic and phononic features of other carbon
allotropes in the Supplementary materials [46]. Our theoretical results not
only find many ideal platform to study the topological electron and phonons in
allotropes of carbon, but also greatly expands the topological features in all
carbon allotropes.

Figure 2: (a) Crystal structure of $uni$ in SG 178 in a primitive cell. (b)
The bulk BZ of $uni$ and the (100) (red square) and (110) (blue square)
surface BZ. (c) The distribution of single pair monopole WPs (green dots with
Chern number -1) and one-nodal surfaces (yellow region) between the $17^{th}$
and $18^{th}$ phonon bands of $uni$ in the first BZ. (d) The phonon spectra of
$uni$ along the high-symmetry paths and the corresponding phononic density of
states (DOS). A red box is marked along the high-symmetry paths as shown in
Fig. 1(e). (f) The evolution of the average position of Wannier centers for
single pair monopole WPs.

2\. Materials and method

The crystallographic data of all allotropes of carbon can be taken from the
international database SACADA [47]. The phononic dispersions and elecronic
bands of allotropes of carbon are calculated by the density functional theory
(DFT) using the Vienna _ab initio_ Simulation Package (VASP) with the
generalized gradient approximation (GGA) in the form of Perdew-Burke-Ernzerhof
(PBE) function for the exchange-correlation potential [48, 49, 50]. An
accurate optimization of structural parameters is employed by minimizing the
interionic forces less than 0.001 $\textrm{eV}/\textrm{{\AA}}$ and an cutoff
energy at 450 eV. The BZ is gridded with 3$\times$3$\times$3 _k_ points. Then
the phononic spectra are gained using the density functional perturbation
theory (DFPT), implemented in the Phonopy Package [51]. The force constants
are calculated using a $2\times 2\times 2$ supercell. To reveal the
topological nature of phonons, we construct the phononic Hamiltonian of tight-
binding (TB) model and the surface local density of states (DOS) with the
open-source software Wanniertools code [52] and the surface Green’s functions
[53]. The Chern numbers or topological charge of WPs are calculated by Wilson
loop method [54]. The electronic surface states have been performed using the
open-source code WANNIERTOOLS based on the Wannier tight-binding model
constructed using the WANNIER90 code. The irreps of the electronic and
phononic bands are computed by the program $ir2tb$ on the electronic and
phononic Hamiltonian of TB model [55].

3\. Results and discussion

_3.1 Research the topological electrons and phonons_

The guiding principle of our search the topological electrons and phonons in
all carbon allotropes as shown in Fig. 1. Firstly, among 703 carbon allotropes
structures in the international database SACADA [47], we will exclude the no
topological features of space groups, according to the irreducible
representations of high-symmetry lines (HLs) and high-symmetry points (HPs) in
space groups [56, 57]. Then, we can obtain the 505 allotropes of carbon which
can exist the topologically electronic or phononic properties. Secondly, their
dynamical stabilities were examined by using the density functional
perturbation theory (DFPT) calculations, which yields 329 dynamical stable
allotropes of carbon. Finally, by using the first-principle calculations, we
can divide 329 dynamical stable allotropes of carbon into 32 perfect
topological electronic materials (TEMs) including topological insulators
(TIs), Dirac semimatals (DSMs), Weyl semimatals (WSMs) and so on, and 321
perfect topological phononic materials (TPMs) including Weyl phonons (WPs),
Dirac phonons (DPs), triple-fold phonons (TPs) and so on as shown in
Supplementary materials [46].

_3.2 Symmetry analysis_

Table 1: The complete list of allotropes of carbon with topological properties in SGs 178 and 230. The first column indicates the space groups (SGs), the second column indicates the number of allotropes of carbon in SACADA, the third column indicates the SG symbol, and the fourth column shows the types of topological features (TFs). The red letters stand for the topological insulators or semimatel.  SGs |  Numbers |  Names |  TFs
---|---|---|---
P6122(No. 178) |  29 56 |  uni unj |  Weyl points and nodal surfaces
Ia$\bar{\text{3}}$d(No. 230) |  13 53 |  lcs pbg |  triple nodal points

To elucidate the topologically nontrivial features of single pair Weyl phonons
and one-nodal surface phonons in SG 178 in Table-I, we can prove that it exist
the single pair Weyl phonons at high-symmetry point K and one-nodal surface
phonons at $k_{z}$= $\pm\pi$ planes by using a two-band $k\cdot{p}$ model as

$\mathcal{H}_{kp}(k)=g_{x}(k)\sigma_{x}+g_{y}(k)\sigma_{y}+g_{z}(k)\sigma_{z},$
(1)

where $k=(k_{x},k_{y},k_{z})$, $\sigma_{x,y,z}$ represents the three Pauli
matrices, and $g_{x,y,z}(k)$ represents the complex functions versus $k_{x}$,
$k_{y}$ and $k_{z}$.

Let us first consider a three screw symmetry $\\{C_{31}^{+}|00\frac{1}{3}\\}$
and a two-fold screw symmetry $\\{C_{21}^{{}^{\prime\prime}}|00\frac{1}{2}\\}$
at K point in SG 178. Under a 2D irreps:R3 [56, 57], the representation
matrixes of $\\{C_{31}^{+}|00\frac{1}{3}\\}$ and
$\\{C_{21}^{{}^{\prime\prime}}|00\frac{1}{2}\\}$ can be described as

$\displaystyle\\{C_{31}^{+}|00\frac{1}{3}\\}=\left[\begin{array}[]{ll}-\frac{1}{2}&\frac{\sqrt{3}}{2}\\\
-\frac{\sqrt{3}}{2}&-\frac{1}{2}\end{array}\right],$
$\displaystyle\\{C_{21}^{{}^{\prime\prime}}|00\frac{1}{2}\\}=\left[\begin{array}[]{ll}0&1\\\
1&0\end{array}\right],$

So,the $k\cdot{p}$-invariant Hamiltonian at point K is derived as

$\displaystyle\mathcal{H}_{kp}=ak_{x}\sigma_{x}+bk_{z}\sigma_{y}+ck_{y}\sigma_{z},$

where a, b and c are constant coefficients. According to the two-band
$k\cdot{p}$ Hamiltonian, we find that it has a $k$ dispersion along all
phononic dispersion’s direction at K point.

Then, we will prove that it can exist the one nodal surfaces at at $k_{z}$=
$\pm\pi$ planes in the nonsymmorphnic SG 178. We continue to consider the skew
axial symmetry $S_{2z}$: $(x,y,z,t)\mapsto(-x,-y,-z+\frac{1}{2},t)$, which
indicates that $S_{2z}$ inverses $k_{x}$ and $k_{y}$ while preserves $k_{z}$
in the _k_ -space. One can drive that $(S_{2z})^{2}$ = $T_{001}$ =
$e^{-ik_{z}}$, where $T_{001}$ is the translation along the $x$ direction by a
lattice constant. As we know that $\mathcal{T}$ is antiunitary and inverses
$k$ with the relation $T^{2}$ = 1. So the compound symmetry $S$ = $T$$S_{2z}$
is antiunitary and only inverses $k_{z}$. Since [$T$, $S_{2z}$] =0, $S$
satisfies

$S^{2}=e^{-ik_{z}}.$ (4)

The crossing on the surface can be comprehended as a result of the Kramer’s
degeneracy. Under the operation $S$, we find that any point on the $k_{z}=\pi$
is invariant. Moreover, from Eq. (2), one can drive that the antiunitary
symmetry satisfies $S^{2}$ = -1 on the whole $k_{z}=\pi$ plane. Thus, the two
Kramer’s degeneracy can arise at the every point in this plane. Away from this
plane, this Kramer’s degeneracy is generically destroyed owing to the loss of
symmetry protections. In a word, a nodal surface should be formed at the
$k_{z}=\pi$ plane.

Figure 3: The surface phonon dispersions and isofrequency surface contours of
$uni$ in SG 178. (a) The surface phonon dispersion on the (100) surface along
the high-symmetry
$\overline{\textrm{Y}}$-$\overline{\textrm{X}}$-$\overline{\textrm{$\Gamma$}}$-$\overline{\textrm{Y}}$.
(b) The isofrequency surfaces contours of (100) surface BZ at 36.155 THz. (c)
The surface phonon dispersion on the (110) surface along the high-symmetry
$\widetilde{\textrm{X}}$-$\widetilde{\textrm{$\Gamma$}}$-$\widetilde{\textrm{Y}}$.
(d) The isofrequency surfaces contours 0f (110) surface BZ at 36.155 THz.

Next, we will derive the $k\cdot{p}$ model of $\Gamma$ point in SG 230 in
Table-I. The point $\Gamma$ belongs to point group $O_{h}$, which includes
five symmetries $\\{S_{61}^{-}|000\\}$,
$\\{\sigma_{x}|\frac{1}{2}\frac{1}{2}0\\}$,
$\\{\sigma_{z}|\frac{1}{2}0\frac{1}{2}\\}$, $\\{C_{2c}|00\frac{1}{2}\\}$ and
time-reversal symmetry $\mathcal{T}$. Under a 3D irreps:$\Gamma_{4}^{+}$, the
representation matrixes of there symmetries can be shown

$\displaystyle\\{S_{61}^{-}|000\\}=\left[\begin{array}[]{lll}0&0&1\\\ 1&0&0\\\
0&1&0\\\ \end{array}\right],$
$\displaystyle\\{\sigma_{x}|\frac{1}{2}\frac{1}{2}0\\}=\left[\begin{array}[]{lll}-1&0&0\\\
0&1&0\\\ 0&0&-1\end{array}\right],$
$\displaystyle\\{\sigma_{z}|\frac{1}{2}0\frac{1}{2}\\}=\left[\begin{array}[]{lll}0&1&0\\\
1&0&0\\\ 0&0&-1\end{array}\right],$
$\displaystyle\\{C_{2c}|00\frac{1}{2}\\}=\left[\begin{array}[]{lll}-1&0&0\\\
0&1&0\\\ 0&0&-1\end{array}\right],$ $\displaystyle\mathcal{T}=\mathcal{K}.$

The $\mathcal{K}$ is a complex conjugate operator. So, the three-band
$k\cdot{p}$-invariant Hamiltonian at point $\Gamma$ is derived as

$\displaystyle\mathcal{H}_{kp}=\left[\begin{array}[]{lll}0&mk_{x}k_{y}&nk_{x}k_{z}\\\
mk_{x}k_{y}&0&qk_{y}k_{z}\\\ nk_{x}k_{z}&qk_{y}k_{z}&0\end{array}\right],$

where m, n and q are real constant coefficients. According to this equation,
we find that it can form quadratic triple nodal point (QDNP) at the $\Gamma$
point in SG 230.

_3.3 Topological phononic features of $uni$ in SG 178 _

The crystal structure of $uni$ in SG 178 is shown in Fig. 2(a), which includes
6 carbon atoms in a unit cell. The bulk BZ, (110) surface BZ (blue square) and
(100) surface BZ (red square) are shown in Fig. 2(b). The phononic bands and
the corresponding phononic density of states (DOS) of $uni$ are drawn in Fig.
2(d), where no imaginary frequencies indicate that this material is
thermodynamically stable. From the first glance, we can find that two single
pair Weyl points are located at the K point formed by 17th and 18th bands, and
their distribution of single pair Weyl phonons (green dots) are shown in Figs.
2(c) and 2(e) with Chern number -1 in Fig. 2(f). And two-fold degeneracy bands
can be observed along the high-symmetry lines A-H-L in Fig. 2(d), which can
form one nodal surfaces at $k_{z}=\pm\pi$ planes (yellow region) in Fig. 2(c).

In order to study the topological features of $uni$, we calculate the (100)
and (110) surface states of $uni$ in Fig 3. As shown in Fig. 1(b), the high-
symmetry points $\Gamma$ and M (L and A) are projected to
$\overline{\textrm{$\Gamma$}}$ ($\overline{\textrm{Y}}$) in the (001) surface
in Fig. 2(b). In (100) surface, one can see that the surface arcs surface
states can be formed by the single Weyl phonons and one nodal surfaces in Fig.
3(a), which is for the first time reported in the phononic systems. We also
calculate the sofrequency surface at frequency 36.155 THz and we observe that
the isolated surface acrs can be formed by a Weyl phonon with charge of -1 in
Fig. 3(b). In (110) surface BZ, the high-symmetry points $\Gamma$ and K (A)
are projected to $\widetilde{\textrm{$\Gamma$}}$ ($\widetilde{\textrm{Y}}$) in
Fig. 2(b). We find that two Weyl phonons with charge of -1 are projected o the
same point which can form a charge-two Weyl phonons at
$\widetilde{\textrm{$\Gamma$}}$, leading to isolated double-helix surface
states in Figs. 3(c). And we also can observe the isolated double-helix
surface arc states in the sofrequency surface at frequency 36.155 THz as shown
in Fig. 3(d). So, these novel physical features confirm that it exit the
topological nontrivial nature in $uni$.

_3.4 Topological electronic features of $pbg$ in SG 230 _

Next, we will disscus the topological electronic features of $pbg$ in SG 230
in Figs. 4 and 5. The unit cell of $pbg$ is shown in Fig. 4(a), which includes
48 carbon atoms. The 3D BZ, (001) (blue square) and (110) (red square) surface
BZ are shown in Fig. 4(b). The electronic bands and density of states (DOS) of
$pbg$ are draw in Fig. 4(c), where blue lines (dashed line) stand for the
result of DFT (Wannier90 fitting). We find that a ideal triple nodal point can
exist near the Fermi level ($E_{f}$) at the point $\Gamma$.

Like to $uni$ in SG 178, we also calculate the (001) and (110) surface states
of $pbg$ in Fig. 5. The $\Gamma$ and H points (P and H) are projected to
$\overline{\textrm{$\Gamma$}}$ ($\overline{\textrm{M}}$) in the (001) surface
in Fig. 4(b). At the first glance, we find that at the point
$\overline{\textrm{$\Gamma$}}$, there are three surface arcs surfaces states
projected by triple nodal points in Fig. 5(a). And we can also observe the
triple nodal point at the Fermi level as shown in Fig. 5(b). In (110) surface,
we can also clearly observe the three surface arcs surfaces states projected
by triple nodal points at point $\widetilde{\textrm{$\Gamma$}}$ in Fig. 5(c)
and triple nodal point at $\widetilde{\textrm{$\Gamma$}}$ in Fig. 5(d). These
nontrivial topological feature all prove that $pbg$ is topological.

Figure 4: (a) Crystal structure of $pbg$ in SG 230 in a primitive cell. (b)
The bulk BZ of $pbg$ and the (001) (blue square) and (110) (red square)
surface BZ. (c) The DFT (blue solid lines) and Wannier90 (red dashed lines)
electronic bands of $pbg$ along the high-symmetry paths and the corresponding
electronic density of states (DOS). Figure 5: The surface electronic
dispersions and isofrequency surface contours of $pbg$ in SG 230. (a) The
surface phonon dispersion on the (100) surface along the high-symmetry
$\overline{\textrm{M}}$-$\overline{\textrm{$\Gamma$}}$-$\overline{\textrm{X}}$.
(b) The isofrequency surfaces contours of (001) surface BZ at Fermi level. (c)
The surface electronic dispersion on the (110) surface along the high-symmetry
$\widetilde{\textrm{K}}$-$\widetilde{\textrm{$\Gamma$}}$-$\widetilde{\textrm{X}}$.
(d) The isofrequency surfaces contours 0f (110) surface BZ at Fermi level.

4\. Conclusion

In summary, performing symmetry analysis and first principles calculations, we
systematically study 703 allotropes of carbon, and discovered 315 ideal
topological phononic materials, including single, charge-two, three, four Weyl
phonons, the Dirac or Weyl nodal lines phonons, and nodal surfaces phonons,
and 32 topological electronic materials, including topological insulator,
(Type-II) Dirac points, triple nodal points, the Dirac (Weyl) nodal lines,
quadratic nodal lines and so on. In $uni$ carbon in SG 178, we find that it is
the coexistence of single pair Weyl phonons and one-nodal surfaces phonons in
the uni in SG 178, which can form the single surface arc in the (100) surface
BZ and isolated double-helix surface states in the (110) surface BZ. Another
example, we find we find that the perfect triple nodal point can be found in
the near Fermi level, and it can form the clear triple surface arc states in
the (001) and (110) surface BZ. And more carbon allotropes are tabulated in
the Supplementary materials [46]. Our work not only greatly expands the
topological features in all allotropes of carbon, but also provide many ideal
platforms to study the topological electrons and phonons.

Conflict of interest

The authors declare that they have no conflict of interest.

Acknowledgements.

This work is supported by the National Science Foundation of China with Grants
Nos. 11774107, 12147113, 12104348 and U20A2077, by the Science and Technology
Department of Hubei Provincial with Grant No. 2022CFD041, and partially by the
National Key R&D Program of China (2021YFC2202300).

Author contributions

Ziyang Yu, Lun Xiong and Hua-Hua Fu proposed and supervised the project. Qing-
Bo Liu carried out the high-throughput calculations. All authors contributed
to writing of the manuscript.

Author contributions

Supplementary materials to this article can be found online at XXXX.

## References

  * [1] B. Bradlyn, J. Cano, Z. Wang, M. G. Vergniory, C. Felser, R. J. Cava, and B. A. Bernevig, Begyond Dirac and Weyl ferimions: unconvertianl quasiparticles in conventional crystals, Science 353, 6299 (2016).
  * [2] X. Huang, L. Zhang, Y. Long, P. Wang, D. Chen, Z. Yang, H. Liang, M. Xue, H. Weng, Z. Fang, X. Dai, and G. Chen, Obversation of the chiral-anomaly-induced negative magnetoresistance in 3D Weyl semimetal TaAs, Phys. Rev. X 5, 031023 (2015).
  * [3] T. T. Zhang, R. Takahashi, C. Fang, and S. Murakami, Twofold quadruple Weyl nodes in chiral cubic crystals, Phys. Rev. B 102, 125148 (2020).
  * [4] S. Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. L. Zhang, R. Sankar, G. Q. Chang, Z. J. Yuan, C. C. Lee, S. M. Huang, H. Zheng, J. Ma, D. S. Sanchez, B. K. Wang, A. Bansil, F. C. Chou, P. P. Shibayev, H. Lin, S. Jia, and M. Z. Hasan, Discovery of a Weyl fermion semimental and topological Fermi arcs, Science 349, 613 (2015).
  * [5] 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, Observation of Weyl nodes in TaAs, Nat. Phys. 11, 724 (2015).
  * [6] N. P. Armitge, E. J. Mele, and A. Vishwanath, Weyl and Dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * [7] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Classification of topological insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008).
  * [8] C. L. Kane and E. J. Mele, Z2 Topological order and the quantum spin Hall effect, Phys. Rev. Lett. 95, 146802 (2005).
  * [9] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * [10] R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Majorana Fermions and a topological phase transition in semiconductor-superconductor heterostructures, Phys. Rev. Lett. 105, 077001 2010).
  * [11] S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani, Observation of Majorana fermions in ferromagnetic atomic chains on a superconductor, Science 346, 602.
  * [12] X. L. Qi and S. C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * [13] L. F. Zhang, J. Ren, J. S. Wang, and B. W. Li, Topological nature of the phonon Hall effect, Phys. Rev. Lett. 105, 225901 (2010).
  * [14] Y. Liu, X. Chen, and Y. Xu, Topological phononics: from fundamental models to real materials. Adv. Funct. Mater., 30, 1904784 2020.
  * [15] L. Lu, L. Fu, J. D. Joannopoulos, and M. Soljacic, Weyl points and line nodes in gyroid photonic crystals, Nat. Photon. 7, 294 (2013).
  * [16] L. Lu, Z. Y. Wang, D. X. Ye, L. X. Ran, L. Fu, J. D. Joannopoulos, and M. Soljaci, Experimental observation of Weyl points, Science 349, 622 (2015).
  * [17] Q. Xie, J. Li, S. Ullah, R. Li, L. Wang, D. Li, Y. Li, S. Yunoki, and X.-Q. Chen, Phononic Weyl points and one-way topologically protected nontrivial phononic surface arc states in noncentrosymmetric WC-type materials, Phys. Rev. B. 99, 174306 (2019).
  * [18] B. W. Xia, R. Wang, Z. J. Chen, Y. J. Zhao, and H. Xu, Symmetry-protected ideal type-II Weyl phonons in CdTe. Phys. Rev. Lett. 123, 065501 (2019).
  * [19] G. Ding, F. Zhou, Z. Zhang, Z. M. Yu, and X. Wang, Charge-two Weyl phonons with type-III dispersion. Phys. Rev. B. 105, 134303 (2022).
  * [20] B. Zheng, B. Xia, R. Wang, Z. Chen, J. Zhao, Y. Zhao, and H. Xu, Ideal type-III nodal-ring phonons, Phys. Rev. B 101, 100303(R) (2020).
  * [21] Q.-B. Liu, Y. Qian, H.-H. Fu, and Z. Wang, Symmetry-enforced Weyl phonons, npj Comput. Mater. 6, 95 (2020).
  * [22] R. Wang, B.-W. Xia, J. Chen, B.-B. Zheng, J. Zhao, and H. Xu, Symmetry-Protected Topological Triangular Weyl Complex. Phys. Rev. Lett. 124, 105303 (2020).
  * [23] X. Wang, F. Zhou, Z. Zhang, Z.-M. Yu, and Y. Yao, Hourglass charge-three Weyl phonons, Phys. Rev. B 106, 214309 (2022).
  * [24] G. Liu, Z. Chen, P. Wu, and H. Xu, Triple hourglass Weyl phonons, Phys. Rev. B 106, 214308 (2022).
  * [25] Q.-B. Liu, Z. Wang, and H.-H. Fu, Charge-four Weyl phonons, Phys. Rev. B 103, L161303 (2021).
  * [26] T.-T. Zhang, Z. D. Song, A. Alexandradinata, H.-M. Weng, C. Fang, L. Lu, and Z. Fang, Double-weyl phonons in transition-metal monosilicides, Phys. Rev. Lett. 120, 016401 (2018).
  * [27] H. Miao, T.-T. Zhang, L. Wang, D. Meyers, A.-H. Said, Y.-L. Wang, Y.-G. Shi, H.-M. Weng, Z. Fang, and M.-P.-M. Dean, Observation of double Weyl phonons in parity-breaking FeSi, Phys. Rev. Lett. 121, 035302 (2018).
  * [28] W. Wang, S. Dai, X. Li, J. Yang, D. J. Srolovitz, and Q. Zheng, Measurement of the cleavage energy of graphite, Nat. Commun. 6, 7853 (2015).
  * [29] E. A. Ekimov, V. A. Sidorov, E. D. Bauer, N. N. Mel’nik, N. J. Curro, J. D. Thompson, and S. M. Stishov, Superconductivity in diamond, Nature 428, 542-545(2004).
  * [30] S. Iijima, Helical microtubules of graphitic carbon, Nature (London) 354, 56 (1991).
  * [31] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * [32] R. C. Haddon, Chemistry of the Fullerenes: The Manifestation of Strain in a Class of Continuous Aromatic Molecules, Science 261, 5128 (1993).
  * [33] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43-50(2018).
  * [34] Y. Cao, D. Rodan-Legrain, J. M. Park, N. FQ Yuan, K. Watanabe, T. Taniguchi, R. M Fernandes, L. Fu, and P. Jarillo-Herrero, Nematicity and competing orders in superconducting magic-angle graphene, Science 372, 264-271(2021).
  * [35] Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Pauli-limit violation and re-entrant superconductivity in moiré graphene, Nature 595, 526-531(2021).
  * [36] H. Weng, Y. Liang, Q. Xu, R. Yu, Z. Fang, X. Dai, and Y. Kawazoe, Topological node-line semimetal in three-dimensional graphene networks, Phys. Rev. B 92, 045108 (2015).
  * [37] Y. Chen, Y. Xie, S. A. Yang, H. Pan, F. Zhang, M. L. Cohen, and S. Zhang, Nanostructured Carbon Allotropes with Weyl-like Loops and Points, Nano Lett. 15, 6974 (2015).
  * [38] J.-T. Wang, C. Chen, and Y. Kawazoe,Topological nodal line semimetal in an orthorhombic graphene network structure, Phys. Rev. B 97, 245147 (2018).
  * [39] J.-T. Wang, H. Weng, S. Nie, Z. Fang, Y. Kawazoe, and C. Chen, Body-Centered Orthorhombic C16: A Novel Topological Node-Line Semimetal, Phys. Rev. Lett. 116, 195501 (2016).
  * [40] Y. Cheng, X. Feng, X. Cao, B. Wen, Q. Wang, Y. Kawazoe, and P. Jena, Body-Centered Tetragonal C16: A Novel Topological Node-Line Semimetallic Carbon Composed of Tetrarings, Small 13, 1602894 (2017).
  * [41] J.-T. Wang, S. Nie, H. Weng, Y. Kawazoe, and C. Chen, Topological Nodal-Net Semimetal in a Graphene Network Structure, Phys. Rev. Lett. 120, 026402 (2018).
  * [42] X. Dong, M. Hu, J. He, Y. Tian, and H.-T. Wang, A new phase from compression of carbon nanotubes with anisotropic Dirac fermions, Sci. Rep. 5, 10713 (2015).
  * [43] C. Zhong, Y. Chen, Y. Xie, S. A. Yang, M. L. Cohen,and S. Zhang, Towards three-dimensional Weyl-surface semimetals in graphene networks, Nanoscale 8, 7232 (2016).
  * [44] Y.-J. Jin, Z.-J. Chen, B.-W. Xia, Y.-J. Zhao, R. Wang, and H. Xu, Ideal intersecting nodal-ring phonons in bcc C8, Phys. Rev. B. 98, 220103(R) (2018).
  * [45] Q.-B. Liu, Z.-Q. Wang, and H.-H. Fu, Topological phonons in allotropes of carbon, Mater. Today Phys. 22, 100694 (2022).
  * [46] Supplementary materials in XXXX
  * [47] R. Hoffmann, A. A. Kabanov, A. A. Golov, D. M. Proserpio, Homo Citans and Carbon Allotropes: For an Ethics of Citation, Angew. Chem. Int. Ed. 55, 10962–10976 (2016).
  * [48] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. B 136, A1133 (1964).
  * [49] P. E. Blöchl, Projector augmentad-wave method, Phys. Rev. B 50, 17953 (1994).
  * [50] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * [51] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * [52] Q.-S. Wu, S.-N. Zhang, H.-F. Song, M. Troyer, and A. A. Soluyanov, WannierTools: an open-source software package for novel topological materials, Comput. Phys. Commun. 243, 110 (2019).
  * [53] M. P. L. Sancho, J. M. L. Sancho, J. M. L. Sancho, and J. Rubio, Highly convergent schemes for the calculations of bulk and surface Green functions, J. Phys. F: Metal Phys., 15, 851 (1985).
  * [54] R. Yu, X.-L. Qi, A. Bernevig, Z. Fang, and X. Dai, Equivalent expression of $Z_{2}$ topological invariant for band insulators using the non-abelian Berry connection. Phys. Rev. B 84, 075119 (2011).
  * [55] J.-C. Gao, Q.-S. Wu, C. Persson, and Z. Wang, Irvsp: To obtain irreducible representations of electronic states in the VASP, Comput. Phys. Commun., 261, 107760 (2021).
  * [56] Z.-M. Yu, Z. Zhang, G.-B. Liu, W. Wu, X.-P. Li, R.-W. Zhang, S.-A. Yang, and Y. Yao, Encyclopedia of emergent particles in three-dimensional crystals, Sci. Bull. 10, 023 (2022).
  * [57] C. J. Bradley and A. P. Cracknell, The mathematical theory of symmetry in solids: representation theory for point groups and space groups, Oxford University Press (2009).