††thanks: These authors contributed equally.††thanks: These authors
contributed equally.

# Anomalous crystalline-electromagnetic responses in semimetals

Mark R. Hirsbrunner  Department of Physics and Institute of Condensed Matter
Theory, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
Oleg Dubinkin  Department of Physics and Institute of Condensed Matter Theory,
University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Fiona J.
Burnell  Department of Physics, University of Minnesota Twin Cities, MN 55455,
USA  Taylor L. Hughes  Department of Physics and Institute of Condensed Matter
Theory, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA

###### Abstract

We present a unifying framework that allows us to study the mixed crystalline-
electromagnetic responses of topological semimetals in spatial dimensions up
to $D=3$ through dimensional augmentation and reduction procedures. We show
how this framework illuminates relations between the previously known
topological semimetals, and use it to identify a new class of quadrupolar
nodal line semimetals for which we construct a lattice tight-binding
Hamiltonian. We further utilize this framework to quantify a variety of mixed
crystalline-electromagnetic responses, including several that have not
previously been explored in existing literature, and show that the
corresponding coefficients are universally proportional to weighted momentum-
energy multipole moments of the nodal points (or lines) of the semimetal. We
introduce lattice gauge fields that couple to the crystal momentum and
describe how tools including the gradient expansion procedure, dimensional
reduction, compactification, and the Kubo formula can be used to
systematically derive these responses and their coefficients. We further
substantiate these findings through analytical physical arguments, microscopic
calculations, and explicit numerical simulations employing tight-binding
models.

##  I Introduction

Topological responses are a key manifestation of electronic topology in
solids. Celebrated examples such as the integer quantum Hall effect [1, 2, 3]
and axion electrodynamics [4, 5] have paved the way for a broader exploration
of topological response phenomena in insulating systems. As of now, a wide-
variety of phenomena that are directly determined by the electronic topology
have been considered, including thermal response [6, 7], geometric response
[8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24], and
electric multipole response [25, 26, 27]. These responses are robust features
of topological insulators (TIs), and topological phases in general, and are
often described by a quantized response coefficient, e.g., the integer Hall
conductance [1, 2, 3], or the quantized magneto-electric polarizability [5,
28, 29].

Interestingly, certain distinctive features of response of topological Weyl or
Dirac semimetals can be described by response theories that are closely
analogous to those of topologically insulating phases, albeit with
coefficients that are determined by the momentum-space and energy locations of
the point or line nodes [30, 31, 32, 33, 34, 35]. For point-node semimetals,
the relevant response coefficients are momentum-energy vectors determined as a
sum of the momentum and energy locations of the point-nodes weighted by their
chirality (or by their helicity, for Dirac semimetals), yielding a momentum-
energy space dipole. For example, the low-energy, nodal contribution to the
anomalous Hall effect tensor of a 3D Weyl semimetal is determined by the
momentum components of this momentum-energy dipole vector.

The quasi-topological response coefficients of topological semimetals are not
strictly quantized since they can be continuously tuned with the nodal
momenta. However, the forms of the responses share many features with
topological insulators in one lower dimension, or perhaps more precisely, with
weak topological insulators in the same dimension [36, 37]. Indeed,
topological semimetals and weak topological insulators both require discrete
translation symmetry to be protected and both are sensitive to translation
defects such as dislocations [38]. Interestingly, the connection to
translation symmetry has motivated recent work which recasts many previously
proposed topological responses of these systems as couplings between the
electromagnetic gauge field $A_{\mu}$ and gauge fields for translations
$\mathfrak{e}_{\mu}^{a},$ where $\mu$ runs over spacetime indices, and $a$
runs over the spatial directions in which translation symmetry exists. This
insight has also led to the development of new response theories that are just
beginning to be understood [18, 20, 17, 19, 24].

Motivated by these previous results and our recent related work on higher rank
chiral fermions [17, 39, 24], here we study the topological responses of 1D,
2D, and 3D topological semimetals coupled to electromagnetic and strain
(translation gauge) fields. In addition to the well-studied dipole case
mentioned above, we also study cases where point-nodes have momentum-energy
quadrupole or octupole patterns. Our approach allows us to make clear
connections between a wide variety of response theories across dimensions, and
clarifies relationships between many of the response theories we discuss. We
find that the chirality-weighted momentum-energy _multipole_ moments of the
semimetals determine new types of quasi-topological responses to
electromagnetic fields and strain. We are able to explicitly derive many of
these responses from Kubo formula calculations (sometimes combined with
dimensional reduction procedures [5]). Using these results we explicitly study
these families of response theories using lattice model realizations. We also
extend our results to the responses of nodal line semimetals (NLSMs) and
construct a new type of NLSM with an unusual crossed, cage-like nodal
structure.

Our article is organized as follows. In Sec.II we provide an overview of and
intuition about the response theories that will be discussed in more detail,
and in model contexts, in later sections. In Sec. III we derive a family of
effective actions that describe mixed crystalline-electromagnetic responses in
various spatial dimensions. From here we proceed in Sec. IV by presenting
concrete lattice models and explicit numerical calculations that realize and
demonstrate the mixed responses in D$=1,2,3$. We conclude in Sec. V by
discussing possible extensions to future work, and potential pathways to
experimental observation of some of the described phenomena.

##  II Overview of Response Theories

Figure 1: (a) A dimensional hierarchy of theories describing responses of
strong topological insulators. The theories are related by dimensional
reduction ($\theta$ symbol, green arrow) [5], taking the boundary response
($(1/2)\partial$ symbol, purple arrow), or adiabatic pumping ( symbol, red
arrow) [40]. (b) A dimensional hierarchy of insulating systems with mixed
crystalline-electromagnetic responses. The theories are related by stacking
(layer symbol, dark red arrow) and cutting (scissor symbol, blue arrow). (c) A
family tree of dimensional hierarchies establishing connections between
responses of strong TIs and insulators with mixed crystalline-electromagnetic
responses. (d) Illustrations representing the nature of the phases
constituting the hierarchy depicted in (c). (i) A single isolated charge. (ii)
A line of charges forming a lattice. (iii) An insulating chain having a
quantized charge polarization. (iv) A two-dimensional lattice of charges. (v)
A two-dimensional weak topological insulator where polarized chains are
stacked transverse to their polarization. (vi) A two-dimensional Chern
insulator having chiral edge states indicated by red arrows. (vii) A three-
dimensional lattice of charges. (viii) A three-dimensional lattice built from
a two-dimensional array of polarized chains; alternatively, a stack of two-
dimensional weak topological insulators. (ix) A three-dimensional stack of
Chern insulators forming a time-reversal breaking weak topological insulator.
(x) A three-dimensional strong topological insulator with surface Dirac cones.

The systems we consider in this article all exhibit $U(1)$ charge conservation
and discrete translation symmetry in at least one spatial direction. In the
presence of these symmetries we can consider the responses to background field
configurations of the electromagnetic gauge field $A_{\mu}$ and a collection
of translation gauge fields $\mathfrak{e}_{\mu}^{a}.$ For example, if the
system exhibits translation symmetry in the $x$-direction, then we can
consider coupling the system to the field $\mathfrak{e}_{\mu}^{x}.$ Our goal
is to study low-energy response theories of electrons coupled to translation
and electromagnetic gauge fields.

Since most readers are likely less familiar with the translation gauge fields
$\mathfrak{e}_{\mu}^{a}$ than the electromagnetic field $A_{\mu}$, we will
briefly review the nature of these fields as they appear in our work. In a
weakly deformed lattice, $\mathfrak{e}^{a}$ is given by

$\mathfrak{e}_{j}^{a}=\delta_{j}^{a}-\frac{\partial u^{a}}{\partial x^{j}},$
(1)

where the Kronecker $\delta_{j}^{a}$ encodes the fixed reference lattice
vectors, $u^{a}$ is the lattice displacement, and $\frac{\partial
u^{a}}{\partial x^{j}}$ is the distortion tensor [41]. The fields
$\mathfrak{e}_{j}^{a}$ in Eq. (1) are reminiscent of gauge fields (see, e.g.,
[42]): from Eq. (1) we immediately see that line integrals of $\mathfrak{e}$
describe lattice dislocations since $\oint\tfrac{\partial u_{a}}{\partial
x_{i}}dx^{i}=b^{a}$, where $b^{a}$ is the net Burgers vector of all the
dislocations inside the loop [41]. This points to an analogy with the
configurations of the usual electromagnetic field. The analog of magnetic
fields derived from $\mathfrak{e}_{\mu}^{a}$ essentially encodes
configurations of dislocations, each with an amount of flux equal to the
corresponding Burgers’ vector. Additionally, electric fields are time-
dependent strains. In earlier work, e.g., Ref. 11, these fields could have
been called frame-fields, but crucially the translation gauge fields encode
only the translation/torsional part of the geometric distortion, whereas the
frame fields also carry rotational information. In keeping with previous
literature, here we will call the set of (Abelian) fields
$\mathfrak{e}_{\mu}^{a}$ translation “gauge” fields, by analogy of their
relationship to translation “fluxes” (i.e. lattice defects). This language is
convenient because, as we will see, actions describing the response to such
lattice fluxes are invariant under (vector-charge) gauge transformations of
the $\mathfrak{e}_{\mu}^{a}$ fields.

A second way in which we will use the close analogy between
$\mathfrak{e}_{\mu}^{a}$ and electromagnetic gauge fields is through the
lattice analog of the usual Aharonov-Bohm effect (holonomy), in which a
charged particle encircles a magnetic flux of the gauge field. In the
electromagnetic case, a charged particle moving around a magnetic flux
generates a $U(1)$ phase factor. For the translation gauge field, taking a
particle around a translation magnetic flux having Burgers vector $b^{a}$
generates a translation operator by the displacement $b^{a}$. For particles
with a fixed translation charge, i.e., a fixed momentum, this generates a
momentum-dependent $U(1)$ phase factor. This will lead us to introduce
momentum-dependent Peierls factors when performing some lattice calculations.
To complement this discussion, in Appendix A we show more explicitly how
translation symmetry can be “gauged” under a teleparallel constraint of the
underlying system geometry. A very similar approach has been used to study the
effects of strain on graphene [43, 44, 45] and other semimetallic systems [46,
47, 48, 49, 50], where strain can play the role of a valley-dependent magnetic
field.

For our purposes there are many ways in which $\mathfrak{e}_{\mu}^{a}$ can be
treated on equal footing with the electromagnetic gauge field. However, there
are some important distinctions. First, the fields $\mathfrak{e}_{j}^{a}$ in
Eq. (1) are not true gauge fields. This becomes important when considering the
possible response actions: while the total charge of a system strictly
conserved, momentum conservation is not similarly inviolable (see e.g. [48,
49, 50] for some interesting physical consequences of this distinction).
Second, responses involving $\mathfrak{e}_{\mu}^{a}$ are predicated on the
existence of translation symmetry. Thus, if the response is characterized by a
boundary effect or a response to a flux/defect, we must be careful to ensure
that (at least approximate) translation symmetry is maintained in order to
connect the coefficient of the response action to explicit model calculations.
Indeed, some responses are not well-defined unless configurations that
maintain translation symmetry are used. This is unlike the electromagnetic
response for which $U(1)$ charge symmetry is maintained independently of the
geometry and gauge field configuration. Other important distinctions have been
discussed in recent literature that has begun putting the gauging of discrete
spatial symmetries on firmer ground [51, 52, 53]. One important distinction is
that the translation gauge fields correspond to a discrete gauge symmetry
$\mathbb{Z}_{N_{a}}$, where $N_{a}$ is the number of unit cells in the $a$-th
direction. This discreteness can play an important role in the topological
response properties [18], but we will not focus on this aspect in our work.

Using this framework, our goal is to consider the low-energy responses of
electrons to the background electromagnetic and translation gauge fields.
Given a translationally invariant Bloch Hamiltonian $H$, the response theories
we consider can, in principle, be derived from correlation functions of the
electromagnetic current

$j^{\mu}=e\frac{\partial H}{\partial k_{\mu}},$ (2)

and the crystal momentum current

$\mathcal{J}^{\mu}_{a}=\hbar k_{a}\frac{\partial H}{\partial k_{\mu}},$ (3)

where the former couples to $A_{\mu}$ and the latter to
$\mathfrak{e}_{\mu}^{a}$ (see App. A for more details for the latter). Indeed,
we take exactly this approach in Section III to derive response actions for 2D
and 3D systems. While our explicit derivations are important for precisely
determining the coefficients of the response actions we study, it will be
helpful to first motivate the overarching structure that connects a large
subset of these response theories. We also note that alternative approaches to
determining some of the response actions we discuss have been proposed in
Refs. 18, 19, 20, and where the results overlap with ours, they agree.

To understand the connections between the response theories we study, it is
useful to begin by reviewing the well-known dimensional hierarchy of response
theories of strong topological insulators [5]. We show the general structure
in Fig. 1(a), where the response terms are built solely from the
electromagnetic gauge field. Furthermore, Chern-Simons and $\theta$-term
response actions appear in even and odd spatial dimensions, respectively.
There are a number of connections between the theories in different
dimensions, and we will now review three of them. First, a Chern-Simons action
in D spatial dimensions can be dimensionally reduced to a $\theta$-term action
in (D-1)-dimensions by compactifying one spatial direction [54, 55]. The
(D-1)-dimensional system can also represent a TI if the value of $\theta$ is
quantized to be $0,\pi$ by a symmetry that protects the (D-1)-dimensional
topological insulator [5]. Second, one can consider the reverse process where
quantized adiabatic pumping [40] in (D-1)-dimensions will convert a
$\theta$-term action to a D-dimensional Chern-Simons action. Finally, a
$\theta$-term action for a (D-1)-dimensional topological insulator exhibits a
half-quantized (D-2)-Chern-Simons response on boundaries where $\theta$ jumps
by $\pi.$ These general relationships are summarized in Fig. 1(a) where each
type of relationship is color- and symbol-coordinated.

Next we can consider a less familiar set of relationships in Fig. 1(b) between
_gapped_ theories with mixed crystalline-electromagnetic responses arising
from effective actions having both $A_{\mu}$ and $\mathfrak{e}_{\mu}^{a}$
fields. We emphasize that the precise relationships we refer to in Fig. 1(b)
are for gapped systems where the coefficients of the actions are quantized. In
contrast, for the majority of this article we will focus on the quasi-
topological responses of _gapless_ systems which take similar forms, but with
non-quantized coefficients. Remarkably, many of the actions we discuss for
insulators can be generalized to the non-quantized case. For semimetals,
however, the dimensional relationships we point out are more akin to physical
guides than a precise prescription for deriving matching coefficients in-
between dimensions.

With this caveat in mind, let us consider the family of theories in Fig. 1(b).
In 0D we can consider the response action for a gapped system of
electrons,$S[A]=Q\int dtA_{0},$ which represents a system with charge
$Q=eN_{e}$ where $N_{e}$ is the (integer) number of electrons. If we imagine
stacking these 0D systems in a discrete, translationally invariant lattice in
the $x$-direction, then we will generate a line of charges. Indeed, stacking
produces the response for a translation invariant line-charge density which is
captured by the next action in the sequence in Fig. 1(b), i.e., $Q\int A\wedge
e^{x}=Q\int dxdt(A_{0}\mathfrak{e}_{x}^{x}-A_{x}\mathfrak{e}_{0}^{x}).$ In
this action the first term represents the charge density along the line, while
the second term represents a current generated if the lattice of charges is
moving. The latter consequence becomes manifest in the weakly distorted
lattice limit since the current is proportional to the displacement rate:
$j\sim\mathfrak{e}_{0}^{x}\sim\tfrac{\partial u^{x}}{\partial t}$.

We can also imagine a reverse process where we are given a translationally
invariant line of charge at integer filling and cut out a single unit cell.
Since the system is gapped and translation invariant, this will result in a
move in the opposite direction in Fig. 1(b), i.e., from
$A\wedge\mathfrak{e}^{x}$ in 1D to $A$ in 0D with the same _integer_
coefficient $Q.$ We can use this example to highlight our caveat about gapped
vs. gapless systems mentioned above. That is, while it is reasonable to have a
1D gapless system with non-quantized (i.e., non-integer) charge (per unit
cell) described by the 1D action, the cutting procedure will not work properly
at non-integer filling since the result will be a 0D point with a fractional
charge.

In comparison to the response sequence for strong TIs, we see that stacking is
the analog of pumping for the translation gauge field 111As mentioned, this
analogy is precise for gapped systems. For gapless systems the analogy
predicts the correct form of the action, but does not uniquely determine the
coefficient.. Indeed, while pumping adds an extra electromagnetic gauge field
factor $A,$ stacking adds an extra translation gauge field
$\mathfrak{e}^{D+1}$ where D+1 is the stacking direction. As a result, given
any action in the strong TI sequence, we can stack copies to get the response
action of a primary weak TI (stacks of co-dimension-$1$ strong TIs, e.g.,
lines stacked into 2D) by adding a wedge product with $\mathfrak{e}^{D+1}.$ We
can push the stacking idea further to generate secondary weak TIs (stacks of
co-dimension-$2$ strong TIs, e.g., lines stacked into 3D) by a wedge product
with $\mathfrak{e}^{D+1}\wedge\mathfrak{e}^{D+2}$ and so on.

The stacking and cutting procedures are not the only relationships between the
response theories in Fig. 1(b). Just as in the strong TI sequence, we can find
connections between the boundary properties of some D-dimensional systems and
the bulk response of a (D-1)-dimensional system. For example, the 2D response
action in Fig. 1(b) represents the response of a stack of Su-Schrieffer-Heeger
chains (SSH) [57], each with a quantized polarization of $e/2.$ The boundary
of such a 2D system is a line of charge on the edge, albeit with a density of
$e/2$ electrons per unit cell on the edge line instead of the integer density
we would get by stacking integer-filled 0D points. As such, the boundary of
the 2D $A\wedge d\mathfrak{e}^{x}$ action represents a line-charge described
by the action $A\wedge\mathfrak{e}^{x},$ but with a half-integer coefficient.

Now we can combine the dimensional relationships in the sequences of both Fig.
1(a) and (b) to make a family tree of related theories. We show a tree in Fig.
1(c) that includes response actions in 0, 1, 2, and 3 spatial dimensions. In
0D we have only an integer electron charge response that couples to $A_{0}.$
For 1D, we can either stack charges to form a line of charge (upper branch),
or consider an electrically polarized TI (lower branch) where the charge is
split in half and moved to opposing ends of the chain while the interior
remains neutral (so to speak). In 2D we can stack line charges to get a plane
of charge (top branch), stack 1D polarized TIs to get a weak TI (middle
branch), or pump charge in a 1D TI to generate a 2D Chern insulator (bottom
branch).

In 3D the set of responses is richer. We can stack plane charges to generate a
3D volume of charges (top branch), stack Chern insulators to get a 3D primary
weak TI (second from bottom branch), or stack 2D weak TIs to get a 3D
secondary weak TI built from 1D polarized wires (second branch from top). The
other well-known possibility is the magneto-electric response for a 3D strong
TI [5, 28] (bottom branch). Although it is not shown, this theory is related
to a 4D quantum Hall system via pumping (3D to 4D) or dimensional reduction
(4D to 3D) [5]. The final option we consider, which is the middle branch
enclosed by a dotted rectangle, is $\int dA\wedge d\mathfrak{e}^{a}.$ This
response theory has not been previously studied in detail. This theory is a
total derivative, and yields a gapped boundary with an electric polarization
(e.g., a stack of SSH chains on the boundary). This is reminiscent of an
electric quadrupole (higher-order) response [27, 58], and we will explore this
connection further in Sec. III.5.

While this discussion has centered on gapped systems, our primary focus is on
gapless topological semimetals. Importantly, each of the actions that contains
a translation gauge field in the family tree in Fig. 1(c) can also represent a
contribution to the response of various types of metals or topological
semimetals [30, 31, 32, 33, 34, 18, 20, 17, 19, 24]. This is because many
semimetals can be generated by translation-invariant stacking of lower
dimensional topological phases. Since the momentum $k^{a}$ in the stacking
direction is conserved, one can consider adding up the set of topological
response terms for each gapped $k^{a}.$ A semimetal represents a scenario
where the coefficients of these topological terms at each $k^{a}$ are
quantized and have discrete jumps where $k^{a}$ hits a nodal point. For
example, the 2D electric polarization response of a stack of 1D TIs becomes
the response of a 2D Dirac semimetal if the wires forming the stack are
coupled strongly enough to close the insulating gap [34]. In the presence of
reflection symmetry, each momentum in the stacking direction has a quantized
charge polarization that jumps when the momentum hits a gapless 2D Dirac
point. Additionally, the 3D response of a stack of Chern insulators becomes
the non-quantized anomalous Hall effect response of a time-reversal breaking
Weyl semimetal where each fixed-$k$ plane that does not intersect a Weyl point
carries a quantized Chern number that jumps at a Weyl point [30, 31, 32, 33],
and so on. While many of these response theories have been discussed in detail
before, only a few works have highlighted the contributions from the
translation gauge fields [59, 12, 60, 18, 20, 17, 19, 24]. As such, a large
fraction of our paper will be devoted to both the explicit derivations of the
response coefficients of the actions in Fig. 1(c) that have couplings to the
translation gauge fields (Sec. III), and to the explicit calculations of the
physical response phenomena in representative model systems (Sec. IV).

Figure 2: (a) Fermi-surfaces of a 3D time-reversal invariant Weyl semimetal
with a quadrupole Weyl node configuration. Red and blue colors denote positive
and negative Berry-curvature respectively. The associated action has a
coefficient matrix $\mathcal{Q}_{ab}$ which is symmetric and proportional to
the Weyl-node quadrupole moment. For this configuration the coefficients
$\mathcal{Q}_{xx}$ and $\mathcal{Q}_{yy}$ are non-vanishing. (b) Similar to
subfigure (a) except it is the Fermi surfaces for a 2D Dirac semimetal having
four Dirac nodes in a quadrupole pattern. The action is described by a
symmetric matrix of coefficients $\mathcal{Q}_{ab}.$ (c) the Fermi surface of
an unusual cage-like nodal line semimetal built from stacking the Dirac node
quadrupole semimetal in subfigure (b). The action has a set of coefficients
$\mathcal{B}_{ab,c}$ which is anti-symmetric in $a$ and $b.$ Heuristically the
action in (b) can generate the action in (a) by adiabatic pumping, or can
generate the action in (c) by stacking.

Before we move on to more explicit derivations, we want to motivate three
additional response theories we will study that lie outside the family tree in
Fig. 1(c). As mentioned above, a remarkable feature of the response actions of
point-node semimetals is that their coefficients are determined from the
energy-momentum locations of the nodal points. Indeed, for the relevant
response actions in Fig. 1(c), the coefficients are obtained as a chirality-
weighted momentum dipole moment of the point-nodes (note that Dirac points do
not have a chirality, nevertheless there is a signed quantity that plays the
same role). Interestingly, recent work on rank-2 chiral fermions and Weyl
semimetals with a chirality-weighted momentum _quadrupole_ moment [17, 18, 19,
24] has unveiled a new set of response theories. This category of theories has
actions that include factors of more than one translation gauge field of the
same type (e.g., $\mathfrak{e}^{a}\wedge d\mathfrak{e}^{b}$, where $a=b$), and
as such, does not appear in the family tree in Fig. 1(c). This also implies
that the translation gauge field factors in these response theories cannot be
obtained by the conventional stacking of lower dimensional systems that we
discussed above, since stacking produces wedge products with distinct
translation gauge fields. We could also construct related higher dimensional
theories (and lower dimensional theories if we considered both space and time
translational gauge fields) to form an additional connected tree of theories,
but we leave further discussion of those extensions to future work.

To give some explicit examples we show three response theories that follow
this pattern in Fig. 2. Fig. 2(a) shows the Fermi surface structure of a 3D
time-reversal invariant Weyl semimetal having a Weyl node quadrupole moment.
The response action of this system is a mixed response between electromagnetic
and translation gauge fields, and the inset in the Fermi-surface figure lists
which coefficients $\mathcal{Q}_{ab}$ are non-vanishing. Some details of this
response were discussed in Refs. 17, 18, 24, the former of which connected the
response to rank-2 chiral fermions on the surface of the 3D Weyl semimetal.
Fig. 2(b) shows the Fermi surface structure of a 2D Dirac semimetal having a
Dirac node quadrupole structure. This response represents a momentum current
responding to a translation gauge field (e.g., a strain configuration). Its
form shares some similarities with the torsional Hall viscosity [11, 61, 59,
62, 63], though a precise connection will be left to future work. Finally, in
Fig. 2(c) we show the Fermi-surface for an unusual nodal line semimetal formed
from stacking the Dirac node quadrupole semimetal of Fig. 2(b). While one
might have naively expected two independent Fermi rings, we instead find a new
type of Fermi-surface structure where the Fermi lines join at two crossing cap
regions to form a cage.

The symbols on the right-hand-side of Fig. 2 indicate the connections between
these theories: (i) the response of the nodal line structure is just a stacked
version of the 2D Dirac node quadrupole semimetal response from Fig. 2(b), and
(ii) one can heuristically consider the four-node Weyl response in Fig. 2(a)
to be a dimensional extension of the response in Fig. 2(b) via pumping.

##  III Effective response actions

Now that we have described the forms of the various response actions of
interest, we will spend this section determining their coefficients. All of
the response actions in Fig. 1(c) that contain only electromagnetic gauge
fields represent insulators, and their coefficients have been studied in
detail (e.g., see Ref. 5). The actions containing translation gauge fields can
represent insulators or gapless systems, and the two can often be
distinguished by the values of the coefficients. That is, for insulators we
expect the coefficients to be quantized in some units (in even spatial
dimensions they are quantized in the presence of some symmetry), while for
topological semimetals we expect the coefficients to be a tunable function of
the momentum and energy locations of the nodal points or lines. Interestingly,
some of the response coefficients for metals/semimetals can take the same
values allowed for an insulator, although this would typically require fine-
tuning, or extra symmetry. For example, a 1D system can have compensating
particle and hole Fermi surfaces such that the total filling is an integer, as
one would find in an insulator, yet the system is still gapless. In such a
case we will show that the system has additional response terms that have
coefficients that are incompatible with a gapped insulator.

Our focus will be on 2D Dirac, 3D Weyl, and 3D nodal line semimetals, and
before we begin our derivations it is important to acknowledge a key
qualitative difference between these types of topological semimetals. Namely,
we recall that 2D topological Dirac semimetals and 3D nodal line semimetals
require symmetry (the composite $\mathcal{T}\mathcal{I}$ symmetry) to
guarantee the local stability of the gapless points/lines in momentum space.
This is inherently different from the case of Weyl semimetals in 3D where Weyl
nodes require no extra symmetry to protect them against perturbations. Indeed,
a Weyl node can be gapped out only by bringing another Weyl node of opposite
chirality to the same point in the Brillouin zone. A similar story applies to
(semi)metallic systems in 1D: each gapless point has a well-defined chirality
defined as the sign of the Fermi velocity, and a gap can be opened only after
overlapping Fermi points of opposite chiralities.

This distinction in symmetry protection is important for the response theories
describing Dirac and Weyl semimetals as it reflects the well-known structure
of anomalies in even and odd spatial dimensions. Furthermore, it will impact
our strategy for deriving the response coefficients for these systems. As an
example, the response properties of 2D Dirac semimetals can be determined
straightforwardly from the Kubo formula if we first apply a symmetry-breaking
perturbation that weakly gaps out the nodes. The resulting insulator response
can then be taken to the semimetallic limit if we tune the perturbation to
zero. Hence, the effective response action for such systems can be obtained by
treating the system as an insulator and applying the Kubo formula, or more
generally, a gradient expansion procedure. This method can be applied to 2D
and 4D Dirac semimetals, and consequently 3D nodal line semimetals since they
are just stacks of 2D Dirac semimetals. For such semimetals we actually have a
choice of what symmetry to break, e.g., inversion or time-reversal. Which one
we need to break depends on the nodal configuration and the action we are
intending to generate. For example, in the case of a 2D Dirac semimetal with a
pair of nodes, breaking time-reversal is well-studied and generates a quantum
Hall response via a Chern-Simons term. However, breaking inversion symmetry is
relatively less-studied and generates a mixed Chern-Simons response between an
electromagnetic gauge field and a translational gauge field. This is
corroborated by the fact that the electromagnetic Chern-Simons action breaks
time-reversal, while the mixed Chern-Simons term with these fields breaks
inversion. We will show that the mixed Chern-Simons term has a well-defined
limit as the gap closes and inversion symmetry is restored, which leads to a
non-trivial response action for the 2D Dirac semimetal.

Alternatively, the response of isolated chiral gapless points in 1D and 3D can
be determined if they are viewed as theories that live on the boundary of a
higher dimensional topological insulator or topological semimetal. In the
presence of gauge fields, the higher dimensional bulk will generate a current
inflow to the boundary to compensate the anomalous response of the gapless
boundary modes. From this perspective, we expect that the effective response
action of Weyl semimetals in odd spatial dimensions can be obtained by taking
the boundary contribution of a higher dimensional system. There are likely
other methods that could be applied to derive these response actions in their
intrinsic spatial dimension, e.g., via the subtle introduction of an auxiliary
$\theta$-field, but we choose our procedure since it reinforces the
dimensional relationships discussed in the previous section and requires fewer
formal tools.

Thus, our strategy for deriving the general form of the coefficients of mixed
crystalline-electromagnetic responses is to begin by deriving effective
response actions in even spatial dimensions, i.e., 2D and 4D. We will do so by
identifying gradient expansion contributions (see Appendix B for a brief
review) that contain an appropriate effective action constructed out of
translational ($\mathfrak{e}^{\lambda}$) and electromagnetic ($A$) gauge
fields. Then the response of semimetals in odd spatial dimensions can be
obtained by looking at the boundary of a response theory defined in one
dimension higher.

###  III.1 Effective responses of 2D semimetals

In this subsection we will derive the coefficients of two 2D response actions
that contain translation gauge fields, namely response action (v) from Fig.
1(c), and the response action in Fig. 2(b). We will find that the coefficients
of these actions are characterized by the dipole and quadrupole moments of the
Berry curvature in the 2D Brillouin Zone, respectively. When we specialize to
2D Dirac semimetals, the distribution of Berry curvature is sharply localized
as $\pm\pi$-fluxes at the Dirac nodes. Hence, the coefficients will become
proportional to the dipole and quadrupole moment of the distribution of Dirac
nodes.

####  III.1.1 Dirac node dipole semimetal

Let us start by considering a gapped $\mathcal{T}$-invariant system having
broken $\mathcal{I}$ symmetry. Under these conditions the electromagnetic
Chern-Simons term, which represents the Hall conductivity, vanishes, and we
can consider the mixed linear response of a momentum current responding to an
electromagnetic field, or vice-versa. Using the Kubo formula, or applying the
gradient expansion procedure described in App. B, we find the contribution to
the effective action (when the chemical potential lies in the insulating gap)
(See also [50]):

$\begin{split}S_{e,A}=-e\int d^{3}r\
\mathfrak{e}^{\alpha}_{\mu}\partial_{\nu}A_{\rho}\int\frac{d\omega
d^{2}k}{(2\pi)^{3}}k_{\alpha}\Omega^{(3)}_{\mu\nu\rho}(\omega,k),\end{split}$
(4)

where

$\Omega^{(3)}_{\mu\nu\rho}(\omega,k)=\text{tr}\left(G_{0}\frac{\partial
G_{0}^{-1}}{\partial k_{\mu}}\frac{\partial G_{0}}{\partial
k_{\nu}}\frac{\partial G_{0}^{-1}}{\partial k_{\rho}}\right),$ (5)

and $G_{0}(k_{\mu})$ is the single-particle Green function. To extract the
coefficient of the $\mathfrak{e}^{\alpha}\wedge dA$ term, we contract
$\Omega^{(3)}_{\mu\nu\rho}$ with the totally antisymmetric tensor
$\tfrac{1}{3!}\varepsilon^{\mu\nu\rho}$. This gives the coefficient

$c_{\alpha}=e\frac{\varepsilon^{\mu\nu\rho}}{3!}\int\frac{d\omega
d^{2}k}{(2\pi)^{3}}k_{\alpha}\ \Omega^{(3)}_{\mu\nu\rho}(\omega,k)$ (6)

of the response action

$S_{e,A}=c_{\alpha}\int\mathfrak{e}^{\alpha}\wedge dA.$ (7)

We note that Eq. 6 is very similar to the response coefficient of the standard
electromagnetic Chern-Simons term apart from the factor of $k_{\alpha}$ in the
integrand. As such, assuming $\alpha=x,y$, we can use a well-established
result to evaluate the frequency integral to obtain [64]:

$\begin{split}\frac{\varepsilon^{\mu\nu\rho}}{24\pi^{2}}\int d\omega d^{2}k\
\Omega^{(3)}_{\mu\nu\rho}(\omega,k)=\frac{1}{2\pi}\int_{BZ}dk_{x}dk_{y}\
\mathcal{F}^{xy}(k_{x},k_{y}),\end{split}$ (8)

where $\mathcal{F}^{xy}$ is the Berry curvature. Hence, we can rewrite
$c_{\alpha}$ as an integral over the BZ by substituting this relationship into
Eq. 6 to find:

$c_{\alpha}=\frac{e}{(2\pi)^{2}}\int_{BZ}dk_{x}dk_{y}\
k_{\alpha}\mathcal{F}^{xy}(k_{x},k_{y}).$ (9)

We have thus arrived at the result that $c_{\alpha}$ is proportional to the
$\alpha$-th component of the dipole moment of the distribution of Berry
curvature. This coefficient can be non-zero since it is allowed by broken
$\mathcal{I}$ and preserved $\mathcal{T}$, i.e.,
$\mathcal{F}^{xy}(\textbf{k})=-\mathcal{F}^{xy}(-\textbf{k}).$ We also note
that $c_{\alpha}$ is independent of the choice of zone center, and shifts of
$k$ in the integrand in general, because the Chern number (Hall conductivity)
vanishes in the presence of $\mathcal{T}$.

In a gapped $\mathcal{T}$-invariant system, restoring $\mathcal{I}$-symmetry
forces $c_{\alpha}$ to vanish, since $\mathcal{F}^{xy}(\textbf{k})=0.$
However, in gapless systems this need not be the case. To see this, we apply
our result from Eq. 9 to a 2D Dirac semimetal by first introducing a weak
perturbation $V_{\mathcal{I}}$ which breaks $\mathcal{I}$ and opens up a small
gap, and then taking the limit $V_{\mathcal{I}}\rightarrow 0$, in which
inversion symmetry is restored. In the gapped system the Berry curvature
$\mathcal{F}^{xy}$ is distributed smoothly across the entire 2D BZ. In the
gapless limit, however, the Berry curvature distribution will develop sharp
peaks of weight $\pi$ localized at the positions of the Dirac points:

$\mathcal{F}^{xy}=\sum_{a=1}^{N_{D}}\pi\chi_{a}\delta(\textbf{k}-\textbf{k}^{a}),$
(10)

where $a$ runs over all Dirac nodes at momenta $\textbf{k}^{a},$ and
$\chi_{a}=\pm 1$ is an integer indicating the sign of the $\pi$-Berry phase
around the Fermi surface of the $a$-th Dirac point at a small chemical
potential above the node [34]. Ultimately, we find the effective response
action of a Dirac node dipole semimetal is given by:

$S_{DD}=\frac{e\mathcal{P}_{\alpha}}{4\pi}\int\mathfrak{e}^{\alpha}\wedge dA,$
(11)

where

$\mathcal{P}_{\alpha}=\sum_{a=1}^{N_{D}}\chi_{a}k^{a}_{\alpha},$ (12)

is the dipole moment of the Dirac nodes.

Note that if the Dirac nodes meet at the zone boundary, they can be
generically gapped even in the presence of $\mathcal{T}\mathcal{I}$ symmetry.
The resulting insulating phase represents a weak TI having
$\mathcal{P}_{\alpha}=G_{\alpha}$, where $G_{\alpha}$ are the components of a
reciprocal lattice vector. In this case, the action in Eq. 11 describes a
stack (i.e., a family of lattice lines/planes corresponding to $G_{a}$) of 1D
polarized TI chains aligned perpendicular to $G_{a}$. To see this explicitly,
take $G_{x}=\frac{2\pi}{a_{x}}$, and set
$e^{\alpha}_{\beta}=\delta^{\alpha}_{\beta}$ in Eq. 11 to obtain the action

$\frac{e}{2}\int\frac{dx}{a_{x}}\left(\int
dydtE_{y}\right)=N_{x}\frac{e}{2}\int dydtE_{y},$ (13)

where $N_{x}$ is the number of unit cells in the $x$-direction. This action is
just $N_{x}$ copies of the usual $\theta$-term action for 1D, electrically-
polarized topological insulators $(\theta=\pi)$ parallel to the $y$-direction,
stacked along $\hat{x}.$

We have now derived Eq. 11 as a quasi-topological contribution to the response
of a 2D Dirac semimetal where the nodes have a dipolar configuration. However,
there is another important subtlety that we will now point out. Earlier work
has shown that the electromagnetic response of 2D Dirac semimetals with both
$\mathcal{T}$ and $\mathcal{I}$ symmetry is an electric polarization
proportional to the Dirac-node dipole moment [34]. Even more recently,
connections have been made between mixed translation-electromagnetic responses
and the electric polarization [52]. Since we have a clear derivation of the
response term we can use our results to understand the precise connection
between the electric polarization and the coefficient $c_{\alpha}$ of the
$e^{\alpha}\wedge dA$ response action. Using the standard approach of Ref. 26,
the polarization in 2D is

$\begin{split}P_{e}^{\alpha}=\frac{e}{(2\pi)^{2}}i\int_{BZ}d^{2}\textbf{k}\langle
u_{\textbf{k}}|\partial_{k_{\alpha}}u_{\textbf{k}}\rangle\end{split}$ (14)

where $\mathcal{A}^{\alpha}(k)=i\langle
u_{\textbf{k}}|\partial_{k_{\alpha}}u_{\textbf{k}}\rangle$ is the Berry
connection. Hence, we find that the electric polarization $P_{e}^{\alpha}$ is
related to $c_{\alpha}$ by an integration by parts (See Appendix C):

$\begin{split}P^{\alpha}_{e}&=\frac{e}{(2\pi^{2})}\varepsilon^{\alpha\beta}\int
d^{2}k\ k_{\beta}\mathcal{F}^{xy}+\frac{e}{2\pi}W^{\alpha}\\\
&=\epsilon^{\alpha\beta}c_{\beta}+\frac{e}{2\pi}W^{\alpha},\end{split}$ (15)

where we have set the lattice constants equal to unity, and the Wilson loop

$W^{\alpha}=\oint dk_{\alpha}\mathcal{A}^{\alpha}(k_{\alpha},k_{\beta}=\pi),$
(16)

is an integral of the Berry connection $\mathcal{A}^{\alpha}$ along the
$\alpha$-th momentum direction at a fixed, inversion-invariant transverse
momentum $k_{\beta}=\pi$ at the boundary of the BZ.

From this explicit relationship we can immediately draw some conclusions.
First, in the Dirac semimetal limit, we reproduce the result of Ref. 34 where
the polarization is proportional to the Dirac node dipole moment:
$P^{\alpha}_{e}=\frac{e}{2(2\pi)}\varepsilon^{\alpha\beta}\mathcal{P}_{\beta}.$
And second, if we have broken inversion symmetry (while $\mathcal{T}$ is still
preserved), we see that the polarization and the coefficient $c_{\alpha}$ are
not quantized, and _not equal_ to each other. This scenario can be found in
inversion-breaking insulators with a Berry curvature dipole moment. These
insulators will have a charge polarization, and they will also have a mixed
translation-electromagnetic response. However, we find from this calculation,
and explicit numerical checks, that they are generically inequivalent.
Ultimately this boils down to the fact that the Wilson loop at the boundary of
the BZ requires a symmetry to be quantized, e.g., mirror or inversion.
Otherwise, the Wilson loop gives a contribution that distinguishes the
polarization and the mixed crystalline-electromagnetic responses. We leave a
detailed discussion of this subtle distinction to future work.

To summarize, Eq. 11 captures the generic mixed crystalline-electromagnetic
response of the bulk of a 2D system with $\mathcal{T}$-symmetry. In the limit
of a Dirac semimetal, the coefficient of the response coincides with the
electric polarization of the system. We note that in this limit there will be
other non-vanishing response terms since the system is gapless, but Eq. 11
represents a distinct contribution to the total response of the system to
electromagnetic and translation gauge fields. We will study an explicit model
with this response term in Sec. IV.1.

####  III.1.2 Dirac node quadrupole semimetal

Now we will move on to discuss the response of quadrupole arrangements of 2D
Dirac nodes as in Fig. 2(b). If the Chern number and momentum dipole moment
$\mathcal{P}_{\alpha}$ vanish, then our semimetal has a well-defined momentum
quadrupole moment, which is independent of the choice of zone center. We now
show that such systems are described by the response action:

$S_{DQ}=\frac{\hbar\mathcal{Q}_{\alpha\beta}}{8\pi}\int\mathfrak{e}^{\alpha}\wedge
d\mathfrak{e}^{\beta}.$ (17)

From the derivation in the previous section we anticipate that, in the limit
of a Dirac semimetal band structure, the coefficient
$\mathcal{Q}_{\alpha\beta}$ of this response action is related to the momentum
quadrupole moment of the Dirac nodes. To confirm this statement let us
consider the linear response of a momentum current to a translation gauge
field for a gapped system. From the Kubo formula, or gradient expansion, we
find a coefficient of the $\mathfrak{e}^{\alpha}\wedge d\mathfrak{e}^{\beta}$
term:

$\mathcal{Q}_{\alpha\beta}\equiv\frac{1}{2}\frac{\varepsilon_{\mu\nu\rho}}{3!}\int\frac{d\omega
d^{2}k}{(2\pi)^{3}}k_{\alpha}k_{\beta}\ \Omega^{(3)}_{\mu\nu\rho}(\omega,k).$
(18)

We can use the relationship mentioned in Eq. 8 to carry out the frequency
integral to obtain the coefficient of Eq. 17:

$\mathcal{Q}_{\alpha\beta}=\frac{1}{\pi}\int_{BZ}dk_{x}dk_{y}\
k_{\alpha}k_{\beta}\mathcal{F}^{xy}(k_{x},k_{y}).$ (19)

To apply this to the Dirac node quadrupole semimetal shown in Fig. 2(b), we
evaluate the response by first introducing a symmetry-breaking mass term and
then studying the topological response of the resulting gapped system. In this
case the mass term breaks $\mathcal{T}$ but has a vanishing total Chern
number. In the example at hand, this can be done by adding a $k$-independent
term that opens a local mass of the same sign for each of the four Dirac
points in Fig. 2b. Such a mass term preserves $\mathcal{I}$, which in the
gapped system automatically guarantees a vanishing dipole moment of the Berry
curvature. This, together with the vanishing Chern number, is necessary so
that the momentum quadrupole moment is well-defined, independent of the choice
of zone center. For this scenario, in the limit that the perturbative mass
goes to zero,

$\mathcal{Q}_{\alpha\beta}=\sum_{a=1}^{N_{D}}\chi_{a}k^{a}_{\alpha}k^{a}_{\beta},$
(20)

which is the Dirac node quadrupole moment. In Sec. IV.2 we will explicitly
study a model with this Berry curvature configuration and a resulting non-
vanishing $\mathcal{Q}_{\alpha\beta}.$ We will see that while the Dirac node
dipole moment captures the electric polarization (see Appendix C), the nodal
quadrupole moment captures a kind of momentum polarization (see Appendix D)
(this time, without the subtlety of the additional Wilson loop contribution
discussed above). For comparison, the surface charge theorem relates the bulk
electric polarization to a boundary charge, and for the momentum polarization
there will be a boundary momentum.

###  III.2 Effective responses of 1D (semi)metals

Figure 3: (a) One-dimensional band structure of an ordinary metal. The pair of
gapless points is marked by the sign of their respective chiralities,
highlighting the momentum-space dipole characterizing the response of the
system. (b) Band structure of a 1D metal characterized by a momentum
quadrupole moment. The system has an integer (vanishing in this case) charge
filling, but a non-zero momentum. (c) Band structure of a 1D metal
characterized by a momentum octupole moment. The system has an integer
(vanishing) filling, a vanishing momentum, but a non-vanishing expectation
value for the square of the momentum. See Appendix E.

Now that we have derived the responses of 2D systems coupled to
electromagnetic and translation gauge fields, we will use Figs. 1(b) and 2 as
guides to generate related responses in 1D and 3D. To get 1D responses we will
consider the boundary response of the 2D systems (this subsection), and we
will stack the 2D responses to get 3D responses of nodal line semimetals (next
subsection). We note that in the following discussion we will treat
translation as a continuous symmetry (as in Appendix A, as this perspective is
useful for obtaining the correct response actions from our diagram
calculations). One can see Ref. 18, for example, for a discussion that treats
the subtleties associated to having a discrete translation symmetry.

It is well-known that chiral modes in 1D are anomalous, i.e., charge is not
conserved when we apply an electric field. In 1D lattice models this anomaly
is resolved because of fermion doubling, i.e., for every right-moving chiral
mode there is a corresponding left-moving mode that compensates the anomaly.
While it is true that the electromagnetic charge anomaly is resolved with such
a lattice dispersion, the doubled system can still be anomalous in a different
but related sense if we have translation symmetry (see Ref. 18 for a similar
discussion).

To be specific, in the presence of translation symmetry we can consider the
momentum current in Eq. 3: $\mathcal{J}_{x}^{\mu}=\hbar k_{x}j^{\mu}$ where
$j^{\mu}$ is the particle number current. At low energies, current-carrying
excitations lie in the vicinity of Fermi points $k^{F,\alpha}_{x}$ and carry
corresponding particle currents $j^{\mu}_{(\alpha)}.$ The total contribution
to momentum current from these low-lying modes is:

$\mathcal{J}^{\mu}_{x}=\sum_{\alpha}\hbar k^{F,\alpha}_{x}j^{\mu}_{(\alpha)}.$
(21)

In the simplest case of a nearest-neighbor lattice model having a single,
partially-filled band, we have two Fermi points: $k^{F}\equiv
k^{F,R}_{x}=-k^{F,L}_{x}$, with $j^{\mu}_{R}=(\rho_{R},v_{F}\rho_{R})$ and
$j^{\mu}_{L}=(\rho_{L},-v_{F}\rho_{L})$, where $\rho$ is the number density.
Interestingly, the momentum current in this scenario is

$\mathcal{J}^{\mu}_{x}=\hbar k^{F}(j^{\mu}_{R}-j^{\mu}_{L}),$ (22)

which, up to a factor of $\hbar k^{F},$ is just the axial current!

Importantly, even though this lattice model does not have an electromagnetic
charge anomaly, $\partial_{\mu}(ej^{\mu}_{L}+ej^{\mu}_{R})=0$, it does have an
axial anomaly:

$\partial_{\mu}(j^{\mu}_{R}-j^{\mu}_{L})=\frac{eE^{x}}{\pi\hbar}.$ (23)

Taking this point of view, we can reformulate the axial anomaly in this system
as a mixed crystalline-electromagnetic anomaly where an electric field $E_{x}$
violates conservation of the $k_{x}$ momentum current,

$\partial_{\mu}\mathcal{J}^{\mu}_{x}=\frac{e\hbar k^{F}}{\pi\hbar}E^{x}.$ (24)

More generally the anomaly is proportional to the momentum dipole moment of
the Fermi points, which replaces a factor of $2k_{F}$ in Eq. 24 (see App. E).

There is a conjugate effect that occurs in an applied strain field, which can
be implemented as a translation electric field
$\mathcal{E}^{x}_{x}=\partial_{x}\mathfrak{e}_{0}^{x}-\partial_{t}\mathfrak{e}_{x}^{x}.$
Naively such a non-vanishing field will generate violations to the
conservation law for the usual electromagnetic current according to

$\partial_{\mu}(ej^{\mu})=\frac{ek^{F}}{\pi}\mathcal{E}^{x}_{x},$ (25)

(again see App. E for a more general expression in terms of the momentum
dipole). However, this equation is not quite correct if we have an isolated
system with a fixed number of electrons, and hence, we must be careful when
considering time-dependent changes to $\mathfrak{e}_{x}^{x}$ as we will now
describe.

To gain some intuition for Eq. 25, consider increasing the system size by one
lattice constant $a$ during a time $T$ by adding an extra site to the system:
$\int dxdt\ \mathcal{E}_{x}^{x}=a$ (one can also think of threading a
dislocation into the hole of a 1D periodic system). From the anomaly equation
we would find that the amount of charge in the system changes by
$ek^{F}a/\pi,$ as one would expect for adding a unit cell to a translationally
invariant system having a uniform charge density $\rho=ek_{F}/\pi.$ However,
there is a subtlety that we can illustrate by considering a system having a
fixed number of electrons $N_{e}=k^{F}L_{x}/\pi,$ which we strain by uniformly
increasing the lattice constant. Assuming a uniform system, the anomalous
conservation law in this case becomes

$\partial_{t}\rho=\partial_{t}\left(\frac{ek^{F}}{\pi}\mathfrak{e}_{x}^{x}\right).$
(26)

Crucially, we note that if we increase the system size with fixed particle
number, then $k^{F}$ will decrease. Indeed, in the small deformation limit the
momenta are proportional to $(\mathfrak{e}^{x}_{x})^{-1}$ since their finite
size quantization depends inversely on the system size. Using this result, the
conservation law becomes:

$\partial_{t}\rho=\frac{e}{\pi}\left(\mathfrak{e}_{x}^{x}\partial_{t}k^{F}+k^{F}\partial_{t}\mathfrak{e}^{x}_{x}\right)=\frac{ek^{F}}{\pi}(-\partial_{t}\mathfrak{e}_{x}^{x}+\partial_{t}\mathfrak{e}_{x}^{x})=0$
(27)

where we used
$\partial_{t}(\mathfrak{e}_{x}^{x})^{-1}=-(\mathfrak{e}_{x}^{x})^{-2}\partial_{t}\mathfrak{e}_{x}^{x}.$

The outcome that $\partial_{t}\rho=0$ is the result one would expect by
stretching the system uniformly while keeping the number of particles fixed.
To clarify, at a fixed particle number we know the total charge cannot change,
however it perhaps seems counter-intuitive that the _density_ does not
decrease if we stretch the system. The reason is that the quantity $\rho$
above, which is defined as $\tfrac{\delta S}{\delta A_{0}},$ is not a scalar
density. Indeed, for general geometries the scalar charge density would be
defined as

$\bar{\rho}=\frac{1}{\mathfrak{e}^{x}_{x}}\frac{\delta S}{\delta A_{0}},$ (28)

where the $\mathfrak{e}_{x}^{x}$ is essentially playing the role of the
determinant of a spatial metric. To calculate the total charge we would then
use

$Q=\int dx\,\mathfrak{e}_{x}^{x}\bar{\rho}=\int dx\,\rho.$ (29)

Indeed, the scalar charge density $\bar{\rho}$ will decrease as the system is
stretched since $\partial_{t}\bar{\rho}\propto\partial_{t}\mathcal{P}_{x}$
which decreases as the system size increases at fixed electron number.

The effective response action of the 1D system can be derived as a boundary
effective action of an appropriate 2D theory. In fact, we have already seen
such a 2D system when studying the 2D Dirac semimetal with Dirac nodes
arranged in a dipolar fashion. The bulk response for this 2D system with a
weak inversion-breaking gap is Eq. 11. As mentioned above, this bulk theory
implies that the system has an electric polarization. From the surface-charge
theorem for polarization we expect that the boundary will have a charge
density equal the polarization component normal to the boundary. The
contribution to the buondary effective action from Eq. 11 is:

$S_{\partial}=\frac{e}{4\pi}\mathcal{P}_{\alpha}\int\mathfrak{e}^{\alpha}\wedge
A.$

From this we can extract the boundary charge density:
$\rho_{\partial}=\frac{e}{2}\frac{\mathcal{P}_{\partial}}{2\pi}\mathfrak{e}_{\partial}^{\partial}$
where $\mathcal{P}_{\partial}$ is the component along the boundary, and
$\mathfrak{e}_{\partial}^{\partial}$ is the diagonal translation gauge field
component along the boundary that is simply equal to unity in non-deformed
geometries.

While the form of this action is what we expect for a 1D metal, the
coefficient is half the size it should be. The reason is that on the edge of
the 2D Dirac semimetal, the momentum-space projections of the bulk Dirac nodes
in the edge BZ represent points where the edge-filling changes by $\pm e/2,$
[34] not $\pm e$ as would be the case for a 1D Fermi-point in a metal. Hence
for a metal we expect a result twice as large (we will see a similar result in
Sec. III.5 when comparing the boundary response of a 4D system to that of a 3D
Weyl semimetal). Thus the action for the 1D system is

$S_{1D,D}=\frac{e}{2\pi}\mathcal{P}_{\alpha}\int\mathfrak{e}^{\alpha}\wedge
A.$ (30)

From this form we can identify $\mathcal{P}_{\alpha}=(-\Delta\mu/\hbar,\Delta
k_{x})$ such that $\frac{\mathcal{P}_{x}}{2\pi}$ is simply the filling
fraction of the 1D metal and $\frac{\mathcal{P}_{t}}{2\pi}$ measures the
imbalance of left- and right-moving excitations in the system
($\Delta\mu=\mu_{R}-\mu_{L}$).

Introducing a charge current vector

$j^{\mu}=\frac{e}{2\pi}\varepsilon^{\mu\nu}\mathcal{P}_{\nu}=\frac{e}{2\pi}\left(\Delta
k_{x},\Delta\mu/\hbar\right)^{T}$ (31)

we can recast Eq. 31 in the most familiar form: $S_{1D,D}=\int dtdx\
j^{\mu}A_{\mu}$. Thus, we have now generated the action (ii) from Fig. 1(c).
Let us also note that the edge states of the Dirac semimetal can be flat,
while the 1D context we mentioned above has a dispersion. However, the key
feature of both cases is that as momentum is swept across the 1D BZ (1D
surface BZ for the 2D case) the filling of the states changes in discrete
jumps at either the Fermi points in 1D, or the (surface-projected) Dirac
points in 2D. It is this change in the filling that is captured by the
quantity $\mathcal{P}_{x},$ and does not depend on the dispersion in a crucial
way.

Now that we have this example in mind, we can ask what the analogous 1D
boundary system is for the Berry curvature quadrupole action Eq. 17. We
mentioned that this bulk response represents a momentum polarization, which
implies that the boundary should have a momentum density parallel to the edge.
Indeed, we expect that such a 1D system will have a vanishing Fermi-point
dipole moment (i.e., the filling is integer), but a quadrupole moment that is
non-vanishing (see Fig. 3(b)).

From the point of view of the translation gauge fields, such band structures
are chiral since either the right movers or left movers carry larger momentum
charge. To see this, consider a 1D Fermi surface with right-movers at momenta
$\pm K_{F}$, and left-movers at momenta $\pm Q_{F}$. Let us further restrict
our attention to currents for which the net number of right-movers (and of
left-movers) is zero, e.g. $\rho_{R}(K_{F})+\rho_{R}(-K_{F})=0$. Defining
$\delta\rho_{R}=(\rho_{R}(K_{F})-\rho_{R}(-K_{F}))$, and
$\delta\rho_{L}=(\rho_{L}(Q_{F})-\rho_{L}(-Q_{F}))$, we see that the momentum
gauge field couples to

$\mathcal{J}^{\mu}_{x}=K_{F}\delta\rho_{R}+Q_{F}\delta\rho_{L}\ .$ (32)

Thus we see that for $K_{F}\neq Q_{F}$ (as in Fig. 3(b)), the momentum gauge
field couples differently to right- and left- moving density fluctuations. In
the extreme limit that $Q_{F}=0$, the momentum gauge theory is fully chiral.

More generally, in a 1D system with a Fermi-point quadrupole (c.f. Eq. 20)
$\mathcal{Q}_{xx}=\sum_{a=1}^{N_{F}}={\rm{sgn}}(v_{Fa})(k_{x}^{(a)})^{2}$, and
fixed electric charge, this chiral coupling leads to an anomaly in the
presence of a non-vanishing translation gauge field:

$\partial_{\mu}\mathcal{J}^{\mu}_{x}=\frac{\hbar\mathcal{Q}_{xx}}{4\pi}\mathcal{E}_{x}^{x}.$
(33)

This anomaly implies that if we turn on a translation gauge field (e.g., via
strain) then we will generate momentum as shown in App. E 222As shown in the
Appendix, this anomaly has two contributions. One comes from the low-energy
currents that contribute with a factor of $1/2\pi$ and a second from a change
of system-size for a ground state carrying a non-vanishing momentum density
with a factor of $-1/4\pi$..

The response theory describing such a 1D system is similar to that describing
the chiral boundary of a Chern-Simons theory. Indeed, if we start from Eq. 17
and derive the boundary response (and compensate for a similar factor of two
as mentioned above in the momentum-dipole case) we arrive at an effective
action:

$S=-\frac{\hbar}{4\pi}\int dtdx\
\left(\mathcal{Q}_{xx}\mathfrak{e}^{x}_{x}\mathfrak{e}^{x}_{t}+\mathcal{Q}_{xt}\mathfrak{e}^{x}_{x}\mathfrak{e}^{t}_{t}\right).$
(34)

In this effective action the momentum quadrupole moment of the Fermi points
$\mathcal{Q}_{xx}$ encodes the ground state momentum density (see Appendix E).
The quantity $\mathcal{Q}_{xt}$ is the mixed Fermi-point quadrupole moment in
momentum and energy, but we leave a detailed discussion of such mixed moments
to future work.

The arguments of this section can be extended to look at higher moments of the
chirality-weighted Fermi momenta, which are proportional to the ground state
expectation values of higher and higher powers of momenta. To describe these
properties, and related response phenomena, we can introduce gauge fields
$\mathfrak{e}^{abc\ldots}$ that couple to higher monomials of momentum,
$k_{a}k_{b}k_{c}\ldots.$ For example, the fields that couple to zero powers or
one power of momentum are the electromagnetic $A$ and translation gauge fields
$k_{x}\mathfrak{e}^{x}$ respectively, and we could introduce a coupling
$k_{a}k_{b}\mathfrak{e}^{ab}$ to the set of 1-form gauge fields
$\mathfrak{e}^{ab},$ e.g., $k_{x}^{2}\mathfrak{e}^{xx}$. We describe the
hierarchical anomalies associated to these gauge fields in Appendix E.

###  III.3 Effective responses of 3D nodal line semimetals

We can now use our 2D results from Sec. III.1 to generate the responses of two
types of nodal line semimetals in 3D. To generate the two types we imagine
stacking either the action in Eq. 11 or the action in Eq. 17. The action
resulting from the former has been discussed in Refs. [35, 19]; the second is,
to the best of our knowledge, new. From our arguments for gapped systems in
Sec. II, we expect that the form of the actions we obtain from stacking will
contain an extra wedge product with the translation gauge field in the
stacking direction. To be explicit, suppose we are stacking up 2D semimetals
(that are parallel to the $xy$-plane) into the $z$-direction. By stacking
decoupled planes of the responses in either Eq. 11 or Eq. 17, we expect to
find

$S=\frac{e\mathcal{P}_{\alpha}}{4\pi
a_{z}}\int\mathfrak{e}^{z}\wedge\mathfrak{e}^{\alpha}\wedge dA,$

or

$S=\frac{\hbar\mathcal{Q}_{\alpha\beta}}{8\pi
a_{z}}\int\mathfrak{e}^{z}\wedge\mathfrak{e}^{\alpha}\wedge
d\mathfrak{e}^{\beta},$

respectively, where $\alpha,\beta=x,y.$ The forms of these actions match
action (viii) in Fig. 1(c) and the action in Fig. 2(c) respectively. We note
that the stacked, decoupled systems simply inherit the response coefficient of
the 2D system.

We want to consider more general configurations of systems with stacked and
coupled planes, perhaps stacked in several directions. As we have seen, if the
layers we stack are decoupled, then each layer contributes the same amount.
This contribution (for a stack in the $z$-direction) is captured by the
integral $\tfrac{1}{a_{z}}\int\mathfrak{e}^{z}=N_{z}$ where $N_{z}$ is the
number of layers. However, if the layers are coupled, then each fixed-$k_{z}$
plane can have a different amount of Dirac node dipole
($\mathcal{P}_{\alpha}(k_{z})$) or Dirac node quadrupole moment
($\mathcal{Q}_{\alpha\beta}(k_{z})$) respectively. The total coefficient is
then determined by the sum over all values of $k_{z}.$ One can also have
stacks in any direction, not just the $z$-direction. Hence, in this more
generic scenario the actions become

$S_{DD3}=e\mathcal{B}_{\alpha\beta}\int\mathfrak{e}^{\alpha}\wedge\mathfrak{e}^{\beta}\wedge
dA,$ (35)

and

$S_{DQ3}=\hbar\mathcal{B}_{\alpha\beta,\gamma}\int\mathfrak{e}^{\alpha}\wedge\mathfrak{e}^{\beta}\wedge
d\mathfrak{e}^{\gamma},$ (36)

with coefficients

$\mathcal{B}_{\alpha\beta}=\frac{1}{4(2\pi)^{3}}\epsilon^{\alpha\beta\sigma}\int
d^{3}k\ k_{\delta}\mathcal{F}^{\sigma\delta}$ (37)

and

$\mathcal{B}_{\alpha\beta,\gamma}=\frac{1}{6(2\pi)^{3}}\epsilon^{\alpha\beta\sigma}\int
d^{3}k\ k_{\gamma}k_{\delta}\mathcal{F}^{\sigma\delta}.$ (38)

where $\mathcal{F}^{\mu\nu}$ is the Berry curvature of the
$k_{\mu}k_{\nu}$-plane. These forms of the coefficients capture scenarios with
more complicated nodal line geometries. Indeed, as previously shown in Ref. 35
the coefficient $\mathcal{B}_{\alpha\beta}$ is determined by the line nodes
that have non-vanishing area when projected into the $\alpha\beta$-plane.
Additionally, for nodal line semimetals with $\mathcal{T}\mathcal{I}$ symmetry
the coefficient is proportional to the charge polarization in the direction
normal to the $\alpha\beta$-plane [35]. We can see this explicitly by
integrating Eq. 37 by parts with the same caveats mentioned in Sec. III.1.1
surrounding Eq. 15.

Analogously, the coefficient $\mathcal{B}_{\alpha\beta,\gamma}$ can represent
a kind of “momentum”-polarization where the polarization is again normal to
the $\alpha\beta$-plane, and the charge that is polarized is the momentum
along the $\gamma$-direction. We can see this heuristically by integrating by
parts using the derivatives in the $\mathcal{F}^{\sigma\delta}$ to find

$\mathcal{B}_{\alpha\beta,\gamma}\sim-\frac{1}{2(2\pi)^{3}}\int
d^{3}k\left(\epsilon^{\alpha\beta\sigma}k_{\gamma}\mathcal{A}^{\sigma}-\epsilon^{\alpha\beta\gamma}k_{i}\mathcal{A}^{i}\right)$
(39)

where we have used the $\sim$ symbol to indicate that there are boundary terms
we have dropped that can be important if the line nodes span the Brillouin
zone. We can see from this form that the coefficient for the case when
$\alpha,\beta,\gamma$ are not all different, e.g. $\mathcal{B}_{xz,x}$, is
proportional to the polarization in the $y$-direction (i.e. normal to the
$xz$-plane) weighted by the momentum in the $x$-direction.

We note that for $\mathcal{B}_{\alpha\beta}$ to be well-defined, the Chern
number in each plane must vanish. In addition to this constraint,
$\mathcal{B}_{\alpha\beta}=0$ is a necessary constraint for
$\mathcal{B}_{\alpha\beta,\gamma}$ to be well defined. These hierarchical
requirements are analogous to the usual requirements for the ordinary
(magnetic) dipole and (magnetic) quadrupole moments of the electromagnetic
field to be independent of the choice of origin. Here the role of the magnetic
field distribution is being played by $\mathcal{F}^{\sigma\rho}(k)$, and, for
example, the constraint on the vanishing Chern number eliminates the
possibility of magnetic monopoles (i.e., Weyl points).

###  III.4 Effective responses of 4D semimetals

Our next goal is to determine the coefficients for the response actions of 3D
Weyl point-node semimetals. However, because the Weyl nodes in 3D exhibit an
anomaly, the responses are subtle to calculate intrinsically in 3D. Instead,
to accomplish our goal we will first carry out more straightforward
calculations of the responses of 4D semimetals and then return to 3D either by
considering the boundary of a 4D system, or by compactifying and shrinking one
dimension of the bulk. Hence, as a step toward 3D semimetals, in this
subsection we provide the derivation for effective response actions of
semimetals in 4D.

The first action we consider is of the form

$S=c_{\alpha}\int\mathfrak{e}^{\alpha}\wedge dA\wedge dA,$ (40)

where for our purposes $\alpha=x,y,z,w.$ Collecting all terms in the gradient
expansion that have this field content we obtain:

$\begin{split}S&=\frac{e^{2}}{\hbar}\int d^{5}r\
\mathfrak{e}^{\alpha}_{\mu}\partial_{\nu}A_{\rho}\partial_{\sigma}A_{\tau}\\\
&\times\int\frac{d\omega
d^{4}k}{(2\pi)^{5}}k_{\alpha}\Omega^{(5)}_{\mu\nu\rho\sigma\tau}(\omega,k),\end{split}$
(41)

where

$\Omega^{(5)}_{\mu\nu\rho\sigma\tau}(\omega,k)=\text{tr}\left(G_{0}\frac{\partial
G_{0}^{-1}}{\partial k_{\mu}}\frac{\partial G_{0}}{\partial
k_{\nu}}\frac{\partial G_{0}^{-1}}{\partial k_{\rho}}\frac{\partial
G_{0}}{\partial k_{\sigma}}\frac{\partial G_{0}^{-1}}{\partial
k_{\tau}}\right),$ (42)

and $G_{0}(\omega,k)$ is the single-particle Green function. To determine the
coefficient $c_{\alpha}$ we project this coefficient onto the totally
antisymmetric part and then, just as in Eq. 8, we can carry out the frequency
integral [64] to obtain the simpler expression

$\begin{split}&\int\frac{d\omega
d^{4}k}{2\pi}\frac{\varepsilon_{\mu\nu\rho\sigma\tau}}{5!}k_{\alpha}\Omega^{(5)}_{\mu\nu\rho\sigma\tau}(\omega,k)\\\
&=\frac{1}{16}\int_{BZ}d^{4}\textbf{k}\
k_{\alpha}\varepsilon_{ijkl}\mathcal{F}^{ij}\mathcal{F}^{kl}.\end{split}$ (43)

Hence, the response coefficient takes the form

$\begin{split}c_{\alpha}=\frac{e^{2}}{\hbar}\frac{1}{16(2\pi)^{4}}\int_{BZ}d^{4}k\
k_{\alpha}\varepsilon_{ijkl}\mathcal{F}^{ij}\mathcal{F}^{kl}=\frac{e^{2}\mathcal{P}_{\alpha}}{16\pi^{2}\hbar},\end{split}$
(44)

where we introduced

$\mathcal{P}_{\alpha}=\frac{1}{16\pi^{2}}\int_{BZ}d^{4}\textbf{k}\
k_{\alpha}\varepsilon_{ijkl}\mathcal{F}^{ij}\mathcal{F}^{kl}.$ (45)

As we see from this calculation, similar to 2D, the 4D response theories can
be characterized by the distribution of the quantity
$\varepsilon_{ijkl}\mathcal{F}^{ij}\mathcal{F}^{kl}$ across the 4D Brillouin
zone. For our focus, let us consider the case where the 4D system is a
semimetal with a set of isolated Dirac points (linearly dispersing band
touchings where four bands meet). Without symmetry, these Dirac points are
locally unstable in momentum space to the opening of a gap. If we open up an
infinitesimally small energy gap, the quantity
$\varepsilon_{ijkl}\mathcal{F}^{ij}\mathcal{F}^{kl}$ becomes well-defined
across the entire BZ and its distribution takes the following form in the
massless limit:

$\varepsilon_{ijkl}\mathcal{F}^{ij}\mathcal{F}^{kl}=\sum_{a=1}^{N_{D}}16\pi^{2}\chi_{a}\delta(\textbf{k}-\textbf{k}_{a}).$
(46)

If we substitute this into Eq. 45 then we immediately see that
$\mathcal{P}_{\alpha}$ becomes the momentum space dipole of the set of 4D
Dirac nodes. Let us also comment that if we integrate Eq. 45 by parts we see
that $\mathcal{P}_{\alpha}$ can also be interpreted as a set of magneto-
electric polarizabilities [5, 28]. Just as in the case of the polarization of
a 2D Dirac semimetal, the integration by parts will generate a boundary term
that captures the magneto-electric polarizability coming from the 3D
boundaries of the 4D BZ. Hence, the connection between the total magneto-
electric polarizability and the mixed translation-electromagnetic response is
only exact in the symmetric limit when the boundary term is quantized.

In summary, a 4D response of a system characterized by a dipolar distribution
of the $\varepsilon_{ijkl}\mathcal{F}^{ij}\mathcal{F}^{kl}$ quantity reads:

$S=\frac{e^{2}\mathcal{P}_{\alpha}}{16\pi^{2}\hbar}\int\mathfrak{e}^{\alpha}\wedge
dA\wedge dA.$ (47)

Similar to 2D, if the dipolar response vanishes we can obtain a momentum
quadrupole response coefficient for the action:

$S=\frac{e\mathcal{Q}_{\alpha\beta}}{16\pi^{2}}\int\mathfrak{e}^{\alpha}\wedge
d\mathfrak{e}^{\beta}\wedge dA,$ (48)

where $\mathcal{Q}_{\alpha\beta}$ is a symmetric matrix determined by the
momentum space quadrupole moment of the 4D Dirac nodes. Finally, if both the
dipolar and quadrupolar responses vanish we can consider an octupolar
distribution that will give the response coefficient for the action:

$S=\frac{\hbar\mathcal{O}_{\alpha\beta\gamma}}{48\pi^{2}}\int\mathfrak{e}^{\alpha}\wedge
d\mathfrak{e}^{\beta}\wedge d\mathfrak{e}^{\gamma},$ (49)

where $\mathcal{O}_{\alpha\beta\gamma}$ is determined by the momentum space
octupole moment of the 4D Dirac nodes. We will leave the discussion of
octupolar configurations of Dirac and Weyl nodes to future work. We also
mention that, similar to 2D, for these responses to be independent of the
choice of BZ origin we require that the second Chern number of the 4D system
vanishes. Alternatively, if the second Chern number is non-vanishing, then the
boundary of the system will contain a non-vanishing chirality of Weyl nodes.
As such, the anomalous charge response of the chiral boundary will not allow
us to uniquely determine the momentum response on the boundary.

Before moving on to 3D, let us briefly present some physical intuition about
the response in Eq. 47. Consider a 4D time-reversal and inversion invariant
system having two Dirac nodes separated in the $k_{z}$-direction. To simplify
the discussion, let us also assume the system has mirror symmetry $M_{z}.$ The
assumed symmetries imply that each fixed-$k_{z}$ volume can be treated as an
independent 3D insulator having 3D inversion symmetry, and hence the magneto-
electric polarizability of these 3D insulator subspaces is quantized [5, 66,
67]. Now, if we sweep through $k_{z}$ then each bulk 4D Dirac point crossing
changes the magneto-electric polarizability of the fixed-$k_{z}$ volume by a
half-integer (i.e., changes the related axion angle by $\pi$) [5]. Since the
magneto-electric polarizability jumps between its quantized values as we pass
through the two bulk Dirac nodes, the $k_{z}$ Brillouin zone splits into two
intervals: (i) an interval with a vanishing magneto-electric polarizability,
and (ii) an interval with a non-vanishing quantized magneto-electric
polarizability. Indeed, we could have anticipated this result from the form of
the action Eq. 47 when $\alpha=z,$ i.e., the action represents stacks of 3D
topological insulators that each have a non-vanishing magneto-electric
polarizability.

###  III.5 Effective responses of 3D semimetals

From this discussion we see that, in the presence of symmetry, the 4D bulk
Dirac node dipole moment determines the magneto-electric polarizability of
these 4D topological semimetals via Eq. 47. We want to connect this result to
3D semimetals in two ways. First, we will consider the 3D boundary of the 4D
system, and then we will consider the spatial compactification of one spatial
dimension.

Let us begin by considering the boundary response action from Eq. 47. For the
model system described at the end of the previous subsection we know the
system has a $k_{z}$-dependent magneto-electric polarizability. Consider a
boundary in the fourth spatial direction $w.$ Since the magneto-electric
polarizability is changing from inside to outside of the boundary, the
boundary itself will have a non-vanishing Hall conductivity. For our example
system, each fixed-$k_{z}$ slice of this boundary will have a Hall
conductivity $\sigma_{xy},$ which is quantized, but possibly vanishing.
Additionally, since the bulk 4D Dirac nodes are separated in the $k_{z}$
direction, they will project to gapless points in the 3D surface BZ (on
surfaces that have at least one direction perpendicular to the $z$-direction)
where the Hall conductivity discretely jumps by
$\Delta\sigma_{xy}=\pm\tfrac{e^{2}}{2h}.$

From this phenomenology, i.e., discrete Hall conductivity jumps as we sweep
through $k_{z}$ we expect that the boundary response of Eq. 47 captures the
same response as a Weyl semimetal that has a non-vanishing momentum space
dipole moment of the Weyl nodes in the $z$-direction. Indeed the generic
boundary contribution from Eq. 47 has the form:

$S_{WD}=\frac{e^{2}\mathcal{P}_{\alpha}}{8\pi^{2}\hbar}\int\mathfrak{e}^{\alpha}\wedge
dA\wedge A$ (50)

which was proposed by Ref. 33 to describe the response of Weyl semimetals,
though in the more conventional form using an axion field and without the
translation gauge field. Here $\mathcal{P}_{\alpha},$ $\alpha=x,y,z$ is the
momentum dipole of the Weyl nodes in the $\alpha$-th direction. This action is
represented as (ix) in Fig. 1(c). We note that the coefficient in Eq. 50 is
twice as large as the actual boundary term derived from Eq. 47. This is
because when $k_{i}$ passes through a single Weyl point we have
$\epsilon_{ijk}\Delta\sigma_{jk}=\pm\tfrac{e^{2}}{h},$ where the surface the
response of the 4D system has jumps of half the size. This is analogous to the
fact that a 1D metal has an integer jump in the filling as we pass through a
Fermi point, whereas the surface of a 2D Dirac semimetal has a boundary
“filling” that jumps by a half-integer as we pass through a gapless point in
the surface BZ.

We can repeat this analysis for Eq. 48. The coefficient of this term is
proportional to the momentum space quadrupole moment of the nodal points.
Unfortunately the phenomenology of this term is not as easy to analyze in 4D
because it is not generated from a lower dimensional system in a clear way
333Even though there is a $\mathfrak{e}^{\alpha}$ wedge product with a lower-
dimensional action, it is not transverse to the lower-dimensional action since
$\mathcal{Q}_{\alpha\beta}$ is symmetric. For example, there will be terms
where, say, $\mathfrak{e}^{x}$ couples to $d\mathfrak{e}^{x}$, which cannot be
interpreted as a conventional stacked action.. By analogy with the previous
case, the bulk 4D Dirac nodes will project to a quadrupole of 3D Weyl nodes on
the surface. We can extract the form of the 3D action we want by taking the
boundary term generated from Eq. 48. Then accounting for the factor of two as
in the previous case, we arrive at:

$S_{WQ}=\frac{e\mathcal{Q}_{\alpha\beta}}{8\pi^{2}}\int\mathfrak{e}^{\alpha}\wedge
d\mathfrak{e}^{\beta}\wedge A.$ (51)

(Note that since $\mathcal{Q}_{\alpha\beta}$ is symmetric, the related
contribution of the form
$e\mathcal{Q}_{\alpha\beta}/8\pi^{2}\int\mathfrak{e}^{\alpha}\wedge\mathfrak{e}^{\beta}\wedge
dA$ vanishes). This action is the same as that shown in Fig. 2(a). It produces
a mixed crystalline-electromagnetic response and represents a rank-2 vector
charge response when certain mirror symmetries are preserved [17]. Its
response coefficient is determined by the momentum space quadrupole moment of
the Weyl nodes.

Finally, we come to the action (x) in Fig. 1(c). Let us briefly sketch some
salient features of this response, while we leave a detailed discussion to
future work. We can arrive at this action using a formal compactification of
the action in Eq. 47 [5]. First we can integrate that action by parts to
arrive at

$\frac{e^{2}\mathcal{P}_{\alpha}}{16\pi^{2}\hbar}\int A\wedge
d\mathfrak{e}^{\alpha}\wedge dA,$

where we have ignored the boundary term. We now want to dimensionally reduce
the fourth spatial direction $w$, which we accomplish by choosing periodic
boundary conditions in $w$ and letting the size of the system in this
direction shrink toward zero. In this limit any derivatives with respect to
$w$ are (formally in our case) dropped 444Alternatively we can assume the
fields $A_{w},e^{\alpha}_{w}$ are locked to their ground state values and thus
have vanishing derivatives in all directions.. The resulting non-vanishing
contribution is

$\frac{e^{2}\mathcal{P}_{\alpha}}{8\pi^{2}\hbar}\oint A_{w}dw\int
d\mathfrak{e}^{\alpha}\wedge dA,$

where the integral and exterior derivative in the second factor are over only
the remaining four spacetime dimensions. We can now make the definition

$\Theta\equiv 2\pi\frac{e}{h}\int A_{w}dw,$ (52)

to arrive at action (x) from Fig. 1(c):

$\frac{e\mathcal{P}_{\alpha}}{8\pi^{2}}\int\Theta d\mathfrak{e}^{\alpha}\wedge
dA.$ (53)

To illustrate some of the phenomenology of this action let us assume that
$\mathcal{P}_{z}\neq 0.$ Additionally let us assume that we maintain time-
reversal and inversion symmetry. As such, $\Theta=0,\pi$. To begin, we see
that the action in Eq. 53 is a total derivative if $\Theta$ and
$\mathcal{P}_{\alpha}$ are space-time independent. The resulting pure boundary
term is just proportional to the response of a 2D weak TI (or 2D Dirac
semimetal), i.e., Eq. 7. Depending on the symmetry of the surfaces, this
implies that we expect the surface to be gapped except for possibly isolated
Dirac points. Since the boundary terms appear as $\mathfrak{e}^{z}\wedge dA$
we expect that surfaces normal to $\hat{x}$ ($\hat{y}$) will harbor a
$y$-polarization ($x$-polarization), i.e., the polarization is tangent to the
surface.

Importantly, the sign of the polarization depends on the interpolation of
$\Theta$ between its non-trivial bulk value of $\Theta=\pi$ and the trivial
vacuum value $\Theta=0$ outside the system. For neighboring surfaces where the
effective sign of the polarization changes we anticipate hinge charges where
surfaces intersect since the polarizations are converging or diverging from
the hinges. Thus, the response of this system is similar to a stack of 2D
planes of quadrupole moment having component $q_{xy}\neq 0.$ In this scenario,
coupled quadrupole planes could lead to either a higher order weak topological
insulator having a quadrupole moment, or a higher order topological semimetal
with boundary (and possibly bulk) Dirac nodes [70, 71]. To make further
progress it would be advantageous to have a microscopic derivation of the
coefficient in Eq. 53 intrinsically in 3D. Hence, we will leave further
discussion of this action to future work.

##  IV Explicit Response Calculations for Lattice Models

Now that we have completed the derivations of the actions in Figs. 1(c) and 2,
we will provide a series of model examples that manifest these responses.
Using these models we can numerically calculate the various charge and
momentum responses to electromagnetic and translation gauge fields, providing
an independent verification of the coefficients derived in the previous
section. Some of the models and responses we discuss below have appeared
elsewhere in the literature, while others are have not. We will carry out this
analysis in the same order as the previous section, i.e., point-node Dirac
semimetals in 2D, nodal line semimetals in 3D, and then point-node Weyl
semimetals in 3D. Calculations for 1D systems were carried out analytically in
Sec. III.2, and additional discussion can be found in App. E.

###  IV.1 2D Dirac node dipole semimetal and insulator

We begin with the time-reversal invariant 2D systems discussed in Sec. III.1
that exhibit a mixed crystalline-electromagnetic response. Since $\mathcal{T}$
is preserved, the usual Chern-Simons, Hall-effect response of the
electromagnetic field vanishes. Instead, the response action derived in the
Sec. IV.1 takes the form of a mutual Chern-Simons term [52]:

$S[\mathfrak{e}^{\lambda}_{\nu},A_{\mu}]=\frac{e}{4\pi}\mathcal{P}_{\lambda}\int\mathfrak{e}^{\lambda}\wedge
dA.$ (54)

Unlike the purely electromagnetic polarization response action considered in
Ref. 34, this formulation of the low-energy response theory also includes bulk
electromagnetic responses to the translation gauge fields. For example, by
taking a functional derivative with respect to $A_{\mu}$ we have

$\begin{split}&\rho=-\frac{e}{4\pi}\mathcal{P}_{\lambda}\varepsilon^{ij}\partial_{i}\mathfrak{e}^{\lambda}_{j},\\\
&j^{x}=\frac{e}{4\pi}\mathcal{P}_{\lambda}(\partial_{t}\mathfrak{e}^{\lambda}_{y}-\partial_{y}\mathfrak{e}^{\lambda}_{t}),\\\
&j^{y}=-\frac{e}{4\pi}\mathcal{P}_{\lambda}(\partial_{t}\mathfrak{e}^{\lambda}_{x}-\partial_{x}\mathfrak{e}^{\lambda}_{t}).\end{split}$
(55)

We see that the first equation predicts an electric charge density localized
on a dislocation in the bulk of the lattice, which is exactly the
phenomenology we expect for a weak topological insulator [38] or a 2D Dirac
semimetal. The action (54) also predicts a bulk momentum response to the
electromagnetic field when varied with respect to $\mathfrak{e}^{\mu}$,

$\begin{split}&\mathcal{J}^{t}_{\lambda}=-\frac{e}{4\pi}\mathcal{P}_{\lambda}B_{z},\\\
&\mathcal{J}^{i}_{\lambda}=-\frac{e}{4\pi}\mathcal{P}_{\lambda}\varepsilon^{ij}E_{j},\end{split}$
(56)

where $E_{i}$ and $B_{i}$ are the components of electric and magnetic fields
respectively. In the inversion-symmetric limit and in the absence of lattice
defects and deformations, for which the crystalline gauge fields reduce to
$\mathfrak{e}^{\lambda}_{\mu}=\delta^{\lambda}_{\mu},$ Eq. (55) simply
reproduces the boundary charge and current responses of an ordinary 2D Dirac
semimetal or weak topological insulator, which harbors a non-vanishing
electric polarization. However, as we mentioned in Sec. III.1.1, and comment
further on below, we do not expect the coefficient of this action to match the
electric polarization when inversion is strongly broken.

While the electric polarization/magnetization responses of Dirac semimetals
were discussed in detail in Ref. 34, the momentum responses in Eq. 56, and the
charge responses to translation fluxes (i.e., dislocations) in Eq. 55 are less
familiar. Thus, we will explicitly calculate these responses using a minimal
tight-binding model. For simplicity, we employ a two-band Bloch Hamiltonian
that can model both 2D Dirac semimetals and weak topological insulators:

$\begin{split}H(\textbf{k})&=V_{\mathcal{I}}\sigma^{x}+\sin(k_{y}a_{y})\sigma^{y}\\\
&+(m-\cos(k_{x}a_{x})-\cos(k_{y}a_{y}))\sigma^{z}.\end{split}$ (57)

When $V_{\mathcal{I}}=0$, $H$ has both inversion symmetry,
$\mathcal{I}=\sigma^{z}$, and (spinless) time-reversal symmetry,
$\mathcal{T}=K$. In this symmetric regime, $m$ can be chosen to produce a
semimetal with Dirac points located at, for example,
$(k_{x},k_{y})=(\pm\pi/(2a_{x}),0),$ when $m=1$. In the semimetal phase,
turning on $V_{\mathcal{I}}\sigma^{x}$, which breaks inversion while
preserving $\mathcal{T}$, generates a mass term that opens a gap at the Dirac
points. The signs of the Berry curvature localized near the two now-gapped
Dirac points are opposite, as shown in Fig. 4(a), with the sign at a
particular point determined by the sign of the perturbation $V_{\mathcal{I}}$.
Hence the total Berry curvature of the occupied band integrated over the
entire BZ, equivalent to the Chern number, is zero, and hence the Berry
curvature dipole is well-defined.

To confirm our analytic calculations of the response coefficients we will
first calculate the momentum density localized around an out-of-plane magnetic
flux $\Phi_{z}$ using the tight-binding model Eq. (57). In order to determine
the $k_{x}$ momentum density in the lattice model, we must introduce magnetic
flux in a fashion that preserves translation symmetry in the
$\hat{x}$-direction. We show the configuration that we employ in Fig. 4(b).
This configuration keeps the crystal momentum $k_{x}$ as a good quantum number
and allows us to compute the value of $\mathcal{J}^{t}_{x}$ as the probability
density of the occupied single particle states weighted by their momentum
$\hbar k_{x}$. The results of the numerical calculations are presented in Fig.
4(c,d), where we study how the excess $k_{x}$ momentum density bound to
magnetic flux behaves as a function of both the magnetic flux $\Phi_{z}$ at
fixed Berry curvature dipole $\mathcal{P}_{x},$ and and as a function of
$\mathcal{P}_{x}$ at fixed $\Phi_{x}.$ Our numerical results match our
analytic calculations precisely.

We can interpret this result by noting that the momentum current in Eq. 56 can
be obtained in the semiclassical limit by considering the momentum current
carried by electron wavepackets subject to an anomalous velocity [72, 73]. The
equation of motion of an electron wavepacket with momentum $\mathbf{k}$ formed
from a single band is

$v^{i}(\mathbf{k})=\frac{\partial\mathcal{E}}{\hbar\partial
k_{i}}+\frac{e}{\hbar}\epsilon^{ij}E_{j}\mathcal{F}^{xy}(\mathbf{k}),$ (58)

where $v^{i}(\mathbf{k})$ is the wavepacket velocity,
$\mathcal{E}(\mathbf{k})$ is the energy spectrum of the band, $E_{j}$ is the
electric field, and
$\frac{e}{\hbar}\epsilon^{ij}E_{j}\mathcal{F}^{xy}(\mathbf{k})$ is the
anomalous velocity. The momentum current of the occupied states is obtained by
adding up the contributions $\hbar k_{\mu}v^{i}(\mathbf{k})$ in the BZ and
contains a term arising from the anomalous velocity given by

$\displaystyle\mathcal{J}^{i}_{\lambda}$
$\displaystyle=-\frac{e}{(2\pi)^{2}}\epsilon^{ij}E_{j}\int
dk_{x}dk_{y}\,k_{\lambda}\mathcal{F}^{xy}(k_{x},k_{y})$ (59)
$\displaystyle=-\frac{e}{4\pi}\mathcal{P}_{\lambda}\epsilon^{ij}E_{j}.$

We can also numerically probe our response equations by studying the charge
response to the deformation of the lattice. To do so, we introduce a
translation flux to rows of plaquettes located near $y=N_{y}/4$ and
$y=3N_{y}/4$, analogous to the magnetic flux configuration we just considered.
This effectively inserts two rows of dislocations such that if one encircles a
plaquette containing translation flux, the Burgers vector is in the
$x$-direction. This effectively creates opposite translational magnetic fields
$\mathcal{B}^{x}_{z}=\partial_{x}\mathfrak{e}_{y}^{x}-\partial_{y}\mathfrak{e}_{x}^{x}$
penetrating the two rows of plaquettes. Again, we choose this geometry since
it is compatible with translation symmetry in the $x$-direction. In our
lattice model we insert the translation flux by explicitly adding generalized
Peierls’ factors that are $k_{x}$-dependent, i.e.,
$\exp\left(ik_{x}\int\mathfrak{e}^{x}_{i}dx^{i}\right)$ such that the colored
regions in Fig. 4(b) contain non-vanishing translation flux. The resulting
electron charge density localized on the translation magnetic flux has a
dependence on both the $\mathcal{B}^{x}_{z}$ field strength and the Berry
curvature dipole moment $\mathcal{P}_{x}$ as shown in Fig. 4(e),(f). This
again matches the expectation from our analytic response equations.

We emphasize that the effective action (55) describes the _mutual_ _bulk_
response between the electromagnetic and the momentum currents in semimetallic
and insulating systems with vanishing Chern number. We showed in Sec. III.1.1
that one must be careful when comparing this response to the charge
polarization. In particular, our numerics show that, even in the presence of
significant inversion-breaking, the bulk momentum density response to a
magnetic flux tracks the value of the coefficient $c_{\alpha}$ from Eq. 9 as
demonstrated in Fig. 4 (d). In contrast, as shown in Sec. III.1.1, the
expression for the electric polarization, Eq. 15, contains an additional term
that is proportional to the value of a Wilson loop along the boundary of the
BZ. This value is not quantized when inversion symmetry is broken, and, for
large values of $V_{\mathcal{I}}$, this contribution becomes significant
enough that the polarization response clearly deviates from the result one
would expect from a naive interpretation of Eq. 55. However, the mutual
response between the electromagnetic and translation gauge fields described by
this action remains valid. This subtlety is not the focus of our current
article, so we leave further discussions to future work.

Figure 4: (a) Plot of the Berry curvature across the 2D Brillouin zone for the
Dirac node dipole semimetal model (57) for $m=1.1$ with an added inversion-
breaking perturbation with $V_{\mathcal{I}}=-0.5$. We use this model to probe
the $k_{x}$ momentum density response. For that we consider a completely
periodic system and insert the magnetic flux $\Phi_{z}$ thorough two lines of
plaquettes such that the translational symmetry along the $\hat{x}$-direction
is preserved, as shown in panel (b). (c) shows the $k_{x}$ momentum density
localized around one line of plaquettes penetrated by the magnetic field
$B_{z}$ as a function of magnetic flux. (d) shows the $k_{x}$ momentum density
as a function of Berry curvature dipole moment $\mathcal{P}_{x}$ defined in
Eq. (9) which we tune in our model by varying the parameter $m$ between
$m=1.0$ and $m=1.5$. In (e) and (f) we show analogous calculations for the
charge density response to a translation flux with Burgers vector in the
$x$-direction as a function of (e) translation flux at fixed Berry curvature
dipole, and (f) Berry curvature dipole at fixed translation flux. The open
circles in (e) represent Burgers’ vector choices that are not integer
multiples of a lattice constant. The red dashed lines in (c)-(f) are guides to
the eye indicating a slope of 1.

###  IV.2 2D Dirac quadrupole semimetal

Next, we consider the class of 2D semimetallic phases characterized by the
quadrupole moment of the Berry curvature introduced in Section III.1.2. We
know from Section III.1.2 that the low-energy effective response action for
this system takes the form:

$S=\frac{\hbar}{8\pi}\mathcal{Q}_{\alpha\beta}\int\mathfrak{e}^{\alpha}\wedge
d\mathfrak{e}^{\beta}.$ (60)

This action generates a momentum current response

$\mathcal{J}^{\mu}_{\alpha}=-\frac{\hbar}{4\pi}\mathcal{Q}_{\alpha\beta}\varepsilon^{\mu\nu\sigma}\partial_{\nu}\mathfrak{e}^{\beta}_{\sigma}$
(61)

These currents describe both a bulk momentum polarization (e.g., yielding
momentum on the boundary where $\mathcal{Q}_{\alpha\beta}$ changes), and a
bulk energy-momentum response to translation gauge fields. We note that this
response is exactly analogous to that of the Dirac node dipole semimetal
discussed above if we replace the electromagnetic field with a translation
gauge field.

To illustrate and explicitly confirm the responses numerically we use the
following 2-band square lattice Bloch Hamiltonian with next-nearest-neighbor
hopping terms:

$\begin{split}H(\textbf{k})&=V_{\mathcal{T}}\sigma^{x}+\sin(k_{x}a)\sin(k_{y}a)\sigma^{y}\\\
&+(m-\cos(k_{x}a)-\cos(k_{y}a))\sigma^{z}.\end{split}$ (62)

This model has an inversion symmetry (i.e., $C_{2}^{z}$ symmetry) that is
realized trivially on-site with $\mathcal{I}=\mathbb{I}$, mirror symmetry
along the $k_{x}=k_{y}$ axis, and, when $V_{\mathcal{T}}=0$, time-reversal
symmetry $\mathcal{T}=\sigma^{z}K$. This model can be tuned to a semimetal
phase as well, for example, setting $m=1$ we find four gapless Dirac points
located at $(k_{x},k_{y})=(\pm\pi/2a,0)$ and $(k_{x},k_{y})=(0,\pm\pi/2a)$.

To confirm the response action is correct, we first need to calculate the
Dirac-node quadrupole moment. To see that the Berry curvature quadrupole
moment is well-defined, we first note that the choice of $V_{\mathcal{T}}$ as
a mass perturbation forces $\mathcal{P}_{\alpha}$ to vanish. We also need the
Chern number to vanish, which is guaranteed by the mirror symmetry. With these
symmetries, the Berry curvature peaks at Dirac points that are related by
inversion symmetry have the same sign, while the peaks related by mirror
symmetry carry opposite signs, resulting in a quadrupolar distribution of the
Berry curvature, as in Fig. 5(b). Since the Chern number and
$\mathcal{P}_{\alpha}$ both vanish, the quadrupolar distribution is well-
defined and signals the presence of a well-defined elastic response in this
model (see also Ref. 63). The diagonal elements of the Dirac-node quadrupole
moment of our model are equal and opposite,
$\mathcal{Q}_{xx}=-\mathcal{Q}_{yy},$ and the off-diagonal elements are zero.
Since the sign of the Berry curvature flux for 2D Dirac points with
$\mathcal{TI}$-symmetry is ambiguous, we once again treat our system in the
insulating regime with non-zero $V_{\mathcal{T}}$ first and then recover the
semimetallic case by taking the limit $V_{\mathcal{T}}\to 0$.

Using this model, let us first focus on the momentum polarization response and
highlight the difference with the 2D Dirac node dipole semimetal case from
Section IV.1. If the bulk has a momentum polarization we expect translation-
symmetric edges to have a bound momentum density. We will first make a general
argument for the existence of the boundary momentum and then confirm the
results numerically for our model. Let us suppose our system has a boundary
normal to the $y$-direction. We expect such a boundary will carry $k_{x}$
momentum if $\mathcal{Q}_{xx}\neq 0.$ To show this, let us make a gauge
transformation on the fields in Eq. 60:
$\mathfrak{e}_{\mu}^{a}\to\mathfrak{e}_{\mu}^{a}+\partial_{\mu}\lambda^{a}$
for some vector function $\lambda^{a}.$ Since there is a boundary, the
response action is not gauge invariant and we find the variation
$\delta_{\lambda}S=-\frac{\hbar\mathcal{Q}_{ab}}{8\pi}\lambda^{a}(\partial_{0}\mathfrak{e}_{x}^{b}-\partial_{x}\mathfrak{e}_{0}^{b}).$
Our system has no translation-twisting of the boundaries, i.e.,
$\mathfrak{e}_{x}^{y}=\mathfrak{e}_{y}^{x}=0$, so we find the variation
reduces to
$\delta_{\lambda}S=-\frac{\hbar\mathcal{Q}_{xx}}{8\pi}\lambda^{x}(\partial_{0}\mathfrak{e}_{x}^{x}-\partial_{x}\mathfrak{e}_{0}^{x}).$
This variation can be canceled by adding an action of the form Eq. 34. That
is, we expect to have 1D degrees of freedom on the boundary that harbor a non-
vanishing $k_{x}$-momentum density captured by an effective 1D quadrupole
moment $\mathcal{Q}_{xx}$ that matches the value of the 2D quadrupole moment.
Interestingly, we note that the coefficient of Eq. 34 is twice that of the
variation we need to cancel. Hence, the edge of our system has a fractional
momentum density, i.e., a 1D system with the same $\mathcal{Q}_{xx}$ would
have twice as much momentum. This is analogous to the fractional boundary
charge density one finds from the half-quantized electric charge polarization.

We confirm this response numerically by studying the model (62) on a lattice
in a ribbon geometry that is open in the $\hat{y}$-direction and periodic in
$\hat{x}$. Figure 5 (a) shows the resulting band structure, for which a gap is
opened by a non-vanishing $V_{\mathcal{T}}$ and the occupied states have two
symmetrically positioned sets of flat band states: one in an interval having
$k_{x}<0$ and the other in an interval having $k_{x}>0$. The occupied boundary
states with $k_{x}<0$ (red) are localized near the top ($y=N_{y}$) boundary,
while the occupied boundary states with $k_{x}>0$ (blue) are localized near
the bottom ($y=1$) boundary. At half filling we find that the excess/deficit
charge near the boundary depends on $k_{x}$ as shown in Fig. 5(c). We see that
the states at positive and negative $k_{x}$ are imbalanced, indicating a non-
vanishing $k_{x}$ momentum density on the edge. Indeed, each state between the
Dirac nodes contributes an amount to the total edge momentum equal to $k_{x}$
weighted by a factor of $\pm 1/2$, since the particle density on the edge at
each $k_{x}$ in this range is $\pm 1/2.$ Because states at opposite $k_{x}$
have opposite excess/deficit probability density, the total sum is non-
vanishing and depends on $\mathcal{Q}_{xx}$ as shown in Fig. 5(f). We find
that the bulk momentum polarization
$P^{y}_{k_{x}}=\frac{\hbar\mathcal{Q}_{xx}}{8\pi}$ matches the numerically
calculated boundary momentum density, as expected for a generalized surface
charge theorem 555We comment that even though the Chern-Simons term for the
translation gauge fields shares some properties with the electromagnetic
Chern-Simons term, there is a key distinction: The translation gauge fields
have a constant background. This allows the Dirac node quadrupole system to
have a static momentum polarization, whereas the electromagnetic Chern Simons
term in a Chern insulator would predict generating an electric polarization as
one tunes the vector potential.. To further probe the response equations, we
subject the Dirac node quadrupole semimetal to the same linear array of
dislocations employed in the previous subsection (c.f. Fig. 4(b)). From Eq.
61, we expect to find momentum density localized on dislocations. Since our
geometry preserves translation in the $\hat{x}$-direction, we can compute the
amount of $k_{x}$ momentum bound to dislocations, similar to how we computed
the amount of charge bound to dislocations in the previous subsection. We show
our results in Fig. 5(d)(e) where we first plot momentum density as a function
of $\mathcal{Q}_{xx}$ for fixed translation flux $\mathcal{B}_{z}^{x},$ and
then plot momentum density as a function of $\mathcal{B}_{z}^{x}$ for fixed
$\mathcal{Q}_{xx}.$ Both results match the analytic value from the response
action.

Finally, let us briefly consider a case when the mixed energy-momentum
quadrupole moment $\mathcal{Q}_{ta}$ is non-vanishing. In this scenario the
effective action (60) implies the existence of a bulk orbital _momentum
magnetization_ of

$M^{z}_{k_{\mu}}=-\frac{\hbar}{8\pi}\mathcal{Q}_{t\mu},$ (63)

that will manifest as boundary momentum currents, even in equilibrium (note we
assume $\mathfrak{e}_{t}^{t}=1)$. To generate a non-vanishing
$\mathcal{Q}_{t\mu}$ in our model (62), we turn on an additional perturbation

$\Delta H(\textbf{k})=\epsilon\sin(k_{x})\mathbb{I}_{2\times 2}.$ (64)

When $m=1$ and $V_{\mathcal{T}}\to 0_{-}$, this induces
$\mathcal{Q}_{tx}=-\pi\epsilon$ and $\mathcal{Q}_{tt}=\epsilon^{2}$, leading
to momentum $k_{x}$ magnetization,
$M^{z}_{k_{x}}=-\frac{\hbar}{8\pi}\mathcal{Q}_{tx}$, and bulk energy
magnetization, $M^{z}_{k_{t}}=-\frac{\hbar}{8\pi}\mathcal{Q}_{tt}$, following
from Eq. (63). In Fig. 5(g) we plot the boundary energy current response
$\Delta\mathcal{J}^{x}_{t}$ as a function of $\mathcal{Q}_{tt}.$ We calculate
this quantity by summing the particle current $\frac{1}{\hbar}\frac{\partial
H}{\partial k_{x}}$ weighted by the energy $\epsilon(k)$ of each state. The
slope of the plot confirms the coefficients predicted in Eq. (63).

Figure 5: (a) Spectrum of the 2D Dirac node quadrupole semimetal (62) in a
ribbon geometry ($y$-direction open, $x$-direction periodic) for $m=1$, the
$\mathcal{T}$-breaking perturbation set to $V_{\mathcal{T}}=-0.2,$ and the
energy tilt in Eq. 64 $\epsilon=0.1$. At half filling, the ground state of the
model is momentum-polarized: occupied states localized near $y=1$, which are
indicated by the blue color, carry a positive value of the $k_{x}$ momentum,
while the occupied states near $y=N_{y}$ have a negative value of $k_{x}$. (b)
Berry curvature distribution across the Brillouin zone for a small gapping
perturbation $V_{\mathcal{T}}=-0.2.$ (c) The boundary charge distribution as a
function of momentum. (d), (e) $k_{x}$ momentum bound to a row of dislocations
(c.f. Fig. 4(b)) as a function of $\mathcal{Q}_{xx}$ at fixed
$\mathcal{B}_{z}^{x}$ in (d) and as a function of $\mathcal{B}_{z}^{x}$ at
fixed $\mathcal{Q}_{xx}$ in (e). (f) Plot of momentum polarization
$P^{y}_{k_{x}}$ obtained from computing $k_{x}$-momentum bound to an edge
normal to $\hat{y}.$ (g) As a consequence of non-zero $\epsilon$ we see that
the velocities of single-particle states in (a)localized on opposite edges
have the same sign, while the energy and $k_{x}$ momentum charges are exact
opposite. This leads to boundary energy currents as illustrated in panel (g)
as a function of $\mathcal{Q}_{tt}.$

###  IV.3 3D nodal line dipole semimetal

Heuristically we can consider nodal 3D semimetals as arising from stacks of 2D
Dirac node dipole semimetals. Furthermore, similar to the 2D case, with
inversion symmetry the bulk response action

$S[\mathfrak{e}^{\lambda},A]=e\mathcal{B}_{\mu\nu}\int\mathfrak{e}^{\mu}\wedge\mathfrak{e}^{\nu}\wedge
dA$ (65)

can be interpreted as a charge magnetization $M_{i}$ and electric polarization
$P_{e}^{i}$:

$e\mathcal{B}_{ta}=M^{i}\mathfrak{e}_{i}^{a},\quad
e\mathcal{B}_{ab}=\varepsilon_{ijk}P^{k}_{e}\mathfrak{e}_{i}^{a}\mathfrak{e}_{j}^{b}$
(66)

where we have taken functional derivatives of Eq. 65 with respect to the
magnetic and electric fields respectively, and used $\mathfrak{e}_{t}^{t}=1.$
For an unmodified geometry we recover the results of Ref. 35, i.e.,

$e\mathcal{B}_{ta}=M^{a},\,\,\,e\mathcal{B}_{ab}=\varepsilon_{abk}P^{k}_{e}.$
(67)

Microscopically, the coefficient $\mathcal{B}_{ab},$ where $a,b=1,2,3,$ is
proportional to the area of the line nodes that project onto surfaces normal
to the $ab$-plane as illustrated in Fig. 6(a).

Figure 6: (a) Fermi line of a 3D NLSM (70) with $V_{\mathcal{I}}=0$, $m=2$
that is tiled in the energy-momentum space $\\{k_{z},k_{x},E\\}$ by the
perturbation (72) where we set $\epsilon=1$. The projections of this curve
onto the $\\{k_{x},k_{z}\\}$ and $\\{k_{z},E\\}$ planes give the exact values
of the $\mathcal{B}_{xz}$ and $\mathcal{B}_{tz}$ coefficients respectively.
(b) A screw dislocation characterized by a Burgers’ vector $b^{z}=a_{z}$
creates an internal boundary carrying a current circulating around the
magnetization vector $M_{z}$. Note that the currents’ direction is
perpendicular to the Burgers’ vector and the Magnetization vector $M_{z}$, as
predicted by Eq. 73.

The bulk action also implies a non-vanishing momentum response to
electromagnetic fields:

$\mathcal{J}^{\mu}_{\lambda}=2e\mathcal{B}_{\lambda\eta}\varepsilon^{\mu\nu\rho\sigma}\mathfrak{e}^{\eta}_{\nu}\partial_{\rho}A_{\sigma},$
(68)

and a conjugate electromagnetic response to translation gauge fields:

$j^{\mu}=2e\mathcal{B}_{\lambda\eta}\varepsilon^{\mu\nu\rho\sigma}\mathfrak{e}^{\lambda}_{\nu}\partial_{\rho}\mathfrak{e}^{\eta}_{\sigma}.$
(69)

To illustrate how these responses manifest in an explicit model, we can
construct a Hamiltonian for a 3D nodal line dipole semimetal by stacking
copies of the 2D Dirac node dipole semimetal in Eq. (57) in the
$\hat{z}$-direction. When there is no hopping between the 2D layers, such a
system will have two lines of gapless states spanning the BZ along the $k_{z}$
direction, located at $(k_{x},k_{y})=(\pm K,0)$ (for our model). Adding
hopping terms in the $\hat{z}$-direction leads to a Bloch Hamiltonian:

$\begin{split}H(\textbf{k})&=V_{\mathcal{I}}\sigma^{x}+\sin(k_{y}a_{y})\sigma^{y}\\\
&+(m-\cos(k_{x}a_{x})-\cos(k_{y}a_{y})-\cos(k_{z}a_{z}))\sigma^{z}.\end{split}$
(70)

Taking $V_{\mathcal{I}}\to 0$ and $m=2$, we find a single loop of gapless
states located in the $k_{y}=0$ plane, described by the equation
$\cos(k_{x}a_{x})+\cos(k_{z}a_{z})=1$. The stack of 2D Dirac node dipole
semimetals will naturally endow the 3D nodal line system with electric
polarization (and/or magnetization). Correspondingly, this model has a single
non-zero component of the antisymmetric tensor $\mathcal{B}_{xz}$ defined in
Eq. (37), which encodes a charge polarization in the $\hat{y}$-direction. From
Eq. 68, a non-vanishing $\mathcal{B}_{xz}$ also implies a $k_{x}$ momentum
line-density localized on a magnetic flux tube oriented in the
$\hat{z}$-direction:

$\mathcal{J}^{t}_{x}=2e\mathcal{B}_{xz}\varepsilon^{tzij}\mathfrak{e}^{z}_{z}B_{z}=2e\mathcal{B}_{xz}B^{z},$
(71)

similar to a stack of decoupled 2D Dirac semimetallic layers (in the last
equality we replaced $\mathfrak{e}_{z}^{z}=1)$. This is the 3D analog of the
response shown in Figs. 4(c) and (d) for the 2D Dirac semimetal.

We can see an example of a charge response if we tilt the nodal line to
introduce a non-zero value of $\mathcal{B}_{tz}$ as illustrated in Fig. 6(a).
In our model we can tilt the node by adding an extra dispersion

$\Delta H(\textbf{k})=\epsilon\sin(k_{x}a_{x})\mathbb{I}_{2\times 2},$ (72)

to the Hamiltonian. This term breaks $\mathcal{T}$ and induces a net
magnetization $M_{z}=e\mathcal{B}_{tz}$, setting up the corresponding
circulating boundary currents in the system [35].

Now, when $\mathcal{B}_{tz}$ is non-vanishing, Eq. 69 implies that a screw
dislocation with Burgers vector $b^{z}\hat{z}$ hosts a bound electromagnetic
current. Indeed, if we assume the screw dislocation is located at
$(x,y)=(0,0)$ and runs along the $z$-axis we find

$j^{z}=-2e\mathcal{B}_{tz}\varepsilon^{tzjk}\mathfrak{e}^{t}_{t}\partial_{j}\mathfrak{e}^{z}_{k}=-2e\mathcal{B}_{tz}b^{z}\delta(x)\delta(y),$
(73)

where we used $\mathfrak{e}_{t}^{t}=1$ and
$\nabla\times\mathfrak{e}^{z}=b^{z}\delta(x)\delta(y).$

We can illustrate the origin of this current by considering the magnetization
$M_{z}$ (and associated boundary currents) induced by $\mathcal{B}_{tz}.$ A
screw dislocation with Burgers vector $b^{z}\hat{z}$ can be constructed by
cutting a seam into layers normal to $\hat{z}$ and re-gluing them along the
seams with neighboring layers above or below. When cut, the boundary current
associated to $M_{z}$ will appear, and after re-gluing this current will be
routed vertically along the screw-dislocation line, i.e., along the
$z$-direction as shown in Fig. 6 (b). The magnetization $M_{z}$ gives rise to
a surface bound current $j_{\partial}=M_{z}$ circulating around the
$\hat{z}$-axis in each layer. The effective number of current loops winding
around the dislocation line per unit length is equal to the Burgers vector
$b^{z}$. Thus the total current in the $\hat{z}$-direction is:

$j^{z}=-b^{z}j_{\partial}=-2e\mathcal{B}_{tz}b^{z},$ (74)

which reproduces the result obtained from the response action. Furthermore, we
can understand the sign of the current from Fig. 6(b) where we see that the
current on the dislocation has an opposite orientation to the current
generated by $M_{z}.$ Another interesting consequence of Eq. (65) is the
topological piezoelectric effect discussed in Ref. [50].

###  IV.4 3D nodal line quadrupole semimetal

In Sec. III.3, we derived the effective response action:

$S[\mathfrak{e}^{\lambda}]=\hbar\mathcal{B}_{\lambda\eta,\alpha}\int\mathfrak{e}^{\lambda}\wedge\mathfrak{e}^{\eta}\wedge
d\mathfrak{e}^{\alpha}.$

for the nodal line quadrupole semimetal. The response action implies the
energy-momentum currents:

$\mathcal{J}^{\mu}_{\lambda}=2\hbar\left(\mathcal{B}_{\lambda\eta,\alpha}-\mathcal{B}_{\eta\alpha,\lambda}\right)\varepsilon^{\mu\nu\rho\sigma}\mathfrak{e}^{\eta}_{\nu}\partial_{\rho}\mathfrak{e}^{\alpha}_{\sigma},$
(75)

where we have used that $\mathcal{B}_{\lambda\eta,\alpha}$ is anti-symmetric
under exchange of the first two indices.

In analogy with the 2D Dirac node dipole and Dirac node quadrupole semimetals,
we expect that most of the responses from the Dirac nodal line dipole
semimetal in Sec. IV.3 can be translated to describe some of the responses of
this action if we replace charge currents and densities with momentum currents
and densities etc. Indeed, we showed in Eq. 39 that when $\lambda$ and $\eta$
are both spatial indices, $\mathcal{B}_{\lambda\eta,\lambda}$ implies a
momentum polarization in a direction perpendicular to $\lambda$ and $\eta,$
and carrying momentum parallel to $\lambda$. By analogy, the mixed temporal-
spatial components $\mathcal{B}_{it,j}$ describe a momentum magnetization in
the $i$-th direction carrying momentum in the $j$-th direction. The momentum
magnetization is further responsible for generating bound-currents on screw-
dislocations, i.e., the momentum magnetization will have circulating boundary
momentum currents and a momentum current along screw dislocations similar to
the charge bound currents on dislocations shown in Section IV.3.

Figure 7: Fermi Lines of the model (76) with $m=2$ and $V_{\mathcal{T}}\to
0_{-}$. Resolving this structure as a pair of loops with _fixed orientation_
we can project them onto the $k_{x}k_{z}$ or $k_{y}k_{z}$ surfaces to
determine the momentum polarization. The colored regions of the projected
nodes indicate flat drumhead states that would appear in open boundary
conditions on one boundary (red) or the opposing boundary (blue). By looking
at the relative positions of the two areas bounded by the projected loops in
the surface BZ, we see that one surface will have one sign of the $k_{x}$ or
$k_{y}$ momentum, and the other surface will have the other. For example, for
the $k_{x}k_{z}$ surface BZ the the projections indicate a dipole moment of
$k_{x}$ momentum polarized along the $y$ direction captured by the response
coefficient $\mathcal{B}_{zx,x}$. Inset: Cage-like nodal Fermi surface in the
model (76) with $E_{F}=0.2$.

To be more explicit, we can illustrate the momentum polarization in a model by
showing the analog of the surface charge theorem, i.e., momentum polarization
will yield surface momentum densities. To obtain a Hamiltonian for the nodal
line quadrupole semimetal, we begin by stacking 2D Dirac node quadrupole
semimetals (see Fig. 5 (b)) along the $\hat{z}$-direction. When the planes are
completely decoupled, this construction produces a set of four straight Fermi
lines stretching in the $k_{z}$-direction. If we couple the two-dimensional
planes, then we arrive at the following Bloch Hamiltonian:

$\begin{split}H(\textbf{k})&=V_{\mathcal{T}}\sigma^{x}+\sin(k_{x}a)\sin(k_{y}a)\sigma^{y}\\\
&+(m-\cos(k_{x}a)-\cos(k_{y}a)-\cos(k_{z}a_{z}))\sigma^{z}.\end{split}$ (76)

For a wide range of parameters this model has a pair of nodal line loops that
form a cage structure as shown in Figs. 2 and 7 with $m=2$ and
$V_{\mathcal{T}}=0.$ In general, the local gaplessness of the nodal loops can
be protected by the product $\mathcal{TI}.$ The cage structure created by the
joined, intersecting loops can be split apart by, for example, breaking mirror
symmetry along the $k_{x}=k_{y}$ axis while preserving $\mathcal{TI}$.
However, even in this case the nodal loops still produce a non-vanishing
contribution to the response coefficient $\mathcal{B}_{\alpha\beta,\gamma}.$
Hence, the response is more general than the specific cage-like nodal
configuration. Calculating the response coefficient for the action in the
limit $V_{\mathcal{T}}\rightarrow 0_{-}$, we find that
$\mathcal{B}_{xz,x}=-\mathcal{B}_{zx,x},$
$\mathcal{B}_{yz,y}=-\mathcal{B}_{yz,y}$ are non-vanishing, as shown in Fig.
7.

Using this model we can illustrate the origin of the boundary momentum
resulting from the bulk momentum polarization. The discussion is analogous to
the calculation of the boundary momentum of the 2D Dirac node quadrupole
semimetal in Sec. IV.2. Indeed, the analogy is clear since the cage nodal
structure is just arising from a family of 2D Dirac node quadrupoles
parameterized by $k_{z}.$ To specify an unambiguous momentum polarization we
turn on a small $\mathcal{T}$-breaking perturbation $V_{\mathcal{T}}.$ After
doing this, and as shown in Fig. 7, we see that the two nodal loop segments
that lie in the $k_{y}=0$ plane (one for $k_{x}>0$ and one for $k_{x}<0$)
carry the same Berry flux in the $k_{z}$-direction (red arrows in Fig. 7).
Similarly, the two loop segments in the $k_{x}=0$ plane carry the same Berry
flux (blue arrows), which is opposite to that carried by the $k_{y}=0$
segments. Consequently, the loop segment in the $k_{y}=0$, $k_{x}>0$ half-
plane must connect with a loop segment in the $k_{x}=0$ plane in order to form
a closed nodal loop with a consistent helicity/flux sign.

To clarify the consequences of this nodal configuration let us consider the
$k_{x}k_{z}$ plane in Fig. 7. We can calculate a Berry-Zak phase [25] in the
$k_{y}$ direction parameterized by $(k_{x},k_{z}),$ and for our model we find
a Berry phase of magnitude $\pi$ inside the projected nodal region in the
$k_{x}k_{z}$ plane. When $V_{\mathcal{T}}$ is turned on, the sign of the $\pi$
Berry-Zak phases are no longer ambiguous, and are opposite for the projected
areas at $k_{x}>0$ and $k_{x}<0.$ If we calculate the total polarization in
the $y$-direction when summed over all $k_{x}$ and $k_{z}$ it will vanish.
However, the polarization weighted by the $k_{x}$ momentum will be non-zero.
The occupied drumhead surface states in the $k_{x}k_{z}$ surface-BZ (see Fig.
7 and c.f. Fig. 5(a,b,c)) will have an imbalanced $k_{x}$ momentum, but, when
combined with the bulk charge density, a vanishing charge (c.f. Fig. 5(c)).
This is a reflection of the surface charge theorem for a vanishing charge
polarization, and non-vanishing momentum polarization. We numerically
calculated the magnitude of the bound surface momentum, finding it to be in
agreement with the value predicted by the response action,
$2\hbar\mathcal{B}_{xz,x}.$ We see from this picture that to have a non-zero
response $\mathcal{B}_{xz,x}$, we want two oppositely oriented nodal loops
with identical, non-vanishing areas when projected in the $k_{x}k_{z}$-plane,
but positioned so that the sums of all $k_{x}$ inside each nodal loop are
different from each other, e.g., in our model they are opposite values.

As an additional explicit example of a non-vanishing response allowed in our
model we can consider the momentum density

$\mathcal{J}^{0}_{x}=2\hbar\mathcal{B}_{xz,x}\epsilon^{ijk}(2\mathfrak{e}_{i}^{z}\partial_{j}\mathfrak{e}_{k}^{x}-\mathfrak{e}_{i}^{x}\partial_{j}\mathfrak{e}_{k}^{z})$
(77)

generated by a geometric deformation. To generate a non-vanishing response let
us consider an $xz$-planar interface. Since we must preserve translation
symmetry along $x$ to calculate $k_{x}$ momentum, and we want to preserve
translation in $z$ for convenience, we have the following terms:

$\mathcal{J}^{0}_{x}=2\hbar\mathcal{B}_{xz,x}\left(2\mathfrak{e}_{x}^{z}\partial_{y}\mathfrak{e}_{z}^{x}-2\mathfrak{e}_{z}^{z}\partial_{y}\mathfrak{e}_{x}^{x}-\mathfrak{e}_{x}^{x}\partial_{y}\mathfrak{e}_{z}^{z}+\mathfrak{e}_{z}^{x}\partial_{y}\mathfrak{e}_{x}^{z}\right).$

If we cut the system at $y=0$, both sides of the interface will carry a
surface $k_{x}$-momentum density $\mathcal{J}^{0}_{x,surf}=\pm
2\hbar\mathcal{B}_{xz,x},$ since the system has a $k_{x}$ momentum
polarization along $\hat{y}$ with this magnitude. Since each interface carries
an opposite sign of the momentum density, if we glue them back together there
will be no momentum at the interface. Now, for $y>0$ let us perturb away from
the background translation gauge field configuration to
$\mathfrak{e}_{i}^{a}=(1+\epsilon^{a})\delta_{i}^{a}$ where
$\epsilon^{a}=(\epsilon^{x},\,0,\,\epsilon^{z})$ is a small deformation. The
momentum density response to leading order in $\epsilon^{a}$ is

$\mathcal{J}^{0}_{x}=2\hbar\mathcal{B}_{xz,x}\left[-2\epsilon^{x}\delta(y)-\epsilon^{z}\delta(y)\right],$
(78)

which we see is localized at the interface $y=0.$

We can interpret this response by noting that changing $\mathfrak{e}_{x}^{x}$
or $\mathfrak{e}_{z}^{z}$ effectively changes the area of one side of the
interface ($y>0$) relative to the other ($y<0$). Since the total $k_{x}$
momentum on both sides of the interface should be unchanged by this
deformation (we maintain translation symmetry in $x$ during the process), then
increasing the area for $y>0$ must lower the momentum _density_. Indeed, the
surface $k_{x}$ momentum density on $\hat{y}$ surfaces must be inversely
proportional to $L_{x}$ and $L_{z}.$ Finally, since we are considering
$k_{x}$-momentum density, the quantization of which depends on $L_{x}^{-1}$,
$\mathcal{J}^{0}_{x}$ actually depends on $L_{x}^{-2}$, hence the difference
between the coefficients of $\epsilon^{x}$ and $\epsilon^{z}$ in Eq. 78.

###  IV.5 3D Weyl node dipole semimetal

The electromagnetic and geometric response of time-reversal breaking 3D Weyl
semimetals have been discussed extensively in the literature [33, 59, 34, 75,
17, 18, 13, 76, 77, 78, 79, 80, 81, 60, 12, 15, 16, 19, 20, 21, 23, 82, 83,
84, 85]. Here we focus on a few particular consequences of the mixed
crystalline-electromagnetic response and the matching between the response
field theory and microscopic lattice model calculations. Recall that the
response action for a 3D Weyl semimetal with a non-vanishing Weyl-node dipole
moment $\mathcal{P}_{\lambda}$ is

$S[\mathfrak{e}^{\lambda}_{\nu},A_{\mu}]=\frac{e^{2}\mathcal{P}_{\lambda}}{8\pi^{2}\hbar}\int\mathfrak{e}^{\lambda}\wedge
A\wedge dA.$ (79)

This response implies the following bulk electromagnetic and momentum
currents:

$\displaystyle
j^{\mu}=-\frac{e^{2}\mathcal{P}_{\lambda}}{4\pi^{2}\hbar}\varepsilon^{\mu\nu\rho\sigma}\mathfrak{e}^{\lambda}_{\nu}\partial_{\rho}A_{\sigma}+\frac{e^{2}\mathcal{P}_{\lambda}}{8\pi^{2}\hbar}\varepsilon^{\mu\nu\rho\sigma}A_{\nu}\partial_{\rho}\mathfrak{e}^{\lambda}_{\sigma}\
,$ (80)
$\displaystyle\mathcal{J}^{\mu}_{\lambda}=\frac{e^{2}\mathcal{P}_{\lambda}}{8\pi^{2}\hbar}\epsilon^{\mu\nu\rho\sigma}A_{\nu}\partial_{\rho}A_{\sigma}.$
(81)

In the presence of dislocations the translational flux is non-vanishing, and
hence the bulk electromagnetic current is anomalous:

$\begin{split}\partial_{\mu}j^{\mu}=-\frac{e^{2}\mathcal{P}_{\lambda}}{8\pi^{2}\hbar}\varepsilon^{\mu\nu\sigma\rho}\partial_{\mu}\mathfrak{e}^{\lambda}_{\nu}\partial_{\sigma}A_{\rho}.\end{split}$
(82)

This reflects the fact that the action Eq. 79 is not gauge-invariant in the
presence of dislocations. Indeed, in our explicit tight-binding model
calculations below we find the spectrum on a single screw dislocation line
contains a pair of chiral modes of the same chirality (one near each bulk Weyl
node momentum as shown in Fig. 9(b)). These modes are responsible for the
anomalous current on dislocation lines, as was first described by Ref. 38.

To verify the electromagnetic response to the applied crystalline gauge field
we consider a simple 2-band model of a 3D Weyl semimetal with a pair of
gapless nodes:

$\begin{split}H(\textbf{k})&=\sin(k_{z}a_{z})\sigma^{x}+\sin(k_{y}a_{y})\sigma^{y}\\\
&+(2-m-\cos(k_{x}a_{x})-\cos(k_{y}a_{y})-\cos(k_{z}a_{z}))\sigma^{z}.\end{split}$
(83)

The Weyl node with the positive chirality $\chi=+1$ is located at
$\textbf{k}=(\arccos(-m),\,0,\,0)$ and the node with $\chi=-1$ is at
$\textbf{k}=(-\arccos(-m),\,0,\,0)$. The Weyl node dipole moment therefore has
only one non-zero component $\mathcal{P}_{x}=2\arccos(-m)$ and the resulting
response action is

$S[\mathfrak{e}^{\lambda}_{\nu},A_{\mu}]=\frac{e^{2}\mathcal{P}_{x}}{8\pi^{2}\hbar}\int
d^{4}x\epsilon^{\mu\nu\rho\sigma}\mathfrak{e}^{x}_{\mu}A_{\nu}\partial_{\rho}A_{\sigma}.$
(84)

Let us first consider the response arising from the constant background
translation fields $\mathfrak{e}_{x}^{x}=1$ and
$\mathfrak{e}_{y}^{x}=b^{x}/L_{y}$, which describe a twist such that a
particle traversing the lattice in the $y$-direction translates by $b^{x}$ in
the $x$-direction. We note that such a configuration is volume preserving
since $\det({\bf{e}})=1$, where the matrix ${\bf{e}}$ has matrix elements
$e_{ij}=\mathfrak{e}_{i}^{j}.$ When $b^{x}=0$ the response action is

$\frac{e^{2}\mathcal{P}_{x}}{8\pi^{2}\hbar}\int\mathfrak{e}_{x}^{x}dx\int
dtdydz\epsilon^{x\nu\rho\sigma}A_{\nu}\partial_{\rho}A_{\sigma}.$

Using the relation $\int dx\mathfrak{e}_{x}^{x}=L_{x}$ we find an anomalous
Hall effect in the $yz$-plane such that
$\sigma_{yz}=\frac{e^{2}}{h}\frac{\mathcal{P}_{x}L_{x}}{2\pi},$ which is the
standard result [33, 32]. Now, if we turn on $b^{x}$ we will still have the
same $\sigma_{yz},$ but we will also have the additional term

$\frac{e^{2}\mathcal{P}_{x}}{8\pi^{2}\hbar}\int\mathfrak{e}_{y}^{x}dy\int
dtdxdz\epsilon^{y\nu\rho\sigma}A_{\nu}\partial_{\rho}A_{\sigma}.$

Because of the different index on the $\epsilon$-symbol, this term represents
an anomalous Hall effect in the $xz$-plane with
$\sigma_{zx}=\frac{e^{2}}{h}\frac{\mathcal{P}_{x}b_{x}}{2\pi}.$ We can find a
simple interpretation for this effect: when we turn on $\mathfrak{e}_{y}^{x}$,
the minimal coupling $k_{x}\to k_{x},k_{y}\to k_{y}+k_{x}\mathfrak{e}_{y}^{x}$
shifts the bulk Weyl nodes,
$(\pm\mathcal{P}_{x}/2,\,0,\,0)\to(\pm\mathcal{P}_{x}/2,\,\pm\mathcal{P}_{x}b_{x}/(2L_{y}),\,0).$
Hence an effective $\mathcal{P}_{y}=\frac{\mathcal{P}_{x}b_{x}}{L_{y}}$ is
generated when the Weyl momenta are sheared. Indeed, we expect that, at least
for uniform, traceless translation gauge field deformations, the response
phenomena can be simply interpreted as transformations of the Weyl node dipole
$\mathcal{P}_{i}\to\mathfrak{e}^{j}_{i}\mathcal{P}_{j}.$ We show an explicit
example of this in the first and third surface-BZ panels of Fig. 8(a) where
the bulk nodes and their connected Fermi arcs have been rotated in the
deformed geometry relative to the undeformed geometry. We note that if the
deformation is not volume-preserving, then we must be careful when considering
what is held fixed while volume is changing in order to interpret the
resulting phenomena.

Figure 8: (a) The three panels show numerically calculated Fermi arcs in
(left) the surface BZ with un-deformed geometry, (right) the surface BZ with
$\mathfrak{e}_{y}^{x}$ non-vanishing, and (center) the arcs localized at the
interface formed by gluing the two sides of the interface together. The
colored circles in the first and third panels represent the surface BZ
projections of the bulk Weyl nodes on either side of the interface. The color
is a guide to show the connectivity/orientation of the Fermi arcs, not the
chirality of the bulk nodes. On both sides of the interface the bulk nodes
have the same chirality, but since they are effectively projected onto
surfaces having opposite normal vectors they generate Fermi arcs having
opposite chirality. (b) Illustrations of (left) the un-deformed geometry and
(right) the deformed geometry with $\mathfrak{e}_{y}^{x}$ non-vanishing. (c)
The numerically calculated current localized at the interface between un-
deformed and deformed geometries as a function of the chemical potential shift
$A_{0}$. Figure 9: (a) The bulk spectrum of a Weyl semimetal with two nodes on
the $k_{x}$-axis (b) The spectrum of the same Weyl semimetal with periodic
boundary conditions and two screw dislocations with opposite Burgers vectors
threaded along the $x$-direction. Red and blue coloration indicates on which
dislocation the chiral modes are localized. Each dislocation has a net
positive (red) or negative (blue) chirality. (c) Numerical calculation of the
charge density bound to a screw dislocation as $A_{x}$ is tuned.

In addition to these cases of fixed background translation fields, let us
consider varying those fields in space. We are interested in the
electromagnetic response to applied translational _magnetic_ fields
$\mathcal{B}_{i}^{a}=\epsilon^{ijk}\partial_{j}\mathfrak{e}_{k}^{a}$. Since
the nodes in our model are separated in $k_{x}$, we will consider geometries
where the Burgers vector of the translation magnetic field also points along
the $x$-direction, $\mathcal{B}_{i}^{x}\neq 0$.

First let us consider a system containing a domain wall as a function of $z$,
such that at $z=0$ the field $\mathfrak{e}_{y}^{x}$ jumps from 0 to
$b_{x}/L_{y}.$ For $z<0$ we have bulk Weyl nodes that project onto the
$z$-surface at $(\pm\mathcal{P}_{x}/2,0),$ while for $z>0$ the bulk Weyl nodes
have been transformed and sit at
$(\pm\mathcal{P}_{x}/2,\pm\mathcal{P}_{x}b_{x}/(2L_{y})).$ We show the
numerically calculated Fermi arcs for our un-deformed and deformed models in
the left and right surface BZ panels of Fig. 8(a).

Now let us glue the $z<0$ and $z>0$ sides to each other to make a domain-wall
interface. We schematically illustrate the interface geometry in Fig. 8(b).
Since the normal vector on each side of the interface is opposite, we expect
the Fermi arcs for $z<0$ to have the opposite chirality to their corresponding
arcs for $z>0.$ Indeed, as shown in the center surface BZ panel of Fig. 8(a),
the Fermi arcs on both sides can hybridize because of their opposite
chiralities and form new arcs in the 2D subsystem of the interface. These new
Fermi arcs encode the fact that the Hall conductivity $\sigma_{xz}$ is varying
at this interface. These effects are all manifestations of the fact that the
Weyl node dipole moment $\mathcal{P}_{i}$ is changing at the interface, and
hence we expect Fermi arcs to be trapped generically at the interfaces of this
type. We note that a similar strain geometry, and the corresponding Weyl node
configuration, was discussed in [49].

From Eq. 80 we see that applying a uniform, non-vanishing $A_{0}$ to the
system described above should generate a charge current in the $x$-direction.
We can see the microscopic origin of this current as follows. If we increase
$A_{0},$ each linearly-dispersing point on the Fermi arc will have an excess
charge density $\delta n({\bf{k}})=\frac{eA_{0}}{2\pi\hbar|v_{F}({\bf{k}})|}$
where $v_{F}({\bf{k}})$ is the Fermi velocity at the Fermi arc located at
${\bf{k}}$ in the surface-BZ. Hence, the contribution to the current of such a
point on the Fermi arc is $j^{x}({\bf{k}})=ev_{F}({\bf{k}})\delta
n({\bf{k}}).$ For our model and geometry, the contributions to the $j^{x}$
current that are linear in the deformations of $\mathfrak{e}_{i}^{a}$ arise
from the Fermi arcs stretching between $(K,0)\to(K,K\mathfrak{e}_{y}^{x})$ and
$(-K,0)\to(-K,-K\mathfrak{e}_{y}^{x}).$ Each of these arcs has a fixed value
$k_{x}=\pm K$ and each arc has an opposite Fermi velocity. Hence

$\begin{gathered}j^{x}=ev_{F}(K,k_{y})\delta
n\frac{K\mathfrak{e}_{y}^{x}}{2\pi}\\\ +ev_{F}(-K,k_{y})\delta
n\frac{K\mathfrak{e}_{y}^{x}}{2\pi}=\frac{e^{2}\mathcal{P}_{x}\mathfrak{e}_{y}^{x}A_{0}}{4\pi^{2}\hbar}{\rm{sgn}}(v_{F}),\end{gathered}$

where $K\mathfrak{e}_{y}^{x}/2\pi$ counts the density of states on the Fermi
arc in the $k_{y}$ direction, ${\rm{sgn}}(v_{F})$ is sign of the velocity on
the $k_{x}=+K$ arc, and $\mathcal{P}_{x}=2K^{2}$ is the un-deformed value.