# Crystalline axion electrodynamics in charge-ordered Dirac semimetals

Julian May-Mann  Department of Physics, Stanford University, Stanford, CA
94305, USA  Mark R. Hirsbrunner  Department of Physics, University of Illinois
Urbana-Champaign, Urbana, IL 61801, USA  Lei Gioia  Walter Burke Institute for
Theoretical Physics, Caltech, Pasadena, CA, USA  Department of Physics,
Caltech, Pasadena, CA, USA  Perimeter Institute for Theoretical Physics,
Waterloo, Ontario N2L 2Y5, Canada  Taylor L. Hughes  Department of Physics,
University of Illinois Urbana-Champaign, Urbana, IL 61801, USA

###### Abstract

Three-dimensional Dirac semimetals can be driven into an insulating state by
coupling to a charge density wave (CDW) order. Here, we consider the quantized
crystalline responses of such charge-ordered Dirac semimetals, which we dub
Dirac-CDW insulators, in which charge is bound to disclination defects of the
lattice. Using analytic and numeric methods we show the following. First, when
the CDW is lattice-commensurate, disclination-line defects of the lattice have
a quantized charge per length. Second, when the CDW is inversion-symmetric,
disclinations of the lattice have a quantized electric polarization. Third,
when the CDW is lattice-commensurate and inversion-symmetric, disclinations
are characterized by a “disclination filling anomaly” — a quantized difference
in the total charge bound to disclination-lines of Dirac-CDW with open and
periodic boundaries. We construct an effective response theory that captures
the topological responses of the Dirac-CDW insulators in terms of a total
derivative term, denoted the $R\wedge F$ term. The $R\wedge F$ term describes
the crystalline analog of the axion electrodynamics that are found in Weyl
semimetal-CDW insulators. We also use the crystalline responses and
corresponding response theories to classify the strongly correlated
topological phases of three-dimensions Dirac-semimetals.

##  I Introduction

Three-dimensional (3D) Dirac semimetals (DSMs) and Weyl semimetals (WSMs) are
quintessential examples of gapless phases of matter that harbor topological
features [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. These systems are
characterized by linear bulk band crossings (nodal points) that are well-
described at low energies by the 3D Dirac and Weyl equations, respectively.
The gaplessness of Weyl nodes can be viewed as a consequence of band topology.
Each node is a monopole of Berry curvature in momentum space, and, since the
total number of Berry curvature monopoles must vanish, single Weyl nodes
cannot be added or removed. This reasoning also indicates that Weyl nodes come
in pairs having opposite monopole charge [13, 2, 3, 14]. The gaplessness of
Dirac nodes is also protected, but requires additional symmetry, e.g., time-
reversal, inversion, and spatial rotations [15, 16]. While the individual
nodes of the DSM and WSM are locally stable in momentum space, pairs of nodes
can acquire a gap. Normally, this cannot occur since the nodes are located at
different points in the Brillouin zone. However, if, for example, translation
symmetry is broken (such as by charge density wave order) then nodes at
different momenta can couple to each other and form a correlated insulator.

It has been shown that when an inversion-symmetric WSM is gapped out via a
translation symmetry breaking charge density wave (CDW), the resulting
insulator exhibits axion electrodynamics [17, 18, 19, 20, 21]. The axion
electrodynamics is described by a 3D $\Theta$ term [22, 23, 24], and when the
WSM has inversion symmetry, $\Theta$ is quantized to values of $0$ or $\pi$
[25, 26, 27]. In addition to axion electrodynamics, the inversion-symmetric
Weyl-CDW insulator also exhibits a quantized 3D anomalous Hall response [28,
29, 30]. When inversion-symmetric open boundaries are present, the $\Theta=0$
and $\Theta=\pi$ Weyl-CDW insulators can be distinguished by their total Hall
conductances, which differ by $e^{2}/h$ [21].

DSMs can similarly be gapped out via a CDW to produce a Dirac-CDW insulator
[31]. The resulting insulator can be understood by treating the DSM as two
copies of a WSM related by time-reversal symmetry [6, 7]. This intuition
reveals that the Dirac-CDW insulator cannot display a quantized 3D Hall
response, as the 3D Hall response is odd under time-reversal symmetry. The
Dirac-CDW insulator also has trivial axion electrodynamics (i.e., a vanishing
3D $\Theta$ term) since $\Theta$ is defined $\mod(2\pi)$. Despite the absence
of these effects, there is more to the story of the Dirac-CDW insulator
topological responses. In this work, we consider the crystalline-
electromagnetic (CEM) responses of Dirac-CDW insulators. Similar to how
electromagnetic responses describe charge fluctuations induced by
electromagnetic fields, crystalline-electromagnetic responses describe how
charge fluctuations are induced by lattice deformations (e.g., shears,
strains, and defects) [32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
46, 47, 48, 49, 50].

Since Dirac-CDW insulators have rotation symmetry, they inherit a natural
sensitivity to disclination defects, which are line-like fluxes of rotation
symmetry. However, the exact nature of the CEM response to disclinations
depends on: (i) if the CDW is in/commensurate with the lattice and (ii) if the
Dirac-CDW insulator has inversion symmetry. When the CDW is commensurate with
the lattice, the Dirac-CDW insulator has residual discrete translation
symmetry determined by the CDW period. In this case we will show that
disclination loops in the 3D Dirac-CDW insulator bind a quantized charge per
length. The charge per length is determined by a combination of the Frank
angle of the disclination and the length of the CDW period. This response can
be considered as a CEM analog of the 3D Hall effect, as the former describes
how charge per length is bound to disclination lines, while the latter
describes how charge per length is bound to magnetic flux lines. In fact, the
disclination response is essentially a layered version of the 2D discrete Wen-
Zee response that describes how charge is bound to point-like disclinations of
2D lattices [51, 52, 53, 54, 55, 46, 56, 47, 57, 58, 59]. Hence, we refer to
this response of the 3D Dirac-CDW insulator as the 3D discrete Wen-Zee (dWZ)
response. Furthermore, we will show that in the effective response theory, the
3D dWZ response is captured by a layered version of the 2D Wen-Zee term.

Just like the case for Weyl-CDW insulators, we find that there are two
distinct classes of Dirac-CDW insulators when inversion symmetry is preserved.
The difference between the two classes of insulators manifest as a difference
in the parity of charge bound to disclination lines that terminate on the open
boundaries of a system. In analogy to axion electrodynamics, we find that the
effective response theories for the two classes of Dirac-CDW insulators differ
by a quantized total derivative term. This term, called the $R\wedge F$ term,
was previously discussed by some of us in Ref. 60. Similar to the 3D $\Theta$
term that describes axion electrodynamics, the $R\wedge F$ term is a total
derivative term that leads to anomalous boundary physics that we discuss in
detail below.

Generically we show that when the CDW is commensurate and preserves the
inversion symmetry of the DSM, it is possible to define a “disclination
filling anomaly” [61, 62]. A disclination filling anomaly occurs when it is
impossible to change from periodic to (gapped) open boundary conditions while
preserving both inversion symmetry and conserving charge on the disclination
line. Alternatively put, we show that when the disclination filling anomaly is
present, disclination lines of an inversion-symmetric insulator must bind
different amounts of charge for open and periodic boundary conditions. This
difference in charge is quantized and cannot be removed by boundary effects,
provided they do not break inversion symmetry. This serves to further classify
crystalline insulators by considering filling anomalies in the presence of
topological defects such as a disclination.

The rest of this paper is organized as follows: in Section II we provide an
overview of Dirac semimetals, including lattice models, protecting crystalline
symmetries, and their unquantized quasi-topological responses. In Section IV
we describe how adding a CDW to the DSM model can give rise to correlated
topological crystalline insulators. In Section III we discuss relevant
topological response terms, such as the discrete Wen-Zee and $R\wedge F$
terms. Equipped with this knowledge, in Section V we analytically and
numerically demonstrate the topological crystalline-electromagnetic responses
of Dirac-CDW insulators. Finally, in Section VI we conclude with a discussion
and outlook, and include technical details in the appendices.

##  II Dirac semimetals: lattice model, topology, and responses

In this section we will review the relevant properties of DSMs. We will pay
particular attention to the symmetries that protect the DSMs. We will also
discuss the anomalous responses of the gapless DSM in the presence of probe
gauge fields. Physically, we will demonstrate that these anomalous responses
manifest as charge bound to disclination defects of the lattice.

###  II.1 Lattice model

In this subsection we briefly review some of the band theory of DSMs.
Generically, DSMs are classified into two groups: symmorphic DSMs whose low-
energy gapless features are protected by symmorphic symmetries and translation
symmetry; and non-symmorphic DSMs where Dirac nodes arise from non-symmorphic
symmetry-enforcement of band crossings with four-dimensional irreducible
representations at high-symmetry points at the Brillouin zone edge [1, 7, 63,
6]. We will focus on symmorphic DSMs, and study only a minimal modle hosting a
pair of Dirac nodes that are protected by time-reversal, inversion, rotation,
and translation symmetry. Such DSMs are known to possess unquantized anomalous
responses [64, 42].

Our starting point for the symmorphic DSM will be the following Bloch
Hamiltonian for spin-1/2 fermions on an orthorhombic lattice [7],

$\begin{split}\mathcal{H}_{\text{DSM}}(\bm{k})&=\sin(k_{x}a_{x})\Gamma_{1}+\sin(k_{y}a_{y})\Gamma_{2}\\\
&+[\cos(k_{z}a_{z})-\cos(\mathcal{K}a_{z})]\Gamma_{3}\\\
&-b_{xy}[2-\cos(k_{x}a_{x})-\cos(k_{y}a_{y})]\Gamma_{3},\end{split}$ (1)

where $a_{x,y,z}$ are the lattice constants in the $x$, $y$, and $z$
directions, respectively. The $\Gamma_{a}$ matrices are defined as

$\begin{split}&\Gamma_{1}=\sigma^{x}s^{z},\phantom{=}\Gamma_{2}=-\sigma^{y},\phantom{=}\Gamma_{3}=\sigma^{z},\\\
&\Gamma_{4}=\sigma^{x}s^{x},\phantom{=}\Gamma_{5}=\sigma^{x}s^{y}.\end{split}$
(2)

where $\bm{\sigma}$ and $\bm{s}$ are Pauli matrices acting on the sublattice
and spin degrees of freedom, respectively. We note that this identification
will not be important for our further discussion. Here and throughout we leave
any $2\times 2$ identity matrices implicit, (i.e.,
$\Gamma_{2}=-\sigma^{y}s^{0}$, and $\Gamma_{3}=\sigma^{z}s^{0}$, where $s^{0}$
is the $2\times 2$ identity in spin space).

The (doubly-degenerate) energy bands of this model are given by

$\begin{split}E_{\pm}(\bm{k})=&\pm\Big{[}\sin^{2}(k_{x}a_{x})+\sin^{2}(k_{y}a_{y})\\\
&+\big{(}\cos(k_{z}a_{z})-\cos(\mathcal{K}a_{z})\\\
&-b_{xy}[2-\cos(k_{x}a_{x})-\cos(k_{y}a_{y})]\big{)}^{2}\Big{]}^{1/2}.\end{split}$
(3)

From this equation we see that the parameter $\mathcal{K}$ determines the
location of the Dirac nodes at which the bands become four-fold degenerate:
$\bm{k}_{\text{DN}}\equiv(0,0,\mathcal{K}),$
$E_{\pm}(\bm{k}_{\text{DN}})=E_{\pm}(-\bm{k}_{\text{DN}})=0$. Intuitively,
$\mathcal{K}$ parameterizes a process wherein two pairs of doubly-degenerate
bands are inverted. At $\mathcal{K}=0$ the bands are not inverted and have a
quadratic band touching at $\bm{k}=0$. At $\mathcal{K}=\pi$ the bands are
fully inverted and have a quadratic band touching at $\bm{k}=(0,0,\pi)$. For
other values of $\mathcal{K}$, the two pairs of bands are partially inverted
and have linear crossings at $\bm{k}=(0,0,\pm\mathcal{K})$. Additionally, the
parameter $b_{xy}$ determines the gap away from the high symmetry line
$k_{x}=k_{y}=0$. As expected, the Bloch Hamiltonian near the four-fold
degenerate nodal points takes on the Dirac-form

$\begin{split}\mathcal{H}_{\text{DSM}}(\pm\bm{k}_{\text{DN}}+\bm{q})&\approx
q_{x}a_{x}\Gamma_{1}+q_{y}a_{y}\Gamma_{2}\\\ &\pm
q_{z}a_{z}\sin(\mathcal{K}a_{x})\Gamma_{3}.\end{split}$ (4)

Here, and throughout, we will take the system to be half-filled, i.e., two
electrons per unit cell, such that the chemical potential intersects the Dirac
points.

###  II.2 Symmetries

For generic values of $\mathcal{K}$, the Hamiltonian in Eq. 1 consists of two
doubly degenerate bands that meet at a pair of four-fold degenerate Dirac
nodes. The degeneracies are stabilized by the combination of: C4z symmetry
(four-fold rotation symmetry around the z-axis), $T^{z}$ translation symmetry
(translation symmetry in the z-direction), U$(1)$ charge conservation
symmetry, time-reversal symmetry, and inversion symmetry.

C4z rotation symmetry, $T^{z}$ translation symmetry, and U$(1)$ charge
conservation symmetry are all necessary for the Dirac nodes in Eq. 4 to remain
gapless. The C4z operator is given by

$U_{4}=\exp\left(i\frac{\pi}{2}[\tfrac{1}{2}\sigma^{z}s^{z}-s^{z}]\right).$
(5)

where $(U_{4})^{4}=-1$, since we have spin-$1/2$ fermions. This symmetry
prohibits mass terms $\propto\Gamma_{4}$, and $\propto\Gamma_{5}$. The $T^{z}$
translation symmetry prevents the Dirac node at $+\bm{k}_{\text{DN}}$ from
hybridizing with the Dirac node at $-\bm{k}_{\text{DN}}$. The U$(1)$ charge
symmetry prevents the inclusion of superconducting pairing terms.

The two-fold degeneracy of the bands of Eq. 1 away from the Dirac nodes is
ensured by the combination of time-reversal and inversion symmetry: time-
reversal symmetry guarantees that a state with momentum $\bm{k}$ is
accompanied by a degenerate state at $-\bm{k}$ with opposite spin
($\bm{\sigma}$) projection. Similarly, inversion symmetry guarantees that for
a state with momentum $\bm{k}$, there is a degenerate state at $-\bm{k}$ with
the same spin ($\bm{s}$) projection. The combination of time-reversal and
inversion symmetries therefore guarantees that all states are at least two-
fold degenerate. Furthermore, time-reversal and inversion symmetry also
protect the four-fold degeneracy of the Dirac nodes at
$\bm{k}=\pm\bm{k}_{\text{DN}}$, and prevent the DSM phase from being deformed
into a Weyl semimetal phase. For example, adding a time-reversal breaking
term, such as a magnetic Zeeman term $\propto s^{z}$, will split a Dirac nodes
into one positive-chirality and one negative-chirality Weyl node. If inversion
symmetry is preserved, then this splitting of Dirac nodes must occur in pairs,
and if there is a positive chirality Weyl node at $\bm{k}$, there will be a
negative chirality Weyl node at $-\bm{k}$. This may be understood via the
requirement that the Berry curvature satisfies
$\mathbf{\Omega}(\bm{k})=-\mathbf{\Omega}(-\bm{k})$ because of inversion
symmetry. This Dirac node splitting procedure results in a magnetic Weyl
semimetal that possesses an anomalous Hall conductivity. Alternatively, if
inversion symmetry is broken and time-reversal symmetry is preserved, we may
arrive at a time-reversal invariant Weyl semimetal, where a positive chirality
Weyl node at $\bm{k}$, is accompanied by a second positive chirality Weyl node
at $-\bm{k}$. This naturally arises as a consequence of time-reversal symmetry
which implies $\mathbf{\Omega}(\bm{k})=\mathbf{\Omega}(-\bm{k})$. Analogous to
the magnetic WSM case, such a material is known to possess an unquantized
momentum anomaly response that describes non-trivial charge responses to
dislocations [65, 42, 66, 50] .

The time-reversal operator is given by $\mathcal{T}=-is^{y}\kappa$, where
$\kappa$ is the complex conjugation operator ($\mathcal{T}^{2}=-1$ for
spin-$1/2$ fermions). For the lattice model in Eq. 1 there are multiple
inequivalent choices of inversion symmetry, which correspond to different
choices of inversion center. The different inversion centers are related by a
half-translation in a given direction. Since we are going to consider breaking
translation symmetry in the $z$-direction, we will pay special attention to
the two inversion symmetry definitions that differ by a half-translation in
the $z$-direction:

$\begin{split}&\mathcal{I}_{s}=\sigma^{z},\\\
&\mathcal{I}_{b}=e^{ik_{z}a_{z}}\sigma^{z}.\end{split}$ (6)

We will call $\mathcal{I}_{s}$ the “site-centered” inversion symmetry since it
sends $\bm{r}=(r_{x},r_{y},r_{z})\rightarrow(-r_{x},-r_{y},-r_{z})$ and leaves
the site $(0,0,0)$ invariant. Similarly, we will call $\mathcal{I}_{b}$ the
“bond-centered” inversion symmetry since it sends
$(r_{x},r_{y},r_{z})\rightarrow(-r_{x},-r_{y},-r_{z}+a_{z})$ and leaves the
bond $(0,0,a_{z}/2)$ invariant. There are other choices of inversion symmetry
that differ by half a translation in either the $x$ or $y$-directions from
those above. However, these other inversion symmetry definitions are
unimportant for our forthcoming analysis.

The distinction between the site- and bond-centered inversion symmetries has
important implications when considering systems with open boundary conditions.
Namely, for a system of length $L_{z}=a_{z}N_{z}$ in the $z$-direction (i.e.,
$N_{z}$ lattice sites along the $z$-direction) and open boundary conditions,
site-centered inversion symmetry necessitates that $N_{z}$ is odd, while bond-
centered inversion symmetry necessitates that $N_{z}$ is even. We will discuss
the interplay between inversion symmetry and CDWs in Sec. IV below.

###  II.3 Anomalous topological response

As shown in Ref. [42], the DSM in Eq. 1 possess an anomalous topological
response to fluxes of C4z rotation symmetry. This response can be expressed in
terms of probe (non-dynamic) gauge fields for the U$(1)$ charge symmetry
($A_{\mu}$) and C4z rotation symmetry ($\omega_{\mu}$) as

$\mathcal{L}_{\text{anom}}=\frac{\nu}{2\pi^{2}}G_{z}\epsilon^{ijk}\omega_{i}\partial_{j}A_{k},$
(7)

where $i,j,k$ run over $t,x,y$, $\nu=2\mathcal{K}/2\pi$, and
$G_{z}=2\pi/a_{z}$ is the reciprocal lattice vector along the $z$ direction.
The origin of such a term may be understood via the lowest Landau level (LLL)
picture wherein a uniform magnetic flux is turned on in the $z$-direction:
such a system possesses two LLLs with different C4z charges, giving rise to a
filling-type anomaly that protects the bands from gapping out [42]. We also
note that Eq. 7 is not gauge invariant for arbitrary values of $\nu$, and
gauge invariance is restored by the gapless degrees of freedom at the nodes of
the DSM.

Figure 1: (a) A lattice with C4z symmetry and no disclinations. Upon cutting
out the greyed-out quadrant and re-gluing, we arrive at (b) a disclinated
lattice. (c) A stack of disclinations of 2D layers forms a disclination line
of a 3D crystal, as shown in (d).

Fluxes of the rotation gauge field, $\omega$, encode disclinations with a
Frank-vector $\propto{\hat{z}}$. In 3D lattices, the disclinations are line-
defects that extend parallel to their Frank-vector. Because of this, we will
implicitly take all disclinations to stretch along the $z$-direction in this
work. One can think of these disclination lines as stacks of 2D disclinations
of the $xy$-planes, as shown in Fig. 1. At the level of gauge field
configurations, considering disclinations that stretch along only the
$z$-direction amounts to the condition that
$\oint_{C_{xz}}\omega=\oint_{C_{yz}}\omega=0$, where the loop integrals are
over curves $C_{ab}$ in the $ab$-plane that do not intersect the cores of any
disclinations.

While the (discrete) gauge field $\omega$ is formally characterized only by
its holonomies, previous works [60, 57] have shown that it is possible to
treat $\omega$ as a continuous gauge field, but having only quantized fluxes,
$\oint\omega\in[0,\frac{\pi}{2},\pi,\frac{3\pi}{2})$. This flux quantization
can be thought of as arising from a Higgs mechanism, similar to the
electromagnetic flux quantization in superconductors, where the analogy of the
electric charge is replaced by angular momentum. Here and throughout, we shall
treat $\omega$ as smoothly varying, but with quantized fluxes.

Now that we have discussed the probe fields we can interpret the anomalous
response. We find that Eq. 7 implies that disclinations in a DSM bind charge
per length determined by $\nu$. To demonstrate, let us take a functional
derivative of the action generated by Eq. 7 with respect to $A_{0}$. Hence,
the contribution to the charge density from the DSM anomalous response is

$\rho_{\text{anom}}=\frac{\nu}{2\pi^{2}}G_{z}[\partial_{x}\omega_{y}-\partial_{y}\omega_{x}].$
(8)

For a thin-core disclination parallel to the $z$-axis,
$\partial_{x}\omega_{y}-\partial_{y}\omega_{x}=\Theta_{F}\delta(x)\delta(y)$,
where $\Theta_{F}$ is the Frank angle and is quantized to $\Theta_{F}=2\pi
s/4$ with $s\in\mathbb{Z}$. For a system of length $L_{z}=a_{z}N_{z}$ with
$N_{z}\in\mathbb{Z}$, the disclination has total charge

$\begin{split}Q_{\text{disc}}=\frac{\Theta_{F}}{\pi}\nu N_{z}+\delta
Q.\end{split}$ (9)

Here we have assumed that the system is charge neutral in the absence of any
disclinations. Physically, this amounts to adding negatively-charge background
ions to each site of the lattice. In addition to the topological contribution
$\propto\Theta_{F}$, we have included an additional contribution to the
disclination charge, $\delta Q$. This contribution is $\mathcal{O}(1)$ and
arises from the gapless particles that are necessary to restore the gauge
invariance of the anomalous DSM response [42]. In Fig. 2(a) we plot the
disclination charge as a function of $\nu$ using the DSM tight binding model
in Eq. 1, and find that the disclination charge follows the predicted trend,
$Q_{\text{disc}}/N_{z}\approx\frac{\Theta_{F}}{\pi}\nu$.

A useful way of understanding the disclination response of the DSM is to treat
the momentum $k_{z}$ of Eq. 1 as a tunable parameter, such that gapped planes
with fixed $k_{z}$ are trivial 2D insulators for $|k_{z}|>\mathcal{K}$ and
quantum spin Hall (QSH) insulators for $|k_{z}|<\mathcal{K}$ [67, 68]. The
Dirac nodes at $|k_{z}|=\mathcal{K}$ mark the band crossings that connect the
topologically distinct insulators. In this way, the DSM with $\mathcal{K}\neq
0,\pi$ can be viewed as an intermediate phase between a trivial insulator
($\mathcal{K}=0$), and a weak 3D topological insulator ($\mathcal{K}=\pi$)
[69, 64]. As we shall discuss in Sec. III.1.1, disclinations of the 2D QSH
insulator bind charge $\Theta_{F}/\pi$. For a system of length
$L_{z}=a_{z}N_{z}$, $k_{z}$ is discretized in steps of size $2\pi/N_{z}$. Each
$k_{z}$ slice will bind charge $\Theta_{F}/\pi$ to a disclination if
$|k_{z}|<\mathcal{K}$, and charge $0$ if $|k_{z}|>\mathcal{K}$, as shown in
Fig. 2(b). The disclination will therefore bind a total charge

$Q_{\text{disc}}=\sum^{\kappa}_{k_{z}=-\kappa}\frac{\Theta_{F}}{\pi}\approx\frac{\Theta_{F}}{\pi}2\mathcal{K}\frac{N_{z}}{2\pi}=\frac{\Theta_{F}}{\pi}\nu
N_{z}.$ (10)

If we add in the extra $\mathcal{O}(1)$ contribution, $\delta Q$, that arises
from the gapless modes at $k_{z}=\pm\mathcal{K}$, we arrive at Eq. 9.

Figure 2: (a) The numerically computed charge per layer bound to a
$\Theta_{F}=\pi/2$ disclination as a function of $\nu$ for the Hamiltonian Eq.
1 with $b_{xy}=1$ (black circles). The dashed line for comparison has slope
$\frac{1}{\pi}\Theta_{F}\nu$ and vanishing intercept, indicating the
theoretical prediction. The disclinated lattice is constructed as in Ref. [60]
for $N_{x}=N_{y}=39$ lattice sites in the $xy$-plane, and 250 momentum points
along the $k_{z}$ axis. (b) The bound charge resolved in $k_{z}$ (momentum
parallel to the disclination line) for $\mathcal{K}=\pi/3$. The deviations
from 0 and 1/2 near $k_{z}=\pm\frac{\pi}{3}\frac{1}{a_{z}}$ are the result of
finite size effects and the small size of the gap near the Dirac nodes.

##  III Topological crystalline-electromagnetic responses of crystalline
insulators

In preparation for our discussion of the topological crystalline-
electromagnetic responses of the Dirac-CDW insulator in Sec. V, we will first
discuss the crystalline-electromagnetic responses that are relevant for
insulators having C4z rotation symmetry, as well as $T_{z}$ discrete
translation and/or inversion symmetry. These responses are similar to the
anomalous response of the DSM discussed in Sec. II.3. However, unlike the DSM
responses, the topological responses of insulators are quantized because of
gauge invariance and symmetry (DSMs can have unquantzed responses since the
there are low energy gapless degrees of freedom that restore gauge
invariance). We will mainly consider insulators composed of spin-1/2 fermions
with TRS, which is the relevant case for DSMs. Similar responses also exist
for spinless fermions, and systems without TRS, and the quantization of the
responses will be different [46, 60].

For spin-1/2 systems with TRS and C4z rotation symmetry, we will show that
crystalline-electromagnetic responses lead to: a $\mathbb{Z}_{4}$
classification when the system has additional $T^{z}$ discrete translation
symmetry; a $\mathbb{Z}_{2}$ classification when the insulator has additional
inversion symmetry; and a $\mathbb{Z}_{4}\times\mathbb{Z}_{2}$ classification
when the system has both $T^{z}$ discrete translation and inversion symmetry.
We will also show that for a system with discrete translation
$r_{z}\rightarrow r_{z}+a_{z}$, the $\mathbb{Z}_{4}$ index indicates that
disclinations with Frank angle $\Theta_{F}$ bind a charge per length of
$2n\frac{\Theta_{F}}{2\pi a_{z}}$ ($n\in\mathbb{Z}_{4}$) when the system has
periodic boundaries in the $z$-direction. The effective response theory for
the insulator with a non-trivial $\mathbb{Z}_{4}$ index also contains a 3D
discrete Wen-Zee (dWZ) term. The $\mathbb{Z}_{2}$ index corresponds to the
charge parity of a disclination when inversion-symmetric open boundary
conditions are present. Insulators with different $\mathbb{Z}_{2}$ indices
have effective response theories that differ by a quantized total derivative
term that was referred to as the $R\wedge F$ term in Ref. 60. When the
insulator has both $T^{z}$ discrete translation and inversion symmetry, the
$\mathbb{Z}_{2}$ index can be equivalently understood as the presence or
absence of a filling anomaly [61, 62] of the 1D disclination lines, which we
will define more formally in Sec. III.2.

To prevent confusion we note that here we are discussing only the topological
crystalline-electromagnetic responses of insulators in various symmetry
classes. This is not an exhaustive list of all topological responses. Notably,
we are ignoring any purely electromagnetic responses in the section.

###  III.1 The discrete Wen-Zee response

To begin, we will consider the dWZ response in 2D and 3D. The dWZ response of
a 2D system describes how disclinations bind a quantized charge [51, 52, 53,
54, 55, 46, 56, 47, 57, 58, 59]. The 3D dWZ response is essentially a layered
version of the 2D dWZ response, and describes how disclination lines bind a
quantized charge per length.

####  III.1.1 The discrete Wen-Zee response in 2D

The 2D dWZ response corresponds to the following topological term in the
effective field theory,

$\mathcal{L}_{\text{2D-WZ}}=\frac{\mathcal{S}_{\text{2D}}}{2\pi}\epsilon^{ijk}\omega_{i}\partial_{j}A_{k},$
(11)

where $i,j,k$ run over $t,x,y$. Here, fluxes of $\omega$ (which are point-like
in 2D) are disclinations of the 2D lattice. For spin-$1/2$ fermions with TRS
(which will be the focus of this work) $\mathcal{S}_{\text{2D}}$ is quantized
to be an even integer 111For insulators composed of spinless fermions with TRS
or spin-1/2 fermions without TRS, $\mathcal{S}_{\text{2D}}$ is quantized as an
integer. For insulators composed of spinless fermions without TRS,
$\mathcal{S}_{\text{2D}}=C/2\mod(1)$ where $C$ is the Chern number. and
defined $\mod(8)$ in the case of insulators with C4 symmetry (four-fold
rotations of the plane) 222For spin-1/2 fermions with TRS,
$\mathcal{S}_{\text{2D}}$ is defined $\mod(2n)$ in $C_{n}$ symmetric
insulators. For spinless fermions either with or without TRS, or spin-1/2
fermions without TRS, $\mathcal{S}_{\text{2D}}$ is defined $\mod(n)$ in
$C_{n}$ symmetric insulators. Therefore the allowed inequivalent values are
$\mathcal{S}_{\text{2D}}=0$, $2$, $4$, or $6$ for time-reversal and C4
symmetric insulators.

A representative 2D insulator with $\mathcal{S}_{\text{2D}}=2$, is realized by
the following 4-band Hamiltonian,

$\begin{split}\mathcal{H}_{\text{QSH}}(k_{x},k_{y})&=\sin(k_{x})\Gamma_{1}+\sin(k_{y})\Gamma_{2}\\\
&-[m-\cos(k_{x})-\cos(k_{y})]\Gamma_{3},\end{split}$ (12)

with C4 rotation given by Eq. 5, and TRS given by $\mathcal{T}=-is^{y}\kappa$.
This Hamiltonian is a QSH insulator for $0<|m|<2$ and a trivial insulator for
$2<|m|$ [67, 68]. The QSH insulator has a 2D dWZ response with
$\mathcal{S}_{\text{2D}}=2$, while the 2D dWZ response vanishes
($\mathcal{S}_{\text{2D}}=0$) in the trivial phase333In Ref. [59], a trivial
insulator has a non-zero discrete shift. This difference arises because Ref.
[59] considers the charge contribution from only the electrons. Here we are
considering the contribution from the electrons, and the contribution from the
background ions, which are effectively a trivial insulator composed of
negative charges. [57]. We note that Eq. 12 can be extended to Eq. 1 by fixing
$b_{xy}=1$ and setting $M=2+\cos(k_{z})-\cos(Q).$

The 2D dWZ response indicates that a disclination with Frank angle
$\Theta_{F}$ binds charge

$Q_{\text{disc}}=\mathcal{S}_{\text{2D}}\frac{\Theta_{F}}{2\pi}\mod(2),$ (13)

where we have added additional negative background charges such that the
system is charge neutral in the absence of disclinations. The $\mod(2)$
ambiguity in defining $Q_{\text{disc}}$ arises from the fact that for spin-1/2
insulators with TRS, it is possible to add a Kramers pair of particles to the
disclination core without changing any topological properties. Since
$\Theta_{F}$ is a multiple of $\pi/2$, the disclination charge is always
trivial if $\mathcal{S}_{\text{2D}}=8$. This is why $\mathcal{S}_{\text{2D}}$
is defined $\mod(8)$ as noted before.

####  III.1.2 The discrete Wen-Zee response in 3D

For insulators having $T^{z}$ discrete translation, $z\rightarrow z+a_{z}$,
the 2D dWZ response can be extended to 3D as a layered response. The 3D
layered response can be expressed as the following topological term in the
effective response theory,

$\mathcal{L}_{\text{3D-WZ}}=\mathcal{S}_{\text{3D}}\frac{G_{z}}{4\pi^{2}}\epsilon^{ijk}\omega_{i}\partial_{j}A_{k},$
(14)

where $G_{z}=\frac{2\pi}{a_{z}}$ is the reciprocal lattice vector along the
$z$-direction. Here, as in Sec. II, flux-lines of $\omega$ are disclination
lines with Frank-vector $\propto{\hat{z}}$. Again, we take all disclination
lines to extend along the $z$-direction. This response is a layered version of
the 2D dWZ response having one 2D layer per lattice period, $a_{z}$. The
coefficient $\mathcal{S}_{\text{3D}}$, which we will call the 3D discrete
shift, is defined $\mod(8)$ and quantized as an even integer for spin-$1/2$
fermions with TRS and C4z rotation symmetry, i.e. $\mathcal{S}_{\text{3D}}=0$,
$2$, $4$, or $6$. Because of the this, the value of $\mathcal{S}_{\text{3D}}$
defines a $\mathbb{Z}_{4}$ index for insulators in this symmetry class.

This quantization of $\mathcal{S}_{\text{3D}}$ is directly inherited from the
quantization of $\mathcal{S}_{\text{2D}}$. The 3D dWZ term, Eq. 14, has the
same form as the unquantized anomaly of the DSM, Eq. 7. However, the
coefficient of the 3D dWZ term, $\mathcal{S}_{\text{3D}}$, is quantized, while
the coefficient of the DSM anomaly equation, $\nu$, is not (recall that $\nu$
is not quantized due to the gaplessness of the DSM).

The 3D dWZ response physically manifests as a quantized charge per length
bound to disclination lines as illustrated in Fig. 3a. If we take a system of
length $L_{z}=a_{z}N_{z}$ that includes a disclination line with Frank angle
$\Theta_{F}$, Eq. 14 indicates that the total charge on the disclination line
is

$\begin{split}Q_{\text{disc}}=\frac{\Theta_{F}}{2\pi}\mathcal{S}_{\text{3D}}N_{z}\mod(2N_{z}).\end{split}$
(15)

Here we have again implicitly added negatively charged background ions, such
that the system is charge neutral without any disclinations. The
$\mod(2N_{z})$ ambiguity arises from the fact that it is possible to change
the disclination charge by embedding a 1D insulator in the disclination core
without changing any topological properties of the system. Due to the Kramers
degeneracy and the Lieb-Schultz-Mattis theorem [70], the charge per unit
length must be an even integer for translation invariant spin-1/2 1D
insulators with TRS. From the $\mod(2N_{z})$ ambiguity, and the fact that the
Frank angles of disclinations of a C4z symmetric insulators are quantized in
units of $\pi/2$, we also find that $\mathcal{S}_{\text{3D}}=8$ leads to a
trivial charge of the disclination, confirming our earlier assertion that
$\mathcal{S}_{\text{3D}}$ is defined $\mod(8)$.

Figure 3: (a) The total charge per length, $Q_{\text{disc}}/N_{z}$, of an
insulator with a disclination with Frank angle $\Theta_{F}$ is proportional to
$\mathcal{S}_{3D}$, the coefficient of the dWZ term. (b) The disclination
charge parity, $Q_{\text{disc}}\mod(2\Theta_{F}/\pi)$, of a
$\Theta_{F}$-disclination of an inversion-symmetric insulator with open
boundary conditions is determined by $\Phi(0)=\\{0,2\pi\\}$, the value of the
coefficient of the $R\wedge F$ term at the inversion center.

Without $T^{z}$ discrete translation symmetry, the charge per length of a
disclination line is not a well-defined quantity, as it is possible to add
local charges to arbitrary points along the disclination. An important
consequence of this is that the charge per length of a disclination is not a
topological, quantized quantity for insulators with incommensurate CDWs, as it
is possible to shift the charge per length of a disclination by an arbitrary
amount. However, for systems that are finite in the $z$-direction, the total
charge on a disclination line $Q_{\text{disc}}$ is still a quantized multiple
of $\Theta_{F}/\pi$. This can be understood by treating the $z$-coordinate as
an internal degree of freedom of an effectively 2D system with TRS and C4z
rotation symmetry. A disclination-line of the finite 3D system is also a
disclination of the effective 2D system, and the total charge bound to the
disclination is therefore quantized according to the 2D dWZ response (Eq. 11).

###  III.2 The $R\wedge F$ Term

In Ref. 60 some of us showed that the effective response theory of 3D
crystalline insulators with rotation symmetry around a fixed axis can contain
the following topological term,

$\begin{split}\mathcal{L}_{R\wedge
F}&=\epsilon^{\mu\nu\lambda\eta}\frac{\Phi}{4\pi^{2}}\partial_{\mu}\omega_{\nu}\partial_{\lambda}A_{\eta}.\end{split}$
(16)

where $\mu,\nu,\lambda,\eta$ run over $t,x,y,z,$ and $\Phi$ can generally be a
function of position and time. We denote this response as the $R\wedge F$
term, since it couples the lattice curvature
$R_{\mu\nu}=\partial_{\mu}\omega_{\nu}-\partial_{\nu}\omega_{\mu}$ to the
electromagnetic field strength
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ [60]. The $R\wedge F$
term has a similar form to the $\Theta$ term, $F\wedge F$, that describes
axion electrodynamics in 3D topological insulators [17, 22]. Like the $\Theta$
term, the $R\wedge F$ term has a periodic coefficient and is a total
derivative when the coefficient is constant. For spin-$1/2$ insulators with
TRS, $\Phi$ is $4\pi$ periodic. In general, the periodicity of $\Phi$ depends
on the spin of the fermions and the presence of TRS, as discussed in Ref. 60.

The charge responses associated with the $R\wedge F$ term are

$j^{\mu}=\epsilon^{\mu\nu\lambda\eta}\frac{1}{4\pi^{2}}(\partial_{\nu}\Phi)\partial_{\lambda}\omega_{\eta}.$
(17)

If we fixed $\omega$ to the configuration of a disclination line,
$\partial_{x}\omega_{y}-\partial_{y}\omega_{x}=\Theta_{F}\delta(x)\delta(y)$,
the charge responses simplify to

$\begin{split}&\rho_{\text{disc}}=\frac{\Theta_{F}}{4\pi^{2}}(\partial_{z}\Phi)\delta(x)\delta(y),\\\
&j^{z}_{\text{disc}}=-\frac{\Theta_{F}}{4\pi^{2}}(\partial_{t}\Phi)\delta(x)\delta(y).\end{split}$
(18)

The first lines indicates that spatial fluctuations of $\Phi$ bind charge
along disclination lines. The second line indicates that temporal fluctuations
drive a current along the disclination line. Together, these two responses
indicate a charge polarization in the $z$-direction proportional to
$\Theta_{F}$ localized on the disclination line.

Based on Eq. 18, if we adiabatically increase $\Phi$ homogeneously by an
amount $\delta\Phi$ over some time $T$ (i.e.,
$\Phi(t)=\Phi_{0}+\delta\Phi\frac{t}{T}$), then the polarization on the
disclination line, $P^{z}_{\text{disc}}$, will change by

$\begin{split}\Delta
P^{z}_{\text{disc}}&=\int^{T}_{0}\\!dt\partial_{t}P^{z}_{\text{disc}}=\int^{T}_{0}\\!dtj^{z}_{\text{disc}}=\frac{\delta\Phi}{2\pi}\frac{\Theta_{F}}{2\pi}.\end{split}$
(19)

We therefore find that if the effective response theories of two insulators
differ by an $R\wedge F$ term with a constant coefficient $\delta\Phi$, the
polarization of disclination lines will differ by
$\frac{\delta\Phi}{2\pi}\frac{\Theta_{F}}{2\pi}$ between the two insulators.
444Note that two systems can have a well-defined difference in polarization
even if the systems do not have a well-defined polarization individually. For
example, a system with a net charge does not have a well-defined polarization,
as the polarization will depend on the choice of origin. However, the
difference in polarization between two systems with the same net charge is
independent of the choice of origin.

####  III.2.1 $R\wedge F$ term with discrete translation symmetry

We now turn our attention to the connection between the $R\wedge F$ term and
the 3D dWZ response. We will consider a system having $T^{z}$ discrete
translation, $z\rightarrow z+a_{z}$, and, since the disclination responses in
Eq. 18 primarily involve the dependence of $\Phi$ on $z$ and $t$, we will
suppress any dependence on $x$ and $y$ here. Interestingly, the coefficient of
the $R\wedge F$ term can be non-constant for a system with $T^{z}$ discrete
translation symmetry, as $\Phi(z)=\Phi(z+a_{z})\mod(4\pi)$ has non-constant
solutions. We can characterize the solutions by their winding across a single
unit-cell period $\int^{a_{z}}_{0}dz\partial_{z}\Phi(z)$, which will be an
integer multiple of $4\pi$. Examples of different winding configurations are
shown in Fig. 4.

With this in mind we can compute the total charge bound to a disclination for
a periodic system of length $L_{z}=a_{z}N_{z}$ (using Eq. 18):

$\begin{split}Q_{\text{disc}}&=\frac{\Theta_{F}}{2\pi}\int^{L_{z}}_{0}\partial_{z}\frac{\Phi(z)}{2\pi}dz\mod(2N_{z})\\\
&=\frac{\Theta_{F}}{2\pi}N_{z}\int^{a_{z}}_{0}\partial_{z}\frac{\Phi(z)}{2\pi}dz\mod(2N_{z}),\end{split}$
(20)

where the $\mod(2N_{z})$ ambiguity again arises from the fact that it is
possible to change the disclination charge per length by embedding a 1D
insulator in the disclination core. Comparing to Eq. 15, we find a relation
for the dWZ response coefficient:

$\mathcal{S}_{\text{3D}}=\frac{1}{2\pi}\int^{a_{z}}_{0}\\!dz\partial_{z}\Phi(z),$
(21)

from which we see that the winding of $\Phi$ over one period, $a_{z}$,
generates a 3D dWZ response.

We can make the connection between the $R\wedge F$ and 3D dWZ terms more
concrete by setting

$\Phi(z)=\Phi_{0}+\mathcal{S}_{\text{3D}}G_{z}z,$ (22)

where $\Phi_{0}$ is an arbitrary constant. This solution preserves the $T^{z}$
discrete translation symmetry when $G_{z}=2\pi/a_{z}$. Plugging this value of
$\Phi(z)$ into Eq. 16 leads to the 3D dWZ term in Eq. 14, after an integration
by parts. Because of the $\mod(2N_{z})$ ambiguity, the $R\wedge F$ term
identifies a $\mathbb{Z}_{4}$ classification of insulators with $T^{z}$
discrete translation symmetry, exactly as discussed in Sec. III.1.2.

Figure 4: Two examples of how $\Phi$ can wind across the unit cell. Both lines
wind by $-4\pi$ across the unit cell, but one has $\Phi(0)=0$ (red dashed) and
the other has $\Phi(0)=2\pi$ (black solid).

Although it is not manifest in Eq. 22, it is important to note that the
coefficient of the $R\wedge F$ term, $\Phi$, is a $4\pi$-periodic variable.
Importantly, when $\Phi$ has a winding there is no discontinuity between $z=0$
and $z=L_{z}$, since $\mathcal{S}_{\text{3D}}$ is an even integer and
$\mathcal{S}_{\text{3D}}G_{z}L_{z}=\mathcal{S}_{\text{3D}}2\pi
N_{z}=0\mod(4\pi)$. The expression for $\Phi(z)$ given in Eq. 22 is therefore
smooth everywhere with $\partial_{z}\Phi(z)=\mathcal{S}_{\text{3D}}G_{z}$. The
same would be true if we replaced the linear interpolation in Eq. 22 with a
more complicated function like those shown in Fig. 4.

####  III.2.2 $R\wedge F$ term with inversion symmetry

Let us now consider the $R\wedge F$ term for a system with inversion symmetry.
We again take $\Phi$ to be independent of the $x$ and $y$ coordinates but
allow for $\Phi$ to depend on $z$. We will now show that we can use the
$R\wedge F$ term to resolve a $\mathbb{Z}_{2}$ distinction between inversion-
symmetric insulators. Inversion symmetry requires that
$\Phi(z)=-\Phi(-z)\mod(4\pi)$, i.e., $\Phi$ must be an odd function of $z$.
Here we take the inversion center to be in the $z=0$ plane. It is also
possible to have an inversion center at $z=a_{z}/2$, but the difference
between the two inversion centers can be accounted for by a redefinition of
the coordinate system. In continuum effective field theory, such a coordinate
change is innocuous. However, such a change is not innocuous for lattice
systems, as the discrete lattice sites lead to a unique choice of origin
(modulo lattice translations). At $z=0$, inversion symmetry requires that
$\Phi(0)=-\Phi(0)$. Since $\Phi$ is $4\pi$ periodic, $\Phi(0)=0$ or
$\Phi(0)=2\pi$ are both allowed values satisfying this condition. The choice
of these values defines two classes (i.e., a $\mathbb{Z}_{2}$ classification)
of inversion-symmetric insulators. Alternatively expressed, the response
theories of insulators with different $\mathbb{Z}_{2}$ indices differ by a
$R\wedge F$ term with $\Phi(0)=2\pi$.

Importantly, the value of $\Phi(0)$ has a direct physical interpretation when
switching from periodic boundaries to open boundaries in the $z$-direction. We
will show below that with open boundaries, a disclination will have charge
$Q_{\text{disc}}=0\mod(2\frac{\Theta_{F}}{\pi})$ when $\Phi(0)=0,$ and charge
$\frac{\Theta_{F}}{\pi}\mod(2\frac{\Theta_{F}}{\pi})$ when $\Phi(0)=2\pi$. We
will refer to the quantity $Q_{\text{disc}}\mod(2\frac{\Theta_{F}}{\pi})$ as
the disclination charge parity. The disclination charge parity is a natural
quantity to consider for inversion-symmetric systems with open boundary
conditions, as purely 2D boundary effects can at most change the disclination
charge by an integer multiple of $2\frac{\Theta_{F}}{\pi}$, provided they
respect inversion symmetry. Hence the charge parity is a physical observable
that determines the $\mathbb{Z}_{2}$ classification.

To show why the value of $\Phi(0)$ determines the disclination charge parity,
let us consider a system of length $L_{z}$ with periodic boundaries in the
$z$-direction, and use local perturbations to “cut” the system and generate
open boundaries at $z=\pm L_{z}/2$. We assume that the resulting open boundary
system satisfies the following four conditions:

  1. 1.

Changing from periodic to open boundary conditions preserves inversion
symmetry around $z=0$.

  2. 2.

The system with open boundaries is gapped, such that $\Phi(z)$ can be defined
over the entire system.

  3. 3.

The value of $\Phi$ at the inversion center, $\Phi(0)$ remains constant when
changing from periodic to open boundaries. This condition is satisfied if the
gap at the inversion center remains open when switching the boundary
conditions.

  4. 4.

At the boundaries of the system, $\Phi$ vanishes (i.e., $\Phi(z)=0$ for
$|z|\geq L_{z}/2$). This choice of boundary conditions is equivalent to
requiring that the $R\wedge F$ term fully vanishes outside the boundaries of
the system.

Let us now add a disclination to the system with open boundary conditions.
Using Eq. 18, the total charge bound to the disclination is

$\begin{split}Q_{\text{disc}}&=\frac{\Theta_{F}}{2\pi}\frac{1}{2\pi}\int^{L_{z}/2}_{-L_{z}/2}\partial_{z}\Phi(z)dz\\\
&=\frac{\Phi(0)}{\pi}\frac{\Theta_{F}}{2\pi}+2n\frac{\Theta_{F}}{\pi},\\\
\end{split}$ (23)

where $n\in\mathbb{Z}$ is the number of times $\Phi$ fully winds by $4\pi$,
between $0$ and $L_{z}/2$, and we have used that $\Phi(z)=-\Phi(-z)$, and
$\Phi(L_{z}/2)=0$. This equation reflects how if $\Phi(0)=\Phi(L_{z}/2)=0$,
$\Phi$ must wind by an integer multiple of $4\pi$ between $0$ and $L_{z}/2$,
while if $\Phi(0)=2\pi$ and $\Phi(L_{z}/2)=0$, $\Phi$ must instead wind by a
half-integer integer multiple of $4\pi$. This is illustrated in Fig. 3b, and
two representative configurations of $\Phi(z)$ are shown in Fig. 5, where
$n=3$. We therefore find that

$Q_{\text{disc}}=\frac{\Phi(0)}{\pi}\frac{\Theta_{F}}{2\pi}\mod(2\frac{\Theta_{F}}{\pi}).$
(24)

The difference between the $\Phi(0)=0$ and $\Phi(0)=2\pi$ insulators is
therefore manifest in the value of
$Q_{\text{disc}}\mod(2\frac{\Theta_{F}}{\pi})$ when open boundaries are
present. Specifically $Q_{\text{disc}}=0\mod(2\frac{\Theta_{F}}{\pi})$ and
$\frac{\Theta_{F}}{\pi}\mod(2\frac{\Theta_{F}}{\pi})$ for $\Phi(0)=0$ and
$2\pi$ respectively. We also find that
$Q_{\text{disc}}=0\mod(2\frac{\Theta_{F}}{\pi})$ for $\Phi(0)=4\pi$, in
agreement with $\Phi$ being $4\pi$ periodic.

Figure 5: The disclination charge parity of inversion-symmetric insulators is
determined by the value of $\Phi(z)$ at the inversion center (the $z=0$ plane
here). Figures (a) and (b) show two windings of $\Phi(z)$ with $n=3$ but
different values of $\Phi(0)$. Figure (a) has $\Phi(0)=2\pi$ and thus winds a
total of seven times, producing an odd disclination charge parity. Figure (b)
has $\Phi(0)=0$ and therefore winds only six times, yielding an even
disclination charge parity.

####  III.2.3 The $R\wedge F$ term with discrete translation and inversion
symmetry

If an insulator has both inversion symmetry and $T^{z}$ discrete translation
symmetry ($z\rightarrow z+a_{z}$), then the values of $\Phi(0)$ and
$\mathcal{S}_{\text{3D}}$ from Eq. 21 lead to a
$\mathbb{Z}_{4}\times\mathbb{Z}_{2}$ classification based on the phenomena
associated to the $R\wedge F$ term. We stress that this classification is not
exhaustive, but it is sufficient to characterize the Dirac-CDW insulators we
will discuss in the following sections.

If two insulators have the same 3D dWZ response (i.e., the same value of
$\mathcal{S}_{\text{3D}}$ in Eq. 21), but different values of $\Phi(0)$, the
polarization of disclination lines will differ by
$\frac{\delta\Phi(0)}{2\pi}\frac{\Theta_{F}}{2\pi}$, where the difference
between the value of $\Phi(0)$ between the two insulators, $\delta\Phi(0)$, is
equal to either $0$ or $2\pi$ based on our previous discussions. To show why
different values of $\Phi(0)$ lead to a difference in polarization, we take a
given $\Phi(z)$ that is invariant under $T^{z}$ discrete translation and
inversion symmetry, and define an adiabatic evolution where $\Phi(z)$
increases by $2\pi$ over a time $T$, $\Phi(z)\rightarrow\Phi(z)+2\pi t/T$.
Inversion symmetry is broken for $0<t<T$ but is restored at $t=T$. During this
process, $\Phi(0)$ changes by $2\pi$, and
$\int^{a_{z}}_{0}\\!dz\partial_{z}\Phi(z)$ remains constant (i.e.,
$\mathcal{S}_{\text{3D}}$ remains constant). Based on Eq. 18, and 19 the
polarization of a disclination line changes by $\frac{\Theta_{F}}{2\pi}$
during this process. Recall that disclinations of insulators with a non-zero
$\mathcal{S}_{\text{3D}}$ have a finite charge per length. Therefore, the
polarization of a disclination is not strictly well-defined for insulators
with finite $\mathcal{S}_{\text{3D}}$. However, the difference in the
polarization of disclinations between insulators with the same
$\mathcal{S}_{\text{3D}}$ is well-defined. This is analogous to discussions of
polarization in Chern insulators[71, 72].

For insulators having $\mathcal{S}_{\text{3D}}=2$, the two classes of
inversion-symmetric insulators are related to each other by a half-translation
$z\rightarrow z+a_{z}/2$. For $\mathcal{S}_{\text{3D}}=2$, this is clear if we
set $\Phi(z)=2G_{z}z$, for which $\Phi(0)=0$. Translating the system by
$a_{z}/2$ shifts $\Phi(z)\rightarrow 2\pi+2G_{z}z$, for which $\Phi(0)=2\pi$.
This relationship can be naturally understood if we treat the 3D dWZ response
as arising from a stacking of 2D systems having $\mathcal{S}_{\text{2D}}=2$.
In order to preserve inversion symmetry around $z=0$, the 2D systems must be
centered at either $z=na_{z}$ or $z=na_{z}+a_{z}/2$ with $n\in\mathbb{Z}$. A
disclination line of such a system, will have charge $\Theta_{F}/\pi$ centered
at $z=na_{z}$ or $z=na_{z}+a_{z}/2$ respectively. This is the same charge
configurations that are realized by the two topologically distinct phases of
the 1D inversion symmetric Su–Schrieffer–Heeger (SSH) model [73, 74], if we
set the electron charge to be $\Theta_{F}/\pi$ instead of $1$ in the SSH
model.

We similarly have that the two insulators with $\mathcal{S}_{\text{3D}}=6$ are
related by a half-translation. This can be understood using that
$\mathcal{S}_{\text{3D}}=6=-2\mod(8)$, and hence the
$\mathcal{S}_{\text{3D}}=6$ dWZ response is captured by $\Phi(z)=-2G_{z}z$. A
half translation then shifts
$\Phi(0)\rightarrow\Phi(0)-2\pi=\Phi(0)+2\pi\mod(4\pi)$. We can also invoke a
similar stacking argument to the one used above, but instead stacking 2D
systems with $\mathcal{S}_{\text{2D}}=-2$.

The two inversion symmetric insulators having $\mathcal{S}_{\text{3D}}=4$ are
related to each other by a quarter translation, $z\rightarrow z+a_{z}/4$.
Here, we can imagine the 3D dWZ response as arising from having two 2D systems
between $z=na_{z}$ and $z=(n+1)a_{z}$, each with $\mathcal{S}_{\text{2D}}=2$.
To preserve inversion symmetry, these 2D systems must be stacked in one of two
classes of configurations. In the first class of stacking configurations, one
insulator is at $z=na_{z}+\delta_{z}$ and one is at $z=(n+1)a_{z}-\delta_{z}$
for each $n\in\mathbb{Z}$, where $\delta_{z}\in[0,a_{z})$ is a constant
offset. Since inversion symmetry is preserved for all $\delta_{z}$ here,
stacking configurations having different values of $\delta_{z}$ are
adiabatically connected. In the second class of stacking configurations, one
insulator is at $z=na_{z}$ and one is at $z=na_{z}+a_{z}/2$. It can be
directly confirmed that the second stacking configuration is distinct from the
first stacking configuration for any value of $\delta_{z}$. Furthermore, if we
take the second stacking configuration and perform a quarter translation, we
arrive at the first stacking configuration with $\delta_{z}=a_{z}/4$.

As we have discussed, the value of $\Phi(0)$ determines the disclination
charge parity for a system having open boundaries. However, as we discuss in
Appendix A, the absolute disclination charge parity is not a very practical
quantity to consider since it can depend on the choice of inversion center.
For insulators having $T^{z}$ discrete translation symmetry, a more useful
quantity to consider is the difference in the value of the disclination bound
charge $Q_{\text{disc}}\mod(2\frac{\Theta_{F}}{\pi})$ between a system with
open boundaries and a system of the same size with periodic boundaries. For a
system of size $L_{z}=a_{z}N_{z}$, the difference in charge is

$Q_{\text{FA}}=\frac{\Theta_{F}}{2\pi}\left[\frac{\Phi(0)}{\pi}-\mathcal{S}_{\text{3D}}N_{z}\right]\mod(2\frac{\Theta_{F}}{\pi}).$
(25)

Since $\mathcal{S}_{\text{3D}}$ is an even integer for spin-1/2 insulators
with TRS, $Q_{\text{FA}}=0$ or $\frac{\Theta_{F}}{\pi}$ are the distinct
values. When $Q_{\text{FA}}$ is non-zero, we will refer to the system as
having a disclination filling anomaly. The usual filling anomaly[61, 62, 75,
76] reflects an inability to symmetrically deform from periodic to open
boundary conditions while keeping the charge in the system constant. Here, the
disclination filling anomaly reflects an inability to change from periodic to
open boundary conditions while keeping both inversion symmetry and the charge
on disclination lines constant.

In Appendix A we argue that the disclination-line filling anomaly necessarily
vanishes for any layered insulator. Here we are specifically defining a
layered insulator as an insulator where electrons from different unit cells in
the $z$-direction are fully decoupled. Indeed, cutting a layered system to
change from periodic to open boundary conditions is completely a trivial
procedure under this definition. Hence, insulators that exhibit a disclination
filling anomaly are not layered in this way.

In summary, we find that we can use the $R\wedge F$ term to characterize the
3D dWZ response, and the disclination filling anomaly of insulators. The
response phenomena lead to a $\mathbb{Z}_{4}\times\mathbb{Z}_{2}$
classification for insulators having $T^{z}$ discrete translation symmetry and
inversion symmetry. On the one hand, the $\mathbb{Z}_{4}$ index is determined
by calculating the disclination charge per length for a system with periodic
boundary conditions. On the other, the $\mathbb{Z}_{2}$ index is determined by
calculating and comparing the total disinclination charge for a system having
periodic boundary conditions and one having open boundary conditions (with the
same length).

For systems that do not have translation symmetry, the value of $\Phi(0)$
remains quantized, since its quantization requires only inversion symmetry.
This means that for systems with only inversion symmetry, there should be two
distinct insulators characterized by $\Phi(0)=0$ and $\Phi(0)=2\pi$
respectively. For open boundary conditions, the two classes of insulators will
have different disclination charge parities. However, without translation
symmetry, it is not possible to define the disclination filling anomaly. The
disclination charge parity will therefore be significant when considering DSMs
with an incommensurate CDW, which lack any form of translation symmetry.

##  IV Gapping a Dirac Semimetal with a Charge Density Wave

The gapless nodes of the DSM are protected by translation, rotation, and
U$(1)$ charge conservation symmetries. By spontaneously breaking one of these
symmetries it is possible to drive the system into a gapped phase with charge
density wave (CDW), nematic, or superconducting orders, respectively. Here, we
consider the gapped insulating phase that is generated by breaking translation
symmetry using CDW order (see Ref. [31] for a discussion of the instabilities
of a DSM to CDW formation). We will pay special attention to the inversion
symmetry of the DSM, and show that there are two distinct classes of
inversion-symmetric, Dirac-CDW insulators. The differences between these two
classes of insulators will be further discussed in Sec. V.

To consider the effects of a CDW distortion on the low-energy degrees of
freedom, we start with a DSM with two Dirac nodes at $\pm{\bf{k}}_{DN},$ and
expand the fermion creation operators around the two Dirac points, as

$c_{\mathbf{r}}\approx\sum_{\bm{q}}e^{i\mathbf{r}\cdot\bm{q}+ir_{z}\mathcal{K}}c_{\bm{q},R}+e^{i\mathbf{r}\cdot\bm{q}-ir_{z}\mathcal{K}}c_{\bm{q},L},$
(26)

and introduce a new set of Pauli matrices $\bm{\tau}$, where $\tau^{z}=+1$
($-1$) for the fermions in the $R$ ($L$) valley. In terms of the $R$ and $L$
fermions, the $8\times 8$ low-energy Hamiltonian for the DSM is

$\displaystyle\mathcal{H}=$ $\displaystyle
q_{x}a_{x}\Gamma_{1}+q_{y}a_{y}\Gamma_{2}+q_{z}a_{z}\sin(\mathcal{K}a_{x})\Gamma_{3}\tau^{z}.$
(27)

The crystalline symmetries act on the low-energy Hamiltonian as follows. The
C4z rotation symmetry acts as

$R_{4}=\exp\left(-i\frac{\pi}{4}[\sigma^{z}s^{z}-2s^{z}]\right).$ (28)

As discussed previously, there are two possible inversion symmetries, site-
centered and bond-centered. The site-centered inversion symmetry acts on the
low-energy Hamiltonian as

$\begin{split}&\mathcal{I}_{s}=\sigma^{z}\tau^{x}=\Gamma_{3}\tau^{x}.\end{split}$
(29)

The bond-centered inversion symmetry can be similarly written as

$\begin{split}&\mathcal{I}_{b}=\mathcal{I}_{s}e^{iq_{z}a_{z}+i\mathcal{K}a_{z}\tau^{z}}.\end{split}$
(30)

This expression can be simplified if we perform a unitary transformation,
$U=\exp(-i\mathcal{K}a_{z}\tau^{z})$, which leaves the Hamiltonian unaffected,
but reduces the bond-centered inversion symmetry to

$\begin{split}&\mathcal{I}_{b}=\mathcal{I}_{s}e^{iq_{z}a_{z}}.\end{split}$
(31)

This unitary transformation will simplify our later discussions, and so we
will implicitly assume such a unitary transformation has been performed when
discussing continuum models with bond-centered inversion symmetry.

A CDW can be dynamically generated through interactions, as discussed in Ref.
[31]. In the mean field limit, a CDW order corresponds to a translation
symmetry breaking term. We consider two possible types of translation symmetry
breaking terms here. First there is an onsite term that is given in terms of
the original lattice model by

$\hat{H}_{\text{site}}=\sum_{\mathbf{r}}|\Delta_{s}|\cos(2\mathcal{K}r_{z}+\theta_{s})[c^{\dagger}_{\mathbf{r}}\Gamma_{3}c_{\mathbf{r}}].$
(32)

where $|\Delta_{s}|$ is the amplitude of the onsite CDW, and $\theta_{s}$ is
its phase. Second, there is a bond-order term in the $z$-direction that
modulates the hopping terms by

$\hat{H}_{\text{bond}}=\sum_{\mathbf{r}}|\Delta_{b}|\cos(2\mathcal{K}r_{z}+\theta_{b})[c^{\dagger}_{\mathbf{r}}\Gamma_{3}c_{\mathbf{r+{z}}}]+H.c.,$
(33)

where $|\Delta_{b}|$ is the amplitude of the bond-order CDW, and $\theta_{b}$
is its phase. Note that we are explicitly choosing terms that carry momentum
equal to the separation of the Dirac nodes ($2\mathcal{K}$). This nesting
condition allows the mean field term to couple the low-energy degrees of
freedom at the two Dirac nodes. We also could have chosen different
$\Gamma$-matrix structure in Eqs. 32,33, but the choice of $\Gamma_{3}$
ensures these terms can open a mass gap in the continuum Hamiltonian.

Returning to the low-energy theory, Eq. 27, the CDW mean field terms in Eq. 32
and 33 generically induce the following terms in the continuum Hamiltonian

$\begin{split}\mathcal{H}_{\text{mass}}=M(\bm{q})\Gamma_{3}\tau^{x}+M^{\prime}(\bm{q})\Gamma_{3}\tau^{y}.\end{split}$
(34)

The two mass terms are off diagonal in $\tau$ and, hence, couple the two
Dirac-nodes. The spectrum of the continuum Hamiltonian with the $M$ and
$M^{\prime}$ terms is 4-fold degenerate with energy eigenvalues

$\begin{split}E_{\pm}(\bm{q})=\pm\Big{[}&(q_{x}a_{x})^{2}+(q_{y}a_{y})^{2}+\left(q_{z}a_{z}\sin(\mathcal{K}a_{z})\right)^{2}\\\
&+M(\bm{q})^{2}+M^{\prime}(\bm{q})^{2}\Big{]}^{1/2}.\end{split}$ (35)

For our analysis, we are interested in the situation where the CDW is weak,
and tangibly affects the low-energy modes along only the $k_{x}=k_{y}=0$ high
symmetry line. In terms of the original lattice model, this leads to the
requirement that $|\Delta_{s/b}|\ll b_{xy}$ (recall that $b_{xy}$ determines
the gap away from $k_{x}=k_{y}=0$).

Now let us consider the effects of inversion symmetry. By comparing Eq. 31 and
29 to Eq. 34, we find that in presence of inversion $M$ and $M^{\prime}$ must
be an even and odd functions of $\bm{q}$, respectively. Therefore
$M^{\prime}(\bm{q})$ must vanish at the Dirac point, $\bm{q}=0$ in inversion-
symmetric systems. This leads to two distinct inversion-symmetric insulators,
one with and $M^{\prime}(0)=0$ and $M(0)>0,$ and one with $M^{\prime}(0)=0$
and $M(0)<0$. We now wish to relate the values of the $M$ and $M^{\prime}$ to
the microscopic mean-field CDW parameters, $|\Delta_{s}|$, $|\Delta_{b}|$,
$\theta_{s}$ and $\theta_{b}$. Recall that we are implicitly using a unitary
transformation when considering bond-centered inversion symmetry, such that
the inversion operation takes on the simple form in Eq. 31. Because of this,
we will have to consider the site-centered and bond-centered cases (i.e. those
with and without the unitary transformation) separately when relating $M$ and
$M^{\prime}$ to the mean-field parameters. Plugging Eq. 26 into Eq. 32 and 33,
we find that for site-centered inversion symmetry $M$ and $M^{\prime}$ are
related to the mean-field terms as

$\begin{split}M(0)&=|\Delta_{s}|\cos(\theta_{s})+|\Delta_{b}|\cos(\theta_{b}-\mathcal{K}a_{z}),\\\
M^{\prime}(0)&=|\Delta_{s}|\sin(\theta_{s})+|\Delta_{b}|\sin(\theta_{b}-\mathcal{K}a_{z}).\\\
\end{split}$ (36)

To be compatible with inversion symmetry we need $M^{\prime}(0)$ to vanish,
hence, the onsite term is compatible with site-centered inversion symmetry
when $\theta_{s}=0$, or $\pi,$ and the hopping term is compatible with site-
centered inversion symmetry when $\theta_{b}=\mathcal{K}a_{z}$, or
$\mathcal{K}a_{z}+\pi$. Similarly, for bond-centered inversion symmetry (using
the unitary transformation discussed above) $M$ and $M^{\prime}$ are related
to the mean-field terms as

$\begin{split}M(0)&=|\Delta_{s}|\cos(\theta_{s}+\mathcal{K}a_{z})+|\Delta_{b}|\cos(\theta_{b}),\\\
M^{\prime}(0)&=|\Delta_{s}|\sin(\theta_{s}+\mathcal{K}a_{z})+|\Delta_{b}|\sin(\theta_{b}).\\\
\end{split}$ (37)

Hence, the onsite term is compatible with bond-centered inversion symmetry
when $\theta_{s}=-\mathcal{K}a_{z}$, or $-\mathcal{K}a_{z}+\pi,$ and the
hopping term is compatible with bond-centered inversion when
$\theta_{b}=0,\pi$.

We note that the commensurate CDW with period 2 ($\mathcal{K}=\pi/2a_{z}$) is
a special case, since the corresponding onsite term is compatible with site
centered inversion symmetry for all values of $\theta_{s}$, while the hopping
term is not compatible with site-centered inversion symmetry for any value of
$\theta_{b}$. Similarly, the hopping term is compatible with bond-centered
inversion symmetry for all $\theta_{b}$, while the onsite term is never
compatible with bond-centered inversion symmetry. This can be confirmed by
direct inspection of Eq. 32 and 33.

The two inversion-symmetric Dirac-CDW insulators ($M(0)\lessgtr 0$,
$M^{\prime}(0)=0$) are anisotropic topological crystalline insulators (TCI).
As we shall show in the next section, these have different $\mathbb{Z}_{2}$
indices using the classification scheme discussed in Sec. III.2.2. Although
they have different crystalline-electromagnetic responses, both classes of
insulators share many important physical properties. Namely, both classes of
insulators have symmetry protected gapless surface states on boundaries normal
to the $x$ and $y$ directions. This can be understood by noting that the
Dirac-CDW insulator is adiabatically connected to two TRS-related copies of
the Weyl-CDW insulator. Since the Weyl-CDW insulators have gapless 1D chiral
surface modes on the $x$ and $y$ surfaces that circulate around the $z$-axis
[21], the Dirac-CDW insulator has helical surface modes that counter-propagate
around the $z$-axis [31]. These edge modes match those that are found in a
stack of 2D quantum spin Hall (QSH) layers. The helical surface modes also
indicate that edge and screw dislocations of the CDW bind helical modes [31].
The helical modes of the Dirac-CDW insulator are interesting in their own
right, but they do not have an immediate impact on our forthcoming discussion
of the disclination responses of the Dirac-CDW insulator, provided that
disclinations are far away from any boundaries normal to the $x$ or $y$
directions. We note that while the surfaces normal to the $x$ and $y$
directions are gapless due to time-reversal symmetry, surfaces normal to the
$z$ direction can be gapped while preserving all relevant symmetries.

##  V crystalline-electromagnetic responses of the Dirac-CDW insulator

We now turn our attention to the topological responses of Dirac-CDW
insulators. The usual topological electromagnetic responses, namely the 3D
Hall response and axion electrodynamics, are trivial for a Dirac-CDW insulator
since such an insulator is essentially two time-reversal symmetry related,
Weyl-CDW insulators. However, even though Dirac-CDW insulators have trivial
electromagnetic responses, they can host non-trivial crystalline-
electromagnetic responses. The exact nature of these responses depends on if
the CDW preserves a subgroup of the $T^{z}$ discrete translation symmetry of
the lattice (i.e., if the CDW is commensurate with the lattice) and if the CDW
preserves inversion symmetry.

Below we will show the following results. If the CDW is commensurate, the
Dirac-CDW insulators have a quantized 3D dWZ response. If the CDW preserves
inversion symmetry but is incommensurate, we find there are two distinct
inversion-symmetric Dirac-CDW insulators that differ by a $\Phi=2\pi$ $R\wedge
F$ term. Based on our previous arguments, two insulators differing by an
$R\wedge F$ term will have different disclination charge parities. If the CDW
is both commensurate and preserves inversion symmetry, we find that there are
two distinct inversion-symmetric Dirac-CDW insulators that have the same 3D
dWZ response, but differ by a $\Phi=2\pi$ $R\wedge F$ term. One of these
insulators has a disclination filling anomaly, while the other does not. We
will establish these results analytically, and confirm them through explicit
numerical lattice model calculations.

###  V.1 The 3D discrete Wen-Zee response of the Dirac-CDW insulator

In this subsection, we will show that the Dirac-CDW insulator realizes the 3D
dWZ response when the CDW is commensurate. Such a response cannot occur for an
incommensurate CDW, as a well-defined 3D dWZ response requires discrete
translation symmetry.

To this end, let us consider a Dirac-CDW insulator generated by a commensurate
CDW with a period $\frac{p}{q}a_{z}$, where $a_{z}$ is lattice constant along
the $z$-direction ($\mathcal{K}=\frac{q}{p}\frac{\pi}{a_{z}}$ in Eq. 1). Such
a system has an enlarged unit cell of length $\tilde{a}_{z}=pa_{z}$, and a
reduced $T^{z}$ discrete translation symmetry $z\rightarrow z+pa_{z}$. The 3D
dWZ response of such a system can be written analogously to Eq. 14 as

$\begin{split}\mathcal{L}_{\text{3D-WZ
CDW}}&=\mathcal{S}_{\text{3D}}\frac{\tilde{G}_{z}}{4\pi^{2}}\epsilon^{ijk}\omega_{i}\partial_{j}A_{k},\\\
\tilde{G}_{z}&=\frac{2\pi}{pa_{z}},\\\
\mathcal{S}_{\text{3D}}&=2q.\end{split}$ (38)

Here, $\tilde{G}_{z}$ is the reciprocal lattice vector along the $z$-direction
with respect to the enlarged lattice. This response indicates that a
disclination line of a system of size $L_{z}=a_{z}N_{z}$ binds charge

$\begin{split}Q_{\text{disc}}&=\mathcal{S}_{\text{3D}}\frac{\Theta_{F}}{2\pi}N_{z}/p\mod(2N_{z}/p),\\\
&=2\frac{q}{p}\frac{\Theta_{F}}{2\pi}N_{z}\mod(2N_{z}/p),\end{split}$ (39)

where the $\mod(2N_{z}/p)$ ambiguity comes from adding a local Kramers
degenerate pair of electrons to each enlarged unit cell. Note that we assume
$N_{z}$ is a multiple of $p$ in order for the system to be invariant under the
reduced translation symmetry.

We can arrive at Eq. 38 by starting with the anomaly equation of the DSM, Eq.
7, and enlarging the unit cell,

$\begin{split}&a_{z}\rightarrow\tilde{a}_{z}=pa_{z},\\\
&G_{z}\rightarrow\tilde{G}_{z}=G_{z}/p\end{split}$ (40)

Using that $\nu=q/p$, we then arrive at Eq. 38. Note that after this
redefinition of the unit cell, Eq. 7 has an integer coefficient and is
therefore gauge invariant and anomaly free. The system no longer needs to have
gapless degrees of freedom to compensate. Indeed, the Dirac nodes can be
gapped out after enlarging the unit cell, as discussed previously. The 3D dWZ
response of the Dirac-CDW insulator is therefore the descendant of the
anomalous DSM response, when the anomaly is removed by breaking translation
symmetry.

We can also argue that the commensurate Dirac-CDW insulator realizes the
response in Eq. 38 by treating the DSM (without a CDW) as a family of 2D
insulators that are parameterized by $k_{z}$. We recall from above that in
this interpretation, the 2D insulators with $|k_{z}|<\mathcal{K}$ are QSH
insulators with $\mathcal{S}_{\text{2D}}=2$, the insulators with
$|k_{z}|>\mathcal{K}$ are trivial ($\mathcal{S}_{\text{2D}}=0$), and the
system is gapless at $|k_{z}|=\mathcal{K}$. The total charge on a disclination
in the DSM is then simply the sum of contributions coming from the QSH
insulators, and the gapless degrees of freedom at $k_{z}=\pm\mathcal{K}$. To
get an insulator we can turn on the CDW which will couple the 2D system
indexed by $k_{z}$ to the 2D system indexed by $k_{z}+2\mathcal{K}$. If the
CDW is small but finite, the coupling will affect only the gapless systems at
$k_{z}=\pm\mathcal{K}$. The effect of the CDW on the gapped states that are
initially indexed $|k_{z}|<\mathcal{K}$ and $|k_{z}|>\mathcal{K}$ will be
negligible and will therefore not affect the topological responses of these
states. In particular, at weak coupling the topological responses of the
$|k_{z}|<\mathcal{K}$ (i.e., the QSH insulators) will not be affected.
Therefore, if we add a disclination to the system, each of the
$|k_{z}|<\mathcal{K}$ insulators will still bind charge
$\frac{\Theta_{F}}{\pi}$. Additionally, a direct calculation shows that the
two hybridized states at $k_{z}=\pm\mathcal{K}$, bind a net charge of
$\frac{\Theta_{F}}{\pi}$. The states for $|k_{z}|>\mathcal{K}$ correspond to
trivial 2D insulators and therefore do not bind any charge. For a system of
size $L_{z}=a_{z}N_{z},$ where $k_{z}$ is quantized as a multiple of
$2\pi/N_{z}$ (assuming periodic boundary conditions), there will be a total of
$\frac{q}{p}N_{z}-1$ states with $|k_{z}|<\mathcal{K}$ (recall that $N_{z}$
must be a multiple of $q$). A disclination therefore binds a total charge of
$2\frac{q}{p}\frac{\Theta_{F}}{2\pi}N_{z}$, in agreement with Eq. 39.

###  V.2 The disclination filling anomaly of the Dirac semimetal-charge
density wave insulator

In Sec. IV we showed that for a given $\mathcal{K}$, there are two distinct
inversion-symmetric Dirac-CDW insulators. In this section, we will show that
the effective response theories for these insulators differ by an $R\wedge F$
term with $\Phi=2\pi$. If the CDW is commensurate, such that the system has
$T^{z}$ discrete translation symmetry, one of these insulators will have a
disclination filling anomaly, and the other will not (see Sec. III.2). More
generally, i.e., even in incommensurate CDWs, the difference in the $R\wedge
F$ term indicates that the disclination charges of the two insulators will
differ by $\Theta_{F}/\pi\mod(2\Theta_{F}/\pi)$ when open boundaries are
present.

To show that the two inversion-symmetric insulators differ by a $\Phi=2\pi$
$R\wedge F$ term, we consider the low-energy Lagrangian for the DSM subject to
a mean-field CDW term (see Sec. IV and Eq. 34), and minimally couple the Dirac
fermions to the U$(1)$ charge and C4z rotation gauge fields,

$\begin{split}\mathcal{L}&=\bar{\Psi}[i\bar{\Gamma}^{\mu}D_{\mu}+M+M^{\prime}\bar{\Gamma}^{5}]\Psi,\\\
D_{\mu}&=\partial_{\mu}-iA_{\mu}-i\omega_{\mu}\left[\frac{1}{2}\sigma^{z}s^{z}-s^{z}\right],\end{split}$
(41)

where $\Psi$ is an 8 component spinor,
$\bar{\Psi}=\Psi^{\dagger}\bar{\Gamma}^{0}$, and

$\begin{split}&\bar{\Gamma}^{x}=\sigma^{y}s^{z}\tau^{x},\phantom{=}\bar{\Gamma}^{y}=\sigma^{x}\tau^{x},\phantom{=}\bar{\Gamma}^{z}=\tau^{y},\\\
&\bar{\Gamma}^{0}=\sigma^{z}\tau^{x},\phantom{=}\bar{\Gamma}^{5}=\tau^{z}.\end{split}$
(42)

The covariant derivative $D_{\mu}$ minimally couples the low-energy fermions
to the U$(1)$ gauge field, $A_{\mu}$, and the C4z gauge field, $\omega_{\mu}$.
The latter couples via the C4z angular momentum operator,
$\frac{1}{2}\sigma^{z}s^{z}-s^{z}$ [60]. In this minimal coupling procedure,
we are treating $\omega$ as a continuous gauge field with fluxes that are
quantized in multiples of $\pi/2$ as discussed before. Compared to Eq. 34, we
are suppressing any momentum dependence of the CDW mass terms $M$ and
$M^{\prime}$, as it is either irrelevant at this order, or can be absorbed
into a redefinition of the Fermi velocity.

We are interested in the difference between the topological response of the
two gapped inversion-symmetric phases, i.e., inversion symmetry sets
$M^{\prime}=0$, and we want to compare the insulators having $M>0$ or $M<0$.
To this end, we set $M=\bar{M}\cos(\theta)$ and
$M^{\prime}=\bar{M}\sin(\theta)$, and consider a process where $\theta$ is
smoothly increased from $0$ to $\pi$. The fermions remain gapped during this
process (Eq. 35), and the effective response theory (treated as a function of
$\theta$) is found by integrating out the massive fermions. The resulting
effective response theory contains the topological term,

$\displaystyle\mathcal{L}_{\mathrm{top}}[\theta]=\frac{\theta}{2\pi^{2}}\epsilon^{\mu\nu\lambda\eta}\partial_{\mu}\omega_{\nu}\partial_{\lambda}A_{\eta}\quad,$
(43)

which arises from a triangle diagram with external legs $A_{\mu}$,
$\omega_{\mu}$, and $\theta$. The $\theta=0$ and $\theta=\pi$ inversion-
symmetric insulators therefore differ by a $\Phi=2\pi$ $R\wedge F$ term (Eq.
16).

Figure 6:  Layered construction of the period-2 Dirac-CDW insulators, where a
two adjacent $z$-layers (teal) pair up to form a QSH insulator (dashed box)
with $\mathcal{S}_{\text{3D}}=2$. For a system of length $L_{z}=a_{z}N_{z}$
with open boundaries, this stacking can be done such that there are $N_{z}/2$
QSH insulators (a) or $N_{z}/2-1$ QSH insulator layers (b). This describes Eq.
49 with $|\Delta_{b}|=1$, and $\theta_{b}=0$ and $\pi$ respectively. When
there $N_{z}/2-1$ QSH insulator layers, there are a pair of decoupled
$z$-layers on the top and bottom surfaces (orange).

Based on this, and our discussions in Sec. III.2 and V.1, we therefore
conclude that the two distinct inversion-symmetric Dirac-CDW insulators differ
by a $\Phi=2\pi$ $R\wedge F$ term. For commensurate CDWs, this difference is
related to the presence or absence of a disclination filling anomaly. For
incommensurate CDWs, the difference in $R\wedge F$ term indicates that the two
insulators have different disclination charge parities.

###  V.3 Analytic solution for the period 2 charge density wave

To demonstrate the Dirac-CDW insulator responses we will first consider the
insulators that form when a DSM with $\mathcal{K}=\pi/2a_{z}$ is driven into a
massive phase by a period-2 CDW with inversion symmetry. As noted in Sec. II,
there are two types of inversion symmetries to consider, the site-centered
inversion symmetry that sends $z\rightarrow-z,$ and the bond-centered
inversion symmetry that sends $z\rightarrow-z+a_{z}$. For a system with open
boundaries, the site-centered inversion symmetry requires an odd number of
sites, while the bond-centered inversion symmetry requires an even number of
sites. For an insulator with a period-2 CDW, it is natural to require an even
number of sites for commensurability, and so we will consider a lattice model
with bond-centered inversion symmetry.

The mean field limit of the DSM (Eq. 1) having a period 2 CDW hopping term
(Eq. 33) is

$\begin{split}&H_{\text{DSM-CDW}}=\begin{bmatrix}H^{xy}&H^{z}\\\
H^{z\dagger}&H^{xy}\end{bmatrix}\\\
&H^{xy}=\sin(k_{x}a_{x})\Gamma_{1}+\sin(k_{y}a_{y})\Gamma_{2}\\\
&\phantom{====}-b_{xy}(2-\cos(k_{x}a_{x})-\cos(k_{y}a_{y}))\Gamma_{3}\\\
&H^{z}=(1+|\Delta_{b}|\cos(\theta_{b}))\Gamma_{3}\\\
&\phantom{===}+e^{i\tilde{k}_{z}\tilde{a}_{z}}(1-|\Delta_{b}|\cos(\theta_{b}))\Gamma_{3}\end{split}$
(44)

where $\tilde{k}_{z}=k_{z}/2$ is the momentum of the folded Brillouin zone
with conjugate position, $\tilde{z}=2z.$ The unit cells indexed by $\tilde{z}$
are twice as large as in the gapless DSM, and have length
$\tilde{a}_{z}=2a_{z}$. We will consider a weak CDW where $|\Delta_{b}|\ll
b_{xy}$. The bond-centered inversion symmetry acts on $H_{\text{DSM-CDW}}$ as

$\displaystyle\mathcal{I}_{b}=\sigma^{z}t^{x}$ (45)

where the $\bm{t}$ Pauli matrices act on the $4\times 4$ blocks that make up
$H_{\text{DSM-CDW}}$. The even-$z$ (odd-$z$) sublattices of Eq. 1 correspond
to $t^{z}=+1(-1)$ in Eq. 44. The $\bm{t}$ Pauli matrices should not to be
confused with the $\bm{\tau}$ Pauli matrices that differentiate the two Dirac
nodes in the low-energy Hamiltonian, Eq. 27.

There are two inversion-symmetric gapped phases to consider,
$\cos(\theta_{b})>0$ and $\cos(\theta_{b})<0$. For simplicity, we will
restrict our attention to the cases where $\theta_{b}=0$ or $\theta_{b}=\pi$.
For the period-2 CDW, other values of $\theta_{b}$ are equivalent to one of
these two values by a redefinition of $|\Delta_{b}|$. To demonstrate that the
$\theta_{b}=0$ and $\theta_{b}=\pi$ insulators differ by a $\Phi=2\pi$
$R\wedge F$ term we can add an inversion symmetry breaking onsite term

$H^{\prime}=2\Delta^{\prime}\Gamma_{3}t^{z}$ (46)

which is simply the onsite CDW term, Eq. 32 with
$\Delta^{\prime}=|\Delta_{s}|\cos(\theta_{s})$. Indeed, in the continuum
limit, $H_{\text{DSM-CDW}}+H^{\prime}$ is equivalent to Eq. 41 with
$M\propto|\Delta_{b}|\cos(\theta_{b})$ and $M^{\prime}\propto\Delta^{\prime}$.
Hence, based on Sec. V.2 the $\theta_{b}=0$ and $\theta_{b}=\pi$ phases differ
by a $\Phi=2\pi$ $R\wedge F$ term.

To analyze the full lattice model, it is useful to consider $|\Delta_{b}|=1$
and $\theta_{b}=0$ where the Hamiltonian is independent of $\tilde{k}_{z}$,
and can be treated as a stack of 2D insulating layers. The $t^{x}=\pm 1$
sectors of each layer are equivalent to QSH Hamiltonians, Eq. 12, with masses
$m=2\mp\frac{2}{b_{xy}}$ respectively. Recalling that we have set
$b_{xy}\gg|\Delta_{b}|=1$ , the $t^{x}=-1$ sector is in the QSH phase, while
the $t^{x}=+1$ sector is in the trivial phase. We can therefore conclude that
for $|\Delta_{b}|=1$ and $\theta_{b}=0$ each CDW period hosts a single 2D QSH
insulator, in agreement with our earlier arguments.

To determine the 3D dWZ response in this layered limit, let us consider adding
a disclination to a lattice of length $L_{z}=a_{z}N_{z}$ having periodic
boundary conditions. In order for the system to preserve the reduced
translation symmetry, $z\rightarrow z+2a_{z}$, $N_{z}$ must be an even integer
such that there are $N_{z}/2$ CDW periods. Since each CDW period contains a
single QSH insulator which has a $\mathcal{S}_{\text{2D}}=2$ 2D dWZ response,
the charge bound to a disclination is

$Q_{\text{disc}}=\frac{1}{2}\frac{\Theta_{F}}{\pi}N_{z}\mod(N_{z}).$ (47)

Comparing this to Eq. 39, we find that the full 3D system has a 3D dWZ
response with $\mathcal{S}_{\text{3D}}=2$, as expected for a layered system.

Next let us consider the $|\Delta_{b}|=1$ and $\theta_{b}=\pi$ insulator. This
insulator is related to the $\theta_{b}=0$ insulating phase by a translation
in the $z$-direction by half of the doubled unit cell (i.e., translation by
one unit cell of the original DSM lattice model)

$T_{z-1/2}=\begin{bmatrix}0&e^{ik_{z}}\\\ 1&0\end{bmatrix}.$ (48)

The $|\Delta_{b}|=1$ and $\theta_{b}=\pi$ insulator can therefore also be
treated as a stack of 2D insulating layers. The 2D insulators equally occupy
the $t^{z}=+1$ sublattice of layer $\tilde{z}$ and the $t^{z}=-1$ sublattice
of layer $\tilde{z}+1$. Based on our previous results, the bulk of the
$|\Delta_{b}|=1$ and $\theta_{b}=\pi$ insulator also hosts a single 2D QSH
insulator per CDW period, and therefore has $\mathcal{S}_{\text{3D}}=2$.
Hence, both the $|\Delta_{b}|=1$ and $\theta_{b}=0,\pi$ insulators have an
identical 3D dWZ response with $\mathcal{S}_{\text{3D}}=2$. Since the 3D dWZ
response is quantized for insulators with rotation symmetry and charge
conservation, we conclude that all insulating phases of our model with
$b_{xy}\gg|\Delta_{b}|\neq 0$ have this response.

While these phases are not distinguished by their dWZ response, based our
previous discussions in Sec. III.2, we expect two distinct classes of
inversion-symmetric insulators can be distinguished by another response. In
fact, we have shown that the effective response theories for these insulators
differ by a $\Phi=2\pi$ $R\wedge F$ term. Based on our results from Sec. V.2,
the $\theta_{b}=0$, and $\theta_{b}=\pi$ insulators should fall into different
classes, and one of these insulators should have a disclination filling
anomaly, while the other should not. As we shall show, the $\theta_{b}=\pi$
insulator has a disclination filling anomaly, while the $\theta_{b}=0$
insulator does not.

To show this, let us consider a system described by Eq. 44 of length
$L_{z}=aN_{z}$ ($N_{z}$ even) having open boundaries in the $z$-direction.
Using position space in the $z$-direction, and momentum space in the $x$ and
$y$-directions, the Hamiltonian is

$\begin{split}\hat{H}_{\text{DSM-
CDW}}&=\sum^{N_{z}/2}_{\tilde{z}=1}[c^{\dagger}_{k_{x},k_{y},\tilde{z}}H^{xy}c_{k_{x},k_{y},\tilde{z}}\\\
&+(1+|\Delta_{b}|\cos(\theta_{b}))c^{\dagger}_{k_{x},k_{y},\tilde{z}}\Gamma_{3}t^{x}c_{k_{x},k_{y},\tilde{z}}\\\
&+(1-|\Delta_{b}|\cos(\theta_{b}))c^{\dagger}_{k_{x},k_{y},\tilde{z}}\Gamma_{3}t^{x}c_{k_{x},k_{y},\tilde{z}+1}],\end{split}$
(49)

where $\tilde{z}=2z$ labels the doubled unit cells along the $z$-direction. As
before, we consider the representative insulators with $|\Delta_{b}|=1$ and
$\theta_{b}=0$ or $\pi$.

For $|\Delta_{b}|=1$ and $\theta_{b}=0$, the spectrum is fully gapped, and
different $\tilde{z}$-layers are fully decoupled. Furthermore, each of these
layers has $\mathcal{S}_{\text{2D}}=2,$ and constitutes a QSH insulator, i.e.
one non-trivial QSH per $\tilde{z}$. This stacking configuration is shown in
Fig. 6(a). If we add a $\pi/2$ disclination-line to the insulator, each
$\tilde{z}$-layer will bind charge $\Theta_{F}/\pi$. For a system of length
$L_{z}=a_{z}N_{z}$, the total charge on the disclination line is the sum of
the charges on the 2D layers,

$\begin{split}Q_{\text{disc}}&=\frac{1}{2}\frac{\Theta_{F}}{\pi}N_{z}\mod(2\frac{\Theta_{F}}{\pi}),\end{split}$
(50)

where the $2\Theta_{F}/\pi$ ambiguity comes from the ability to add 2D
insulators having $\mathcal{S}_{\text{2D}}\in 2\mathbb{Z}$ to the top and
bottom surfaces while preserving inversion symmetry (recall that
$\mathcal{S}_{\text{2D}}$ is an even integer for spin-1/2 insulators with
TRS). Comparing Eq. 50 and 47, we find that the disclination charge of the
$|\Delta_{b}|=1$ and $\theta_{b}=0$ insulator changes by
$0\mod(2\Theta_{F}/\pi)$ when changing boundary conditions, indicating that
this system does not have a disclination filling anomaly.

We can similarly find the charge on the disclination line when
$|\Delta_{b}|=1$ and $\theta_{b}=\pi$. For open boundaries, the $t^{z}=-1$
sublattice of the $\tilde{z}=1$ layer is fully decoupled. The Hamiltonian for
the fermions in this sublattice/layer sector is just $H^{xy}$ from Eq. 44,
which has a gapless point at $k_{x}=k_{y}=0$. The fermions on the $t^{z}=+1$
sublattice of the $\tilde{z}=N_{z}/2$ layer at the top are also decoupled from
the other layers/sublattices and are gapless. These sectors can be gapped out,
with inversion symmetry preserving mass terms $\propto\Gamma^{3}$. When such a
perturbation has been applied we can determine the charge distribution on a
$\pi/2$-disclination line. As noted previously, the bulk of the
$|\Delta_{b}|=1$, $\theta_{b}=\pi$ insulator is composed of layers of QSH
insulators having $\mathcal{S}_{\text{2D}}=2$. These insulators have equal
weight on the $t^{z}=+1$ sublattice of $\tilde{z}$ and the $t^{z}=-1$
sublattice of $\tilde{z}+1$. Since these QSH insulators live between
$\tilde{z}$ layers, a system of length $L_{z}=aN_{z}$ with open boundaries in
the z-direction will have $N_{z}/2-1$ QSH layers. This stacking configuration
is shown in Fig. 6(b). To compute the disclination charge we also note that
the gapping perturbation on the surfaces can generate boundary insulators that
also have 2D dWZ responses with $\mathcal{S}_{\text{2D-bnd}}\in 2\mathbb{Z}$.
However, the discrete Wen-Zee shift of the boundary insulators at the top and
bottom must be the same because of inversion symmetry. Counting up the
contributions from the bulk QSH layers and the boundary insulating layers, the
total charge on the disclination line is

$\begin{split}Q_{\text{disc}}=\frac{1}{2}\frac{\Theta_{F}}{\pi}(N_{z}-2+2\mathcal{S}_{\text{2D-bnd}})\mod(2\frac{\Theta_{F}}{\pi}).\end{split}$
(51)

Using Eqs. 51 and 47, we find that the charge on the disclination lines of the
$|\Delta_{b}|=1$, $\theta_{b}=\pi$ insulator changes by
$\Theta_{F}/\pi\mod(2\Theta_{F}/\pi)$ when changing boundary conditions. The
$\theta_{b}=\pi$ insulator therefore has a disclination filling anomaly.
Physically, the difference in charge can be traced back to the fact that if we
place the QSH insulators between $\tilde{z}$ layers (as we did for the
$\theta_{b}=\pi$ insulator), there will be one less QSH insulator than if we
place the QSH insulators on the $\tilde{z}$ layers (as we did for the
$\theta_{b}=0$ insulator). Boundary effects can, effectively, add a QSH
insulator to the top and bottom layers. But even with the effect, the parity
of the number of QSH insulators is fixed in each stacking configuration. This
leads to the quantized difference in disclination charge found above.

In summary, we find two distinct inversion-symmetric insulators can arise from
the coupling a $\mathcal{K}=\pi/2a_{z}$ DSM to a CDW. As expected, both
insulators have a $\mathcal{S}_{\text{3D}}=2$ 3D dWZ response, and one has a
disclination filling anomaly while the other does not.

###  V.4 Numeric analysis

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: (a, b) The charge bound to a disclination as a function of $z$ for $\mathcal{K}=\pi/2a_{z}$, $\Delta_{s}=0$, $\Delta_{b}=0.5$, $n_{x}=15$, and $n_{z}=14$ with (a) open and (d) periodic boundary conditions. (b, e) The charge bound to a disclination as a function of $z$ for $\mathcal{K}=\pi/3a_{z}$, $\Delta_{s}=0.5$, $\Delta_{b}=0$, $n_{x}=15$, and $n_{z}=15$ with (b) open and (e) periodic boundary conditions. (c) The charge bound to a disclination as a function of $z$ for $\mathcal{K}=\sqrt{2}\pi/4$, $\Delta_{s}=0.5$, $\Delta_{b}=0$, $n_{x}=15$, and $n_{z}=15$ with open boundary conditions. In all cases we set $b_{xy}=1.0$ and $\Delta_{\text{surf}}=0.25$. The solid black and dashed red lines indicate $\theta_{b/s}=0$ and $\theta_{b/s}=\pi$, respectively. The alternating gray and white shading indicates unit cells of the CDWs. (f) The charge density around a disclination summed over the $z$-direction with $\mathcal{K}=\pi/2a_{z}$, $\theta_{b}=\pi$, $\Delta_{s}=0$, $\Delta_{b}=0.5$, $b_{xy}=1.0$, $\Delta_{\text{surf}}=0.25$, $n_{x}=15$, and $n_{z}=12$. The charge distribution is qualitatively identical for all other cases considered in this work. $\mathcal{K}a_{z}$ | $|\Delta_{s}|$ | $|\Delta_{b}|$ | $\theta_{s/b}$ | $N_{z}$ |  $Q_{\text{disc}}$ OBC |  $Q_{\text{disc}}$ PBC |  DL Filling Anomaly
---|---|---|---|---|---|---|---
$\pi/2$ | 0 | $1/2$ | 0 | 14 | 7/2 | 7/2 | No
$\pi/2$ | 0 | $1/2$ | $\pi$ | 14 | 4 | 7/2 | Yes
$\pi/3$ | $1/2$ | 0 | 0 | 15 | 5/2 | 5/2 | No
$\pi/3$ | $1/2$ | 0 | $\pi$ | 15 | 3 | 5/2 | Yes
$\pi/2\sqrt{2}$ | $1/2$ | 0 | 0 | 15 | 5/2 | N/A | N/A
$\pi/2\sqrt{2}$ | $1/2$ | 0 | $\pi$ | 15 | 3 | N/A | N/A
Table 1: Summary of the numerical analysis of the crystalline-electromagnetic
responses of the Dirac-CDW insulator. The Hamiltonian for the DSM is given by
Eq. 1 with Dirac nodes at $\bm{k}=(0,0,\pm\mathcal{K})$ and $b_{xy}=1$. The
CDW is included at the via either an inversion-symmetric onsite potential
($\hat{H}_{s}$ in Eq. 32 with $\theta_{s}=0$ or $\pi$) or an inversion-
symmetric hopping amplitude ($\hat{H}_{b}$ in Eq. 33 with $\theta_{b}=0$ or
$\pi$) with period $\pi/\mathcal{K}$. For the commensurate values of
$\mathcal{K}=\pi/2$ and $\pi/3$) the total charge bound to a disclination is
calculated for periodic and open and boundary conditions, the difference of
which gives the disclination filling anomaly, defined in Sec. III.2.3.
However, periodic boundary conditions are not compatible with the
incommensurate value of $\mathcal{K}=\pi/2\sqrt{2}$, precluding the
calculation of the disclination filling anomaly.

In this subsection, we numerically calculate the crystalline-electromagnetic
responses for the inversion-symmetric Dirac-CDW insulators that form when a
DSM is coupled to: a commensurate $\mathcal{K}=\pi/2a_{z}$ (period 2) CDW, a
commensurate $\mathcal{K}=\pi/3a_{z}$ (period 3) CDW, and an incommensurate
$\mathcal{K}=\pi/2\sqrt{2}a_{z}$ CDW. These calculations are performed using
the DSM Hamiltonian in Eq. 1 with $b_{xy}=1$. The $\mathcal{K}=\pi/2a_{z}$ CDW
is implemented via the hopping mean field term in Eq. 33 with
$|\Delta_{b}|=1/2$ and either $\theta_{b}=0$ or $\pi$, both of which preserve
bond-centered inversion symmetry. The $\mathcal{K}=\pi/3a_{z}$ CDW is
implemented via the onsite mean field term in Eq. 32 with $|\Delta_{s}|=1/2$
and either $\theta_{s}=0$ or $\pi$, preserving the site-centered inversion
symmetry. The $\mathcal{K}=\pi/2\sqrt{2}a_{z}$ CDW is similarly implemented
with the onsite mean field term in Eq. 32 with $|\Delta_{s}|=1/2$ and either
$\theta_{s}=0$ or $\pi$, also preserving site-centered inversion symmetry.

For each system, we evaluate the Hamiltonian on a lattice having open boundary
conditions and a $\Theta_{F}=\pi/2$ disclination line. The disclination line
is located at the center of the lattice and stretches straight along the
$z$-direction. To reduce the impact of finite size effects, we also add an
inversion-preserving perturbation of the form
$H^{\prime}_{\text{surf}}=\Delta_{\text{surf}}\Gamma_{3}$
($\Delta_{\text{surf}}=1/4$ on both the top and bottom layers) to increase the
gap on the $z$-normal surfaces. Then we calculate the disclination charge by
summing the charge density over all sites within some radial distance from the
disclination core. To further reduce the impact of finite size effects, we
perform these calculations over a range of system sizes and extrapolate to the
infinite-size limit. The details of how the integration radius is determined
and how the extrapolations are performed, are provided in Appendix B.

As expected from our analysis in previous sections, we find that charge is
bound to the disclination line for all Dirac-CDW insulators, as shown in Fig.
7f. Specifically, in Fig. 7 we plot the layer-resolved disclination charge for
each insulator at the largest accessible system size. For the commensurate
CDWs, $\mathcal{K}=\pi/2a_{z}$ and $\mathcal{K}=\pi/3a_{z}$, we find that the
charge distribution oscillates with a period matching the CDW period, with
deviations at open boundaries. There is no notable periodicity for the
$\mathcal{K}=\pi/2\sqrt{2}a_{z}$ CDW, as it is incommensurate with the
lattice. The total disclination charges for each of these insulators
extrapolated to the infinite-size limit are enumerated in Table 1.

For $\mathcal{K}=\pi/2a_{z}$, the total charge bound to the disclination with
open boundary conditions differs by $1/2$ between the $\theta_{b}=0$ and
$\theta_{b}=\pi$ Dirac-CDW insulators. This difference in disclination charge
is a direct manifestation of the difference in the coefficient of the $R\wedge
F$ term between the two Dirac-CDW insulators. Similarly, the disclination
charge differs by $1/2$ between the $\theta_{s}=0$ and $\theta_{s}=\pi$ Dirac-
CDW insulators with $\mathcal{K}=\pi/3a_{z},$ and
$\mathcal{K}=\pi/2\sqrt{2}a_{z}$. We therefore find that the two distinct
inversion-symmetric Dirac-CDW insulators for each CDW wavevector have
different disclination charge parities, matching our analytic predictions.

For the commensurate $\mathcal{K}=\pi/2a_{z}$ and $\mathcal{K}=\pi/3a_{z}$
CDWs, it is also possible to have periodic boundary conditions in the
$z$-direction. Hence in these cases we can determine the 3D dWZ response of
the Dirac-CDW insulators and whether or not the insulators host a disclination
filling anomaly. We find that for periodic boundaries, both the $\theta_{b}=0$
and $\pi$ insulators with $\mathcal{K}=\pi/2a_{z}$ bind charge $1/2$ per CDW
period ($1/4$ per original lattice layer). The $\theta_{b}=0$ and $\pi$
insulators with $\mathcal{K}=\pi/3a_{z}$ also bind charge $1/2$ per CDW period
($1/6$ per original lattice layer). This is exactly the 3D dWZ response that
we have previously discussed. Furthermore, we find that the disclination
charges differ by $1/2$ with periodic and open boundary conditions for the
$\mathcal{K}=\pi/2a_{z}$, $\theta_{b}=\pi$ and $\mathcal{K}=\pi/3a_{z}$,
$\theta_{s}=\pi$ insulators, indicating disclination filling anomalies. The
presence of the disclination filling anomaly can also be observed in the
charge density profiles plotted in Fig. 7, in which it manifests as a
deviation at the boundaries of the disclination charge per CDW period. The
disclination charges for open and periodic boundary conditions are identical
for the $\mathcal{K}=\pi/2a_{z}$, $\theta_{b}=0$ and $\mathcal{K}=\pi/3a_{z}$,
$\theta_{s}=0$ insulators, marking the absence of disclination filling
anomalies. These results agree with the analytic calculations presented in
Sec. V.3. As discussed in Sec. III.2.3, the disclination filling anomaly is
ill-defined for insulators with incommensurate $\mathcal{K}$, like the
$\mathcal{K}=\pi/2\sqrt{2}a_{z}$ Dirac-CDW insulators, since they are
incompatible with periodic boundaries.

##  VI Conclusion and outlook

In this work, we have considered the topological responses of the insulating
state that arises from coupling a 3D DSM to a CDW. Unlike the related Weyl-CDW
insulators, the Dirac-CDW insulators do not display a Hall effect or axion
electrodynamics. Rather, the Dirac-CDW insulators display novel crystalline-
electromagnetic responses wherein charge is bound to disclination defects of
the lattice. These crystalline-electromagnetic responses are encoded in a 3D
discrete Wen-Zee term and an $R\wedge F$ term. Due to the inversion symmetry
of the DSM, there are two classes of insulating states where $\Phi=0$ and
$\Phi=2\pi$ respectively. These two insulating phases can be differentiated by
considering the total charge bound to disclination lines. These responses
naturally arise from the combination of topology and spatial symmetries that
stabilize the DSMs.

In terms of real materials, a potential material candidate for realizing the
Dirac-CDW insulator is TaTe4, a DSM that develops commensurate CDW order at
finite temperature [77]. The CDW momentum matches the momentum space
separation for a pair of Dirac points. Despite this, the CDW phase of TaTe4 is
not insulating, but rather has an 8-fold degenerate doubled Dirac point in the
folded Brillouin zone. Additionally, TaTe4 has Dirac points at other momenta
that do not strongly couple to the CDW, and are expected to remain gapless.
The gapless Dirac points prevent TaTe4 from realizing the quantized
crystalline-electromagnetic responses discussed, but it does indicate that the
combination of Dirac and CDW physics can occur in real materials.

There is also a much broader family of topological semimetals protected by
non-symmorphic crystal symmetries, for which there are many candidate host
materials and some experimentally confirmed examples [78, 79]. Understanding
the mixed crystalline-electromagnetic response of non-symmorphic topological
semimetals is an open question for future work, as the structure of non-
symmorphic symmetry fluxes is not currently understood. In addition, it would
be interesting to consider the crystalline-electromagnetic responses of higher
order Dirac semimetals, both with and without interactions. Previous studies
of higher order semimetals [80, 81, 82, 83, 49] have already revealed novel
behavior not found in their first order counterparts, but a full understanding
of the topological responses of higher order semimetals remains incomplete,
even for symmorphic symmetries.

###### Acknowledgements.

The authors thank Anton Burkov, Pranav Rao, Ryan Thorngren, and Chong Wang for
helpful discussions. J.M.M. thanks the National Science Foundation Graduate
Research Fellowship Program under Grant No. DGE - 1746047 and a startup fund
at Stanford University. L.G. thanks the Government of Canada through the
Department of Innovation, Science and Economic Development and the Province of
Ontario through the Ministry of Economic Development, Job Creation and Trade.
M.R.H., J.M.M., and T.L.H. thank ARO MURI W911NF2020166 for support. T.L.H.
thanks the U.S. Office of Naval Research (ONR) Multidisciplinary University
Research Initiative (MURI) Grant No. N00014-20-1-2325 on Robust Photonic
Materials with High-Order Topological Protection for support.

## References

  * Armitage  _et al._ [2018] N. Armitage, E. Mele, and A. Vishwanath, Reviews of Modern Physics 90, 015001 (2018).
  * Wan  _et al._ [2011] X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Physical Review B 83, 205101 (2011).
  * Burkov and Balents [2011] A. Burkov and L. Balents, Physical review letters 107, 127205 (2011).
  * Yang  _et al._ [2011] K.-Y. Yang, Y.-M. Lu, and Y. Ran, Physical Review B 84, 075129 (2011).
  * Xu  _et al._ [2011] G. Xu, H. Weng, Z. Wang, X. Dai, and Z. Fang, Physical review letters 107, 186806 (2011).
  * Young  _et al._ [2012] S. M. Young, S. Zaheer, J. C. Teo, C. L. Kane, E. J. Mele, and A. M. Rappe, Physical review letters 108, 140405 (2012).
  * Wang  _et al._ [2012] Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu, H. Weng, X. Dai, and Z. Fang, Physical Review B 85, 195320 (2012).
  * Zyuzin and Burkov [2012] A. Zyuzin and A. Burkov, Physical Review B 86, 115133 (2012).
  * Zyuzin  _et al._ [2012] A. Zyuzin, S. Wu, and A. Burkov, Physical Review B 85, 165110 (2012).
  * Wang  _et al._ [2013] Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang, Physical Review B 88, 125427 (2013).
  * Weng  _et al._ [2015] H. Weng, C. Fang, Z. Fang, B. A. Bernevig, and X. Dai, Physical Review X 5, 011029 (2015).
  * Jia  _et al._ [2016] S. Jia, S.-Y. Xu, and M. Z. Hasan, Nature materials 15, 1140 (2016).
  * Nielsen and Ninomiya [1983] H. B. Nielsen and M. Ninomiya, Physics Letters B 130, 389 (1983).
  * Burkov [2018] A. Burkov, Annual Review of Condensed Matter Physics 9, 359 (2018).
  * Klinkhamer and Volovik [2005] F. Klinkhamer and G. Volovik, International Journal of Modern Physics A 20, 2795 (2005).
  * Murakami [2007] S. Murakami, New Journal of Physics 9, 356 (2007).
  * Wilczek [1987] F. Wilczek, Physical review letters 58, 1799 (1987).
  * Wang and Zhang [2013] Z. Wang and S.-C. Zhang, Physical Review B 87, 161107 (2013).
  * Maciejko and Nandkishore [2014] J. Maciejko and R. Nandkishore, Physical Review B 90, 035126 (2014).
  * You  _et al._ [2016] Y. You, G. Y. Cho, and T. L. Hughes, Physical Review B 94, 085102 (2016).
  * Wieder  _et al._ [2020a] B. J. Wieder, K.-S. Lin, and B. Bradlyn, Phys. Rev. Res. 2, 042010 (2020a).
  * Qi  _et al._ [2008] X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Physical Review B 78, 195424 (2008).
  * Essin  _et al._ [2009] A. M. Essin, J. E. Moore, and D. Vanderbilt, Physical review letters 102, 146805 (2009).
  * Witten [2016] E. Witten, Reviews of Modern Physics 88, 035001 (2016).
  * Hughes  _et al._ [2011] T. L. Hughes, E. Prodan, and B. A. Bernevig, Physical Review B 83, 245132 (2011).
  * Varnava and Vanderbilt [2018] N. Varnava and D. Vanderbilt, Physical Review B 98, 245117 (2018).
  * Wieder and Bernevig [2018] B. J. Wieder and B. A. Bernevig, arXiv preprint arXiv:1810.02373 (2018).
  * Halperin [1987] B. I. Halperin, Japanese Journal of Applied Physics 26, 1913 (1987).
  * Kohmoto  _et al._ [1992] M. Kohmoto, B. I. Halperin, and Y.-S. Wu, Physical Review B 45, 13488 (1992).
  * Haldane [2004] F. Haldane, Physical review letters 93, 206602 (2004).
  * Zhang  _et al._ [2016] R.-X. Zhang, J. A. Hutasoit, Y. Sun, B. Yan, C. Xu, and C.-X. Liu, Physical Review B 93, 041108 (2016).
  * Chaikin  _et al._ [1995] P. M. Chaikin, T. C. Lubensky, and T. A. Witten,  _Principles of condensed matter physics_ , Vol. 10 (Cambridge university press Cambridge, 1995).
  * Ran  _et al._ [2009] Y. Ran, Y. Zhang, and A. Vishwanath, Nature Physics 5, 298 (2009).
  * Teo and Kane [2010] J. C. Teo and C. L. Kane, Physical Review B 82, 115120 (2010).
  * Juričić  _et al._ [2012] V. Juričić, A. Mesaros, R.-J. Slager, and J. Zaanen, Physical review letters 108, 106403 (2012).
  * Barkeshli and Qi [2012] M. Barkeshli and X.-L. Qi, Physical Review X 2, 031013 (2012).
  * Asahi and Nagaosa [2012] D. Asahi and N. Nagaosa, Physical Review B 86, 100504 (2012).
  * Chung  _et al._ [2016] S. B. Chung, C. Chan, and H. Yao, Scientific reports 6, 1 (2016).
  * Teo and Hughes [2017] J. C. Teo and T. L. Hughes, Annual Review of Condensed Matter Physics 8, 211 (2017).
  * Ramamurthy  _et al._ [2017] S. T. Ramamurthy, Y. Wang, and T. L. Hughes, Physical review letters 118, 146602 (2017).
  * Roy and Juričić [2021] B. Roy and V. Juričić, Physical Review Research 3, 033107 (2021).
  * Gioia  _et al._ [2021] L. Gioia, C. Wang, and A. Burkov, Physical Review Research 3, 043067 (2021).
  * Teo and Hughes [2013] J. C. Teo and T. L. Hughes, Physical review letters 111, 047006 (2013).
  * Gopalakrishnan  _et al._ [2013] S. Gopalakrishnan, J. C. Teo, and T. L. Hughes, Physical review letters 111, 025304 (2013).
  * Benalcazar  _et al._ [2014] W. A. Benalcazar, J. C. Teo, and T. L. Hughes, Physical Review B 89, 224503 (2014).
  * Li  _et al._ [2020] T. Li, P. Zhu, W. A. Benalcazar, and T. L. Hughes, Physical Review B 101, 115115 (2020).
  * Manjunath and Barkeshli [2021] N. Manjunath and M. Barkeshli, Physical Review Research 3, 013040 (2021).
  * Schindler  _et al._ [2022] F. Schindler, S. S. Tsirkin, T. Neupert, B. A. Bernevig, and B. J. Wieder, arXiv preprint arXiv:2207.10112 (2022).
  * Hirsbrunner  _et al._ [2023a] M. R. Hirsbrunner, A. D. Gray, and T. L. Hughes, arXiv preprint arXiv:2308.05796 (2023a).
  * Hirsbrunner  _et al._ [2023b] M. R. Hirsbrunner, O. Dubinkin, F. J. Burnell, and T. L. Hughes, arXiv preprint arXiv:2309.10840 (2023b).
  * Wen and Zee [1992] X.-g. Wen and A. Zee, Physical review letters 69, 953 (1992).
  * Rüegg and Lin [2013] A. Rüegg and C. Lin, Physical Review Letters 110, 046401 (2013).
  * Rüegg  _et al._ [2013] A. Rüegg, S. Coh, and J. E. Moore, Physical Review B 88, 155127 (2013).
  * Liu  _et al._ [2019] S. Liu, A. Vishwanath, and E. Khalaf, Physical Review X 9, 031003 (2019).
  * Han  _et al._ [2019] B. Han, H. Wang, and P. Ye, Physical Review B 99, 205120 (2019).
  * Manjunath and Barkeshli [2020] N. Manjunath and M. Barkeshli, arXiv preprint arXiv:2012.11603 (2020).
  * May-Mann and Hughes [2022] J. May-Mann and T. L. Hughes, Physical Review B 106, L241113 (2022).
  * Peterson  _et al._ [2021] C. W. Peterson, T. Li, W. Jiang, T. L. Hughes, and G. Bahl, Nature 589, 376 (2021).
  * Zhang  _et al._ [2022] Y. Zhang, N. Manjunath, G. Nambiar, and M. Barkeshli, arXiv preprint arXiv:2204.05320 (2022).
  * May-Mann  _et al._ [2022] J. May-Mann, M. R. Hirsbrunner, X. Cao, and T. L. Hughes, arXiv preprint arXiv:2209.00026 (2022).
  * Benalcazar  _et al._ [2019] W. A. Benalcazar, T. Li, and T. L. Hughes, Physical Review B 99, 245151 (2019).
  * Khalaf  _et al._ [2021] E. Khalaf, W. A. Benalcazar, T. L. Hughes, and R. Queiroz, Physical Review Research 3, 013239 (2021).
  * Yang and Nagaosa [2014] B.-J. Yang and N. Nagaosa, Nature Communications 5, 4898 (2014).
  * Ramamurthy and Hughes [2015] S. T. Ramamurthy and T. L. Hughes, Physical Review B 92, 085105 (2015).
  * Gioia and Wang [2022] L. Gioia and C. Wang, Phys. Rev. X 12, 031007 (2022).
  * Dubinkin  _et al._ [2021] O. Dubinkin, F. Burnell, and T. L. Hughes, arXiv preprint arXiv:2102.08959 (2021).
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Physical review letters 95, 226801 (2005).
  * Bernevig  _et al._ [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, science 314, 1757 (2006).
  * Fu and Kane [2007] L. Fu and C. L. Kane, Physical Review B 76, 045302 (2007).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Annals of Physics 16, 407 (1961).
  * Coh and Vanderbilt [2009] S. Coh and D. Vanderbilt, Phys. Rev. Lett. 102, 107603 (2009).
  * Vaidya  _et al._ [2023] S. Vaidya, M. C. Rechtsman, and W. A. Benalcazar, arXiv preprint arXiv:2304.13118 (2023).
  * Su  _et al._ [1979] W. Su, J. Schrieffer, and A. J. Heeger, Physical review letters 42, 1698 (1979).
  * Su  _et al._ [1983] W. Su, J. Schrieffer, and A. Heeger, Physical Review B 28, 1138 (1983).
  * Fang and Cano [2021] Y. Fang and J. Cano, Physical Review B 103, 165109 (2021).
  * Rao and Bradlyn [2023] P. Rao and B. Bradlyn, Physical Review B 107, 195153 (2023).
  * Zhang  _et al._ [2020] X. Zhang, Q. Gu, H. Sun, T. Luo, Y. Liu, Y. Chen, Z. Shao, Z. Zhang, S. Li, Y. Sun, _et al._ , Physical Review B 102, 035125 (2020).
  * Bradlyn  _et al._ [2016] B. Bradlyn, J. Cano, Z. Wang, M. Vergniory, C. Felser, R. J. Cava, and B. A. Bernevig, Science 353, aaf5037 (2016).
  * Lv  _et al._ [2021] B. Lv, T. Qian, and H. Ding, Reviews of Modern Physics 93, 025002 (2021).
  * Lin and Hughes [2018] M. Lin and T. L. Hughes, Physical Review B 98, 241103 (2018).
  * Ghorashi  _et al._ [2020] S. A. A. Ghorashi, T. Li, and T. L. Hughes, Phys. Rev. Lett. 125, 266804 (2020).
  * Wang  _et al._ [2020] H.-X. Wang, Z.-K. Lin, B. Jiang, G.-Y. Guo, and J.-H. Jiang, Phys. Rev. Lett. 125, 146401 (2020).
  * Wieder  _et al._ [2020b] B. J. Wieder, Z. Wang, J. Cano, X. Dai, L. M. Schoop, B. Bradlyn, and B. A. Bernevig, Nat. Comms. 11, 627 (2020b).

##  Appendix A Disclination charge parity and the disclination filling anomaly
in layered systems

In this Appendix, we will discuss why the disclination charge parity is not a
meaningful quantity for inversion-symmetric insulators, as it depends on the
choice of inversion center for systems with open boundaries. However, as we
shall see, the difference between the disclination charge parity for open and
periodic boundary conditions (referred to in the main text as the disclination
filling anomaly) does not depend on the inversion center, and is therefore a
more physically relevant quantity. We shall establish this by considering
layered insulators. Although we are only considering specific systems, the
dependence of the disclination charge parity on the choice of inversion center
is general. We recall from the main text that what we mean by a layered
insulator is one that is adiabatically connected to a limit of decoupled
layers where electrons are not hybrdized between different unit cells in the
$z$-direction.

To this end, let us consider an inversion-symmetric 3D insulator composed of
$N_{z}$ layers of 2D insulators, stacked along the $z$-direction. We further
take the system to have open boundary conditions in the $z$-direction. For
such a system, there are two choices of inversion symmetry that differ with
respect to their inversion center. First, there is inversion symmetry where
the inversion center is in a given layer, such that this layer maps to itself
under inversion symmetry. This was referred to as site-centered inversion
symmetry in the main text. Second, there is inversion symmetry where the
inversion center is in-between two adjacent layers. This was referred to as
bond-centered inversion symmetry in the main text. For site-centered inversion
symmetry and open boundaries, $N_{z}$ must be odd, since under inversion
symmetry, one layer maps to itself, while all other layers must map onto an
inversion symmetry related partner. Similarly, for bond-centered inversion
symmetry, $N_{z}$ must be even, since all layers map onto an inversion
symmetry related partner. The two inversion-symmetric stacking configurations
are shown in Fig. 8(a),(b).

Figure 8:  An inversion-symmetric stack of 2D layers (green) for site-centered
inversion symmetry (a) and bond-centered inversion symmetry (b). For site-
centered inversion symmetry, one layer maps onto itself under inversion
symmetry, while other layers map onto a partner, such that the total number of
layers, $N_{z}$ is odd. For bond-centered inversion symmetry, all layers map
onto a partner, such that the total number of layers, $N_{z}$ is even.  Figure
9: The charge bound to a $\Theta_{F}=\pi/2$ disclination for a Dirac-CDW
insulator as a function of the integration radius distance $r$ with $Q=\pi/2$,
$b_{xy}=1$, $\Delta_{\text{surf}}=0.25$, $\Delta_{s}=0$, $\Delta_{b}=0.5$,
$N_{x}=15$, and $N_{z}=14$. The radius $r$ is scaled such that the farthest
point from the disclination core is at $r=1$. Note that the $r$-axis only
extends from $0.1$ to $0.6$.

We now consider the disclination response of the layered insulator. For a 3D
layered insulator, the disclination response of the full system is simply the
sum of responses of the 2D layers. Each 2D layer is described by a 2D dWZ
term, with discrete shift $\mathcal{S}_{\text{2D}}$. For spin-1/2 fermions
with TRS, $\mathcal{S}_{\text{2D}}$ is quantized as an even integer. Because
of the 2D dWZ response, a $\Theta_{F}$ disclination line binds charge
$Q_{\text{disc-2D}}=\mathcal{S}_{\text{2D}}\frac{\Theta_{F}}{2\pi}\mod(2)$ on
each layer, where the $\mod(2)$ factor reflects that it is possible to locally
add a Kramers degenerate pair of electrons to the disclination core. The total
disclination charge of the full 3D system is then
$Q_{\text{disc-3D}}=N\mathcal{S}_{\text{2D}}\frac{\Theta_{F}}{2\pi}\mod(2)$.
For a finite size system of spin-1/2 fermions with time-reversal symmetry
(TRS), the disclination charge parity is equal to the total charge on
disclination line $\mod(2\frac{\Theta_{F}}{\pi})$,

$Q_{\text{disc-3D}}\,\,\text{mod}\left(2\frac{\Theta_{F}}{\pi}\right)=\mathcal{S}_{\text{2D}}N_{z}\frac{\Theta_{F}}{2\pi}\,\,\text{mod}\left(2\frac{\Theta_{F}}{\pi}\right),$
(52)

where we have used the fact that $\Theta_{F}$ is a multiple of $2\pi/n$ for
$C_{n}$ symmetric insulators. The disclination charge parity is therefore zero
when $N_{z}\mathcal{S}_{\text{2D}}/2$ is even, and non-zero when
$N_{z}\mathcal{S}_{\text{2D}}/2$ is odd. From this, we can conclude the
following: when $\mathcal{S}_{\text{2D}}/2$ is even, the disclination charge
parity always vanishes, but when $\mathcal{S}_{\text{2D}}/2$ is odd, the
disclination charge parity vanishes when $N_{z}$ is even and is non-zero when
$N_{z}$ is odd. So when $\mathcal{S}_{\text{2D}}/2$ is odd, the disclination
charge parity always vanishes for bond centered inversion symmetry, and is
always non-zero for site-centered inversion symmetry. This occurs for an
insulator composed of quantum spin Hall layers, each of which have
$\mathcal{S}_{\text{2D}}=2$.

Having established that the disclination charge parity depends on the choice
of inversion center, let us now consider the disclination filling anomaly. To
do this, we will take the layered system and “sew” the top and bottom layers
to one another such that the system has periodic boundaries in the
$z$-direction. For a layered system, the sewing procedure is trivial by
definition and will not change the total disclination charge, or, by
extension, the disclination charge parity. The disclination filling anomaly
will therefore necessarily vanish for all layered systems, regardless of the
choice of inversion center or the value of $\mathcal{S}_{\text{2D}}$.

##  Appendix B Details of the numerics

Here we discuss the details of how we calculate the charge bound to a
disclination and the extrapolations in system size that we performed. To
calculate the charge bound to a disclination, we sum the charge on all lattice
sites in each layer that are within some radius from the disclination core. In
Fig. 9 we plot the charge bound to a disclination as a function of the
integration radius for a Dirac-CDW insulator with $Q=\pi/2$, $N_{x}=15$,
$N_{z}=14$, $b_{xy}=1$, $\Delta_{\text{surf}}=0.25$, $\Delta_{s}=0$,
$\Delta_{b}=0.5$, and both $\theta_{b}=0$ and $\theta_{b}=\pi$. The radius $r$
is scaled such that the farthest point from the disclination core in the
lattice is at $r=1$. The disclination charges approach the theoretical
predictions, marked by dashed horizontal lines, for roughly $r>0.3$,
indicating the exponential localization of the disclination charge. The
disclination charge approaches zero for both small and large $r$, and the
largest values of the disclination charge generally approach the theoretical
predictions, except for very small system sizes. However, the integration
radius at which the disclination charge most closely obtains the theoretical
prediction varies significantly as a function of system size and other
parameters. As such, we always calculate the disclination charge using a range
of integration radius and report the maximal absolute value.

To diminish the impact of finite-size effects, we calculated the charges for
each case enumerated in Table 1 over a range of system sizes and extrapolated
to the infinite system-size limit. We found that effects arising from the
finite extent of the lattice in the $z$-direction were minimal, so for each
calculation we fixed $N_{z}$ and varied $N_{x}$, the side length of the
disclinated lattice. In Table 2 we plot the disclination charges as a function
of $N_{x}$ along with fits to a decaying exponential,
$|Q_{\text{disc}}(N_{x})|=Q_{\infty}+be^{-cN_{x}}$. We plot the infinite
system size limit, $Q_{\infty}$, with red dashed lines and report the values
and uncertainties in the insets. The values reported in Table 1 are the values
of $Q_{\infty}$ rounded to the nearest half integer, which produces less than
a $1\%$ error in all cases.

$\mathcal{K}a_{z}$ | OBC | PBC
---|---|---
|  $\theta_{s/b}=0$ |  $\theta_{s/b}=\pi$ |  $\theta_{s/b}=0$ |  $\theta_{s/b}=\pi$
$\pi/2$ |  |  |  |
$\pi/3$ |  |  |  |
$\pi/2\sqrt{2}$ |  |  |  |
Table 2: System size extrapolations of the charge bound to a
$\Theta_{F}=\pi/2$ disclination for $\mathcal{K}=\pi/2$, $\pi/3$, and
$\pi/2\sqrt{2}$ with $\theta_{s/b}=0$ and $\pi$ for both open and periodic
boundary conditions. For each calculation we used $b_{xy}=1$,
$\Delta_{\text{surf}}=0.25$, and $\Delta_{s}=0$. The system size along the
disclination was set to $N_{z}=14$ for $\mathcal{K}=\pi/2$ and $N_{z}=15$ for
$\mathcal{K}=\pi/3$ and $\pi/2\sqrt{2}$. The solid red lines are fits to a
decaying exponential, $|Q_{\text{disc}}(N_{x})=Q_{\infty}+be^{-cN_{x}}$. The
value and uncertainty of the large-$N_{x}$ limit, $Q_{\infty}$, is reported in
the legend and plotted with a red dashed line.