# Topological Dipole Conserving Insulators and Multipolar Responses

Julian May-Mann  Department of Physics and Institute of Condensed Matter
Theory, University of Illinois at Urbana-Champaign, 1110 West Green Street,
Urbana, Illinois 61801-3080, USA  Taylor L. Hughes  Department of Physics and
Institute of Condensed Matter Theory, University of Illinois at Urbana-
Champaign, 1110 West Green Street, Urbana, Illinois 61801-3080, USA

###### Abstract

Higher order topological insulators (HOTIs) are a novel form of insulating
quantum matter, which are characterized by having gapped boundaries that are
separated by gapless corner or hinge states. Recently, it has been proposed
that the essential features of a large class of HOTIs are captured by
topological multipolar response theories. In this work, we show that these
multipolar responses can be realized in interacting lattice models, which
conserve both charge and dipole. In this work we study several models in both
the strongly interacting and mean-field limits. In $2$D we consider a ring-
exchange model which exhibits a quadrupole response, and can be tuned to a
$C_{4}$ symmetric higher order topological phase with half-integer quadrupole
moment, as well as half-integer corner charges. We then extend this model to
develop an analytic description of adiabatic dipole pumping in an interacting
lattice model. The quadrupole moment changes during this pumping process, and
if the process is periodic, we show the total change in the quadrupole moment
is quantized as an integer. We also consider two interacting $3$D lattice
models with chiral hinge modes. We show that the chiral hinge modes are
heralds of a recently proposed “dipolar Chern-Simons” response, which is
related to the quadrupole response by dimensional reduction. Interestingly, we
find that in the mean field limit, both the $2$D and $3$D interacting models
we consider here are equivalent to known models of non-interacting HOTIs (or
boundary obstructed versions). The self-consistent mean-field theory treatment
provides insight into the connection between free-fermion (mean-field)
theories having vanishing polarization and interacting models where dipole
moments are microscopically conserved.

###### Contents

  1. I Introduction
  2. II Linear Response Theory of Polarization
  3. III Quadrupole Moment in Dipole Conserving Systems
    1. III.1 Exactly Solvable Dipole Pumping Process
    2. III.2 Continuous Dipole Pumping and Mean Field Theory
  4. IV Dipolar Chern-Simons Response Action
  5. V $C_{4}\mathcal{T}$ Symmetric Lattice model
    1. V.1 Mean Field Analysis
    2. V.2 Continuum Analysis
    3. V.3 Effective Response Theory
  6. VI $M\mathcal{T}$ Symmetric Lattice model
    1. VI.1 Mean Field Analysis
    2. VI.2 Dimensional Reduction
  7. VII Conclusion and Outlook
  8. A Details of the Self-Consistent Mean Field Theory Calculations
  9. B Quadrupole Moment of the Mean Field Hamiltonian
  10. C Role of the Phase Fields in the Mean Field Limit
  11. D Continuum Analysis of $M\mathcal{T}$ model

##  I Introduction

Through intense theoretical effort, symmetry protected topological phases
(SPTs) have become a well understood area of condensed matter physicsHasan and
Kane (2010); Senthil (2015); Kane and Mele (2005); Bernevig  _et al._ (2006);
Kitaev (2009); Ryu  _et al._ (2010); Hasan and Moore (2011); Chen  _et al._
(2013); Kapustin  _et al._ (2015). The most notable feature of these phases of
matter is that they have a gapped bulk and robust surface states. As the name
implies, the surface states of an SPT are protected by a symmetry of the
system, and in the absence of this symmetry the system can be smoothly
deformed into a trivial insulator without closing the bulk gap. For example,
in $1$D there exists an SPT protected by inversion symmetry with quantized
half-integer surface chargesSu  _et al._ (1979); Jackiw and Schrieffer (1981),
and in $3$D there exists an SPT protected by time reversal symmetry that hosts
a single surface Dirac coneFu  _et al._ (2007); Moore and Balents (2007).
Additionally, there are also topological phases, where the surface states
remain robust even in the absence of any symmetry. In $2$D such a topological
phase is realized in integer quantum hall insulatorsKlitzing  _et al._ (1980).

In addition to their gapless surfaces states, the topological nature of SPTs
is also manifest through various quantized topological response phenomena. For
example, the aforementioned $1$D SPT has a half-integer quantized polarization
response, the $2$D topological phase has a quantized Chern-Simons response,
and the $3$D SPT has quantized axion electrodynamics King-Smith and Vanderbilt
(1993a); Qi  _et al._ (2008, 2013); Zhang  _et al._ (1989). Additionally,
certain topological responses are also related to each other through a
dimensional hierarchyQi  _et al._ (2008). For example, the $2$D Chern-
Simons/quantum Hall response on a thin torus can be mapped on to the charge
polarization response of a $1$D insulator. In this limit, the Laughlin gauge
argument for Hall current in $2$DLaughlin (1981) maps onto Thouless charge
pumping in $1$DThouless (1983). In SPTs, topological responses such as these
are quantized by the interplay of topology and symmetry, and are robust,
provided the relevant symmetry is preserved, and the bulk of the system
remains gapped. Because of their sharp quantization, topological responses can
be used to predict experimentally relevant characteristics of topological
materials.

While the study of SPTs initially focused on systems protected by internal
symmetries, it has since evolved to show that there are a rich array of
topological systems that are protected by spatial symmetries–referred to as
topological crystalline insulatorsFu (2011); Hughes  _et al._ (2011); Hsieh
_et al._ (2012); Isobe and Fu (2015); Cheng  _et al._ (2016); Bradlyn  _et
al._ (2017). Of particular interest for this work are topological crystalline
insulators having _gapped_ boundaries, but which exhibit gapless corner, or
hinge modes. Systems with these novel surface states are known as higher order
topological insulators (HOTIs)Benalcazar  _et al._ (2017a, b); Langbehn  _et
al._ (2017); Song  _et al._ (2017a); Serra-Garcia  _et al._ (2018); Peterson
_et al._ (2018); Imhof  _et al._ (2018). The first example was provided in in
Ref. Benalcazar  _et al._ , 2017a, where it was shown that there exists a non-
interacting HOTI in $2$D that can be protected by $C_{4}$ spatial rotations
symmetry (or $M_{x}$ and $M_{y}$ mirror symmetries for a boundary obstructed
phaseBenalcazar  _et al._ (2017a); Khalaf  _et al._ (2019)). This model
exhibits a symmetry-quantized quadrupole moment, and when defined a on a
lattice with boundary this HOTI has half-integer charges localized at its
corners. Similar to how a polarized $1$D system can be related to a $2$D Chern
insulator by dimensional reduction, the $2$D quadrupole insulator is related
to a $3$D HOTI with chiral hinge modesBenalcazar  _et al._ (2017b). Due to the
rich phenomena exhibited by the quadrupole insulator and its cousins, HOTIs
have attracted a great deal of theoretical attention. In particular, there
have been several works classifying these systems according to their
symmetriesSong  _et al._ (2017b); Song and Schnyder (2017); Rasmussen and Lu
(2018); Thorngren and Else (2018). Despite these advances, the topological
response of HOTIs still remains largely nebulous. This is in striking contrast
to SPTs with internal symmetries, and topologically ordered systems, which
have well understood topological responses.

Recently, it has been proposed that the topological responses of certain HOTIs
can be written in terms of symmetric rank-2 gauge fieldsYou  _et al._ (2019).
Rank-2 gauge fields are symmetric two index gauge fields $A_{ij}$ that, for
scalar-charge theories, transform as $A_{ij}\rightarrow
A_{ij}+\partial_{i}\partial_{j}\Lambda$ under a gauge transformation
$\Lambda$. In the last several years, rank-2 gauge fields have attracted
attention in condensed matter physics because of their relation to fracton
phasesNandkishore and Hermele (2019a); Pretko  _et al._ (2020). A key feature
of these gauge theories is that the Gauss’s law for rank-2 gauge fields
naturally conserves both charge and dipole Pretko (2017a, b, 2018). In
particular, this inherent dipole conservation leads to modified dynamics for
the excitations in fracton phasesNandkishore and Hermele (2019b); Haah (2011);
Vijay  _et al._ (2016); Hsieh and Halász (2017); Pretko (2017c); Ma  _et al._
(2018); Gromov (2019). Indeed, the exotic conservation laws can force
quasiparticles to be either immobile, or confined to move along sub-
dimensional manifolds, such as lines or planes in $3$D. In lattice models,
this sub-dimensional confinement can be interpreted as arising from a
microscopic dipole conservation lawXu (2006); Xu and Hořava (2010); Pretko
(2017c). Because of this, it has been proposed that the essential physics of
fracton systems can be captured by effective rank-2 gauge theories. Generally
it has been shown that rank-2 gauge fields can couple to matter theories that
conserve dipolePretko and Radzihovsky (2018). In these systems, the dipole
moment can be treated as the conserved charge of a global $1$-form
symmetrySeiberg (2019); Seiberg and Shao (2020); Dubinkin  _et al._ (2020).

Returning to HOTIs, it was shown in Ref. You  _et al._ , 2019, that a rank-2
quadrupole response can describe the quantized corner charges of a $2$D HOTI.
This rank-2 quadrupole response can be considered as a rank-2 analog of the
(rank-1) charge polarization of a $1$D systemQi and Zhang (2011). Similar to
how charge polarization describes the surface charge of a system, the rank-2
quadrupole response describes the boundary polarization and corner charge of a
system. The rank-2 quadrupole response is also related to a topological
response in $3$D by dimensional reduction. This $3$D response is a kind of
dipolar Chern-Simons responseYou  _et al._ (2019), and it can describe a $3$D
HOTI with chiral hinge modes. As noted before, rank-2 gauge fields can couple
only to matter theories that conserve a dipole moment. In order to conserve
dipole, such a theory must not have any single-particle charge dynamics, since
single-particle charge transport changes the dipole moment of a system.
However, it is possible to have multi-particle charge dynamics in a dipole
conserving system, e.g., correlated/pair hopping terms which are represented
by quartic interaction terms. Because of this, the multipolar rank-2
topological responses of Ref. You  _et al._ , 2019 most naturally describe
interacting HOTIs.

In this work, our goal is to develop tools to further investigate these
interacting HOTIs and their topological responses, and then to apply our
techniques to several dipole conserving lattice models of interacting
fermions. To analyze these systems, we formulate a linear response theory of
the rank-2 quadrupole moment of dipole conserving solids. This formulation of
the quadrupole moment is largely analogous to the linear response formulation
of polarization in charge conserving solidsResta and Vanderbilt (2007).
Similar to how a change in polarization is viewed as a pumping of charge
across the system, here the change in the quadrupole moment is viewed as a
pumping of dipoles across the system. For a system with boundaries, the change
in quadrupole moment can also lead to surface charge currents and a change in
corner charges. Using our linear response formalism, we are able to exactly
calculate the quadrupole response of the interacting quadrupole model of Ref.
You  _et al._ , 2019 and provide an alternative confirmation that it is a HOTI
having a quantized quadrupole moment. Interestingly, we are also provide a
protocol for tuning the parameters of the model to fully demonstrate a dipole
pumping process in this strongly interacting system, i.e., we can see how
dipole is pumped across the system as the quadrupole moment of the system is
tuned by varying parameters of the model. Furthermore, we show that analyzing
the model within a self-consistent mean field theory approximation, maps this
interacting HOTI to the non-interacting quadrupole model of Ref. Benalcazar
_et al._ , 2017a. We can use this mean-field approximation to demonstrate a
remarkable connection between non-interacting HOTIs and topological rank-2
responses, the latter of which would naively apply only to systems with exact
dipole conservation.

We then move on to consider two interacting $3$D models. The first model we
study is related to the dipole conserving $2$D interacting quadrupole model by
dimensional reduction. This model conserves the dipole moment along the $x$
and $y$-directions, and is invariant under $C_{4}\mathcal{T}$ symmetry (the
product of $C_{4}$ symmetry and time reversal symmetry $\mathcal{T}$). Since
the $3$D dipole Chern-Simons response is related to the $2$D rank-2 quadrupole
response via dimensional reductionYou  _et al._ (2019), we expect this $3$D
model to exhibit a quantized dipole Chern-Simons response and have chiral
hinge modes. Indeed, we are able to show that the bulk and surface responses
of this $3$D model are identical to those predicted by the dipole Chern-Simons
theory. We also verify that this model supports chiral hinge modes that are
consistent with $C_{4}\mathcal{T}$ symmetry. Similar to 2D case, we show that
within a self-consistent mean field approximation this interacting $3$D HOTI
is equivalent to a known non-interacting chiral hinge insulator in $3$D. The
second $3$D model we consider is a related interacting model that is invariant
under $M_{x}\mathcal{T}$ (the product of $M_{x}$ mirror symmetry and time
reversal symmetry $\mathcal{T}$) and $M_{y}\mathcal{T}$ (the product of
$M_{y}$ mirror symmetry and time reversal symmetry $\mathcal{T}$) instead of
$C_{4}\mathcal{T}$. This model also has chiral hinge modes, however, in the
mean field limit, this model is equivalent to a layered system that harbors
dangling Chern insulators on the top and bottom surfaces having Chern number
$+1$ and $-1$ respectively. Because of this, the mean field model constitutes
a “boundary-obstructed topological insulator” using the terminology of Ref.
Khalaf  _et al._ , 2019.

In summary, our work on the generic formalism and the explicit $2$D and $3$D
models represents a further link between higher order topological phases
protected by spatial symmetries and analogous fracton-like systems having
stricter multipole conservation laws. Our paper is organized as follows. In
Section II we motivate our linear response formulation of the quadrupole
moment by reviewing the linear response theory of polarization, and
considering relevant examples. In Section III we formulate a linear response
definition of the rank-2 quadrupole response in solids. We use this formalism
to calculate the quadrupole response of a dipole conserving $2$D lattice model
and show how adiabatic deformations of the model can lead to a change in the
quadrupole moment via dipole pumping. We also show that this lattice model can
be tuned to a $C_{4}$ symmetric topological phase, with half-integer
quadrupole moment. In Section IV we review the $3$D dipolar Chern-Simons
response theory, which is related to the rank-2 quadrupole response in $2$D by
dimensional reduction. In Section V we present a $3$D $C_{4}\mathcal{T}$
symmetric lattice model that is related to the dipole conserving $2$D lattice
model we discussed earlier by dimensional reduction. We show that this model
realizes the dipole Chern-Simons response and has protected chiral hinge
modes. In Sec. VI we present a related $\mathcal{M}_{y}\mathcal{T}$ and
$\mathcal{M}_{y}\mathcal{T}$ symmetric lattice model in $3$D. We show that
this model also has protected hinge modes, and that it is related to a
boundary obstructed phase in the mean-field limit. We conclude our results in
Section VII. We also have several appendices that contain the technical
details of our calculations.

##  II Linear Response Theory of Polarization

Before considering the linear response theory of the quadrupole moment in
solids, it will be useful to first review the modern linear response theory of
polarizationResta and Vanderbilt (2007) (from here on we will use the terms
“dipole moment” and polarization interchangeably in this paper). For
simplicity we shall focus on systems in $1$D. Our starting point is the
effective response action for a system with polarization $P$:

$S_{P}=\int dtdxP_{x}(\partial_{x}A_{0}-\partial_{t}A_{x}).$ (1)

Here, the gauge fields $A_{x}$ and $A_{0}$ are background fields that serve as
probes for the charge responses of the system. To show that Eq. 1 does in fact
describe the polarization of a $1D$ system, let us consider an infinite line
where $P_{x}=P_{0}$ for $0\leq x\leq L_{x},$ and $P_{x}=0$ everywhere else.
This corresponds to a finite $1$D system with polarization $P_{0}$ embedded in
an unpolarized vacuum. According to Eq. 1, the electric charge is
$j_{0}=\frac{\delta}{\delta A_{0}}S_{P}=-\partial_{x}P_{x}$. This means that a
charge of $-P_{0}$ will be localized at the $x=0$ boundary, and a charge of
$+P_{0}$ localized at $x=L_{x}$ boundary (see Fig. 1). This is exactly the
boundary charge distribution we expect for a $1$D system with polarization
$P_{0}$.

Figure 1: The bound charges of a $1$D sample with polarization $P_{x}=P_{0}$
(gray) embedded in a $1$D vacuum with $P_{x}=0$ (white).

An important feature of Eq. 1 is that $\partial_{x}A_{t}-\partial_{t}A_{x}$ is
a total derivative. Because of this, there are only non-trivial charge
responses when $P_{x}$ varies in space or time, e.g., when boundaries are
present. In other words, only changes in polarization are observable. This was
a key insight in the development of the modern theory of polarizationResta
(1992); King-Smith and Vanderbilt (1993b); Resta (1994). In this framework,
the “polarization” of a given system is properly defined as the polarization
of a system relative to an (unpolarized) reference state trivial insulator
(which is often taken to be a trivial atomic insulator). With this in mind,
let us consider the change in polarization during an adiabatic process.
Physically, this change in polarization can be attributed to the pumping of
charge across the system, i.e., to a charge current. If we parameterize the
adiabatic process by $\theta,$ then the change in polarization for a system
that is minimally coupled to a flat background gauge field $A_{x}$ is given by

$\begin{split}\frac{\partial}{\partial\theta}P=\lim_{\epsilon\rightarrow
0}\frac{i}{\epsilon L}\sum_{n\neq 0}&\Big{[}\frac{\bra{0}\frac{\partial
H}{\partial\theta}\ket{n}\bra{n}\frac{\partial
H}{\partial{A_{x}}}\ket{0}}{\epsilon+E_{0}-E_{n}}\\\
&-\frac{\bra{0}\frac{\partial H}{\partial{A_{x}}}\ket{n}\bra{n}\frac{\partial
H}{\partial\theta}\ket{0}}{\epsilon+E_{n}-E_{0}}\Big{]},\end{split}$ (2)

where $H$ is the Hamiltonian of the system we are considering, and $\ket{n}$
is an energy eigenstate with energy $E_{n}$. After some algebra, this
expression can be rewritten in terms of a Berry curvature as

$\frac{\partial}{\partial\theta}P_{x}=\frac{i}{L}\Big{[}\frac{\partial}{\partial
A_{x}}\bra{0}\frac{\partial}{\partial\theta}\ket{0}-\frac{\partial}{\partial\theta}\bra{0}\frac{\partial}{\partial
A_{x}}\ket{0}\Big{]},$ (3)

where $\ket{0}$ is the ground state wavefunction of the model being
considered, and we are implicitly taking $\ket{0}$ to be a function of both
the adiabatic parameter $\theta,$ as well as the background gauge field
$A_{x}$. Eq. 3 is invariant under shifting $\ket{0}$ by an arbitrary phase
that depends on $\theta$ and $A_{x}$. Because of this, we can choose an
overall phase for the ground state wavefunction such that
$\frac{\partial}{\partial
A_{x}}\bra{0}\frac{\partial}{\partial\theta}\ket{0}=0$. For this choice, the
total change in polarization during an adiabatic process is found to be

$\Delta P_{x}=\Delta\left[\frac{-i}{L}\bra{0}\frac{\partial}{\partial
A_{x}}\ket{0}\right].$ (4)

For band insulators, where the ground state can written as a product of
eigenfunctions of the single particle Hamiltonian, the change in polarization
can written as

$\Delta
P_{x}=\Delta\left[(-i)\int\frac{dk_{x}}{2\pi}\sum_{\alpha\in\text{occ}}\bra{\alpha,k_{x}}\frac{\partial}{\partial
k_{x}}\ket{\alpha,k_{x}}\right],$ (5)

where $\ket{\alpha,k}$ are the eigenfunctions of the single particle
Hamiltonian, $\alpha$ is the band index, and the sum is over the occupied
bands. This is the celebrated Berry phase formulation of the polarization of
band insulatorsKing-Smith and Vanderbilt (1993b); Resta (1994).

An important feature of this framework is that $\Delta P_{x}$ is only defined
modulo an integer. In Eq. 4 this ambiguity is due to the fact that the ground
state wavefunction can be multiplied by the gauge invariant Wilson loop
$\ket{0}\rightarrow\exp(i\int dxA_{x})\ket{0}$ (we set $e=\hbar=1$). In Eq. 5
this ambiguity is due to the fact that the eigenfunctions of the single
particle Hamiltonian can be shifted by the $k_{x}$ dependent phase
$\ket{\alpha,k_{x}}\rightarrow e^{ik_{x}}\ket{\alpha,k_{x}}$. In both cases,
the redefinition of the ground state wavefunction does not change any physical
properties of the ground state, but it will shift the polarization by $+1$.
Physically, this phase shift corresponds to moving every particle in the
system over by one unit cell. For a system with periodic boundary conditions,
such a translation is trivial. For a system with boundaries, translating each
particle over by one unit cell will add or remove a single unit of charge from
the boundaries. This effect can be canceled by adding an integer number of
electrons to the boundaries and hence does not meaningfully affect the bulk
polarization of the system.

The integer ambiguity in the polarization enables the identification of an SPT
protected by inversion symmetry. The polarization response in Eq. 1 is odd
under inversion symmetry, so in an inversion symmetric system $P_{x}=-P_{x}$.
Since $P_{x}$ is only defined modulo an integer, $P_{x}=-P_{x}$ is satisfied
when $P_{x}=0$ or $P_{x}=1/2$. The former corresponds to a trivial insulator
with no fractional boundary charges, while the latter corresponds to a $1$D
SPT with half-integer boundary charges. It is well known that this SPT can be
realized in the Su-Schrieffer–Heeger (SSH) model for polyacetyleneSu  _et al._
(1979).

Having established the linear response theory of polarization, it will be
useful to apply this formalism to a lattice model. We would also like to
verify that the polarization is only defined modulo an integer, and that the
polarization is quantized in units of $1/2$ in an inversion symmetric system.
To do this, we shall consider the SSH Hamiltonian with an additional onsite
potential term:

$\begin{split}H^{\text{SSH}}&=\sum_{k_{x}}\bm{c}^{\dagger}(k_{x})h^{\text{SSH}}(k_{x})\bm{c}(k_{x}),\\\
h^{\text{SSH}}(k_{z})&=(u+v\cos(k_{z}))\tau_{1}+v\sin(k_{z})\tau_{2}+\mu\tau_{3},\end{split}$
(6)

where $\bm{c}=(c_{1},c_{2})$, and $\tau_{i}$ are the Pauli matrices. Here, $u$
and $v$ are the amplitudes for intracell and intercell hopping respectively,
and $\mu$ is the strength of the onsite potential. When $\mu=0$, the model has
inversion symmetry. Here, we are interested in the change in polarization
during an adiabatic process, and so we will make the parameters $u$, $v$ and
$\mu$ functions of an adiabatic parameter $\theta$. Specifically, we will set

$\begin{split}&u=\max(\cos(\theta),0),\\\ &v=\max(-\cos(\theta),0),\\\
&\mu=-\sin(\theta).\end{split}$ (7)

This process is periodic with respect to $\theta$, and the Hamiltonian at
$\theta=0$ is the same as the Hamiltonian at $\theta=2\pi$. Additionally, the
system has inversion symmetry when $\theta=0(2\pi)$ and $\pi$. Using Eq. 5 the
polarization of this system is given by

$\Delta P_{x}=\begin{cases}0&\text{for }0\leq\theta\leq\pi/2\\\
\frac{1}{2}[1-\sin(\theta)]&\text{for }\pi/2\leq\theta\leq 3\pi/2\\\
1&\text{for }3\pi/2\leq\theta\leq 2\pi\end{cases},$ (8)

where $\Delta P_{x}\equiv P_{x}(\theta)-P_{x}(0)$. The polarization as a
function of $\theta$ during this process is shown in Fig. 2. From this
calculation we can clearly see that when $\theta$ is increased from $0$ to
$2\pi$, the polarization of the system increases from $0$ to $1$. As noted
before, this process is periodic, confirming that the polarization is only
defined modulo an integer. Additionally, the model we are considering has
inversion symmetry when $\theta=0(2\pi)$ and $\theta=\pi$, and at these points
the polarization is $0(1)$ and $1/2$ respectively. At $\theta=0(2\pi)$ the
system only has intracell terms, and hence is a trivial inversion symmetric
insulator. At $\theta=\pi$ the we therefore expect that the system is in the
SPT phase with polarization $P_{x}=1/2$. In this phase, the model has half-
integer charges localized at it boundaries.

Figure 2: The change in polarization $\Delta P_{x}\equiv
P_{x}(\theta)-P_{x}(0)$ of Eq. 6 as a function of the adiabatic parameter
$\theta$ in Eq. 7.

It is also useful to consider the fully continuous parameterization of the SSH
chain (as opposed to the piecewise continuous parameterization in Eq. 7),

$\begin{split}&u=\gamma[1+\cos(\theta)],\\\ &v=v\text{(const.)},\\\
&\mu=-\gamma\sin(\theta).\end{split}$ (9)

For $v/\gamma<2$, this pumping process is topologically equivalent to the one
given in Eq. 7. In particular, at $\theta=\pi$, the parameterization in Eq. 9
leads to the same inversion symmetric SPT discussed earlier, and as $\theta$
is increased from $0$ to $2\pi$, the polarization $P_{x}$ increases by $+1$.
Additionally, since the parameterization depends smoothly on $\theta$, we can
relate the $1$D model to a $2$D model by identifying $\theta$ with the lattice
momentum $k_{y}$. This $2$D model is an insulator with Chern number $+1$. For
a more detailed discussion of this mapping see Ref. Qi  _et al._ , 2008.

##  III Quadrupole Moment in Dipole Conserving Systems

Having reviewed the linear response formulation of polarization, we can now
turn our attention to determining the rank-2 quadrupole response of a system.
This formulation will largely parallel the formulation of polarization we
reviewed in Sec. II. Here, we shall start by considering the quadrupole
response term. For a $2D$ system coupled to a background rank-2 gauge field,
the $Q_{xy}$ quadrupole response term is given byYou  _et al._ (2019):

$S_{Q}=\int
dtd\bm{r}Q_{xy}[\partial_{x}\partial_{y}A_{0}-\partial_{t}A_{xy}].$ (10)

This term can naturally be interpreted as the rank-2 generalization of the
polarization response given in Eq. 1. In $2$D, there are also $xx$ and $yy$
quadrupole responses, which can be written analogously to Eq. 10. Here, we
shall focus on the $xy$ quadrupole response. It is straightforward to
generalize the results of this section to other quadrupole responses.

To show that Eq. 10 does indeed give the desired physics of a model with non-
vanishing quadrupole moment, let us consider an infinite 2D plane, where
$Q_{xy}=Q_{0}$ for $0\leq x\leq L_{x}$ and $0\leq y\leq L_{y},$ and $Q_{xy}=0$
everywhere else. This corresponds to a rectangular system with quadrupole
moment $Q_{0}$ embedded in a $2$D vacuum. Using $j_{0}=\frac{\delta}{\delta
A_{0}}S_{Q}=\partial_{x}\partial_{y}Q_{xy}$, we find that this system has
charge $+Q_{0}$ localized at $(x,y)=(0,0)$ and $(L_{x},L_{y})$, and charge
$-Q_{0}$ localized at $(x,y)=(0,L_{y})$ and $(L_{x},0)$ (see Fig. 3). This is
exactly the corner charge distribution we expect for a system with quadrupole
moment $Q_{0}$.

Figure 3: The bound charges of an open boundary $2$D sample with quadrupole
moment $Q_{xy}=Q_{0}$ (gray) embedded in a $2$D vacuum with $Q_{xy}=0$
(white).

Since $\partial_{x}\partial_{y}A_{0}-\partial_{t}A_{xy}$ is a total
derivative, Eq. 10 is only non-trivial when $Q_{xy}$ varies in space or time,
e.g., when boundaries are present. Because of this, only changes in the
quadrupole moment are observable. In this framework, the “quadrupole moment”
of a given system is properly defined as the quadrupole moment of a system
relative to a reference state. Based on this, we will consider the change in
the quadrupole moment during an adiabatic process. Analogously to how a change
in polarization is caused by pumping charges, here the change in the
quadrupole moment can be attributed to pumping dipoles. If we parameterize the
adiabatic process by $\theta$, then the change in quadrupole moment $Q_{xy}$
for a system that is coupled to a flat background gauge field $A_{xy}$ is
given by

$\begin{split}\frac{\partial}{\partial\theta}Q_{xy}=\lim_{\epsilon\rightarrow
0}\frac{i}{\epsilon L_{x}L_{y}}\sum_{n\neq
0}&\left[\frac{\bra{0}\frac{\partial
H}{\partial\theta}\ket{n}\bra{n}\frac{\partial
H}{\partial{A_{xy}}}\ket{0}}{\epsilon+E_{0}-E_{n}}\right.\\\
&-\left.\frac{\bra{0}\frac{\partial
H}{\partial{A_{xy}}}\ket{n}\bra{n}\frac{\partial
H}{\partial\theta}\ket{0}}{\epsilon+E_{n}-E_{0}}\right].\end{split}$ (11)

As before, this expression can be simplified and written in terms of a Berry
curvature as

$\partial_{\theta}Q_{xy}=\frac{i}{L_{x}L_{y}}[\partial_{A_{xy}}(\bra{0}\partial_{\theta}\ket{0})-\partial_{\theta}(\bra{0}\partial_{A_{xy}}\ket{0})],$
(12)

where $\ket{0}$ is the ground state of the microscopic model we are
considering, and $L_{i}$ is the length of the system in the $i$-direction.
Here, we are implicitly treating the ground state wavefunction as a function
of both the adiabatic parameter $\theta$ and the rank-2 gauge field $A_{xy}$.
Eq. 12 is invariant under shifting the wavefunction by an arbitrary phase that
depends on $\theta$ and $A_{xy}$. Because of this, we can choose a overall
phase for the ground state wavefunction such that
$\partial_{A_{xy}}(\bra{0}\partial_{\theta}\ket{0})=0$. For this choice of
phase Eq. 12 can reduce to

$\Delta
Q_{xy}=\Delta\left[\frac{-i}{L_{x}L_{y}}\bra{0}\partial_{A_{xy}}\ket{0}\right].$
(13)

A similar expression for the quadrupole moment was presented in Ref. Dubinkin
_et al._ , 2019.

Based on Eq. 13, we can infer that the quadrupole moment $Q_{xy}$ is only
defined modulo $1$. This is because the ground state wavefunction can be
multiplied by the gauge invariant “Wilson surface”
$\ket{0}\rightarrow\exp(i\int dxdyA_{xy})\ket{0}$. Shifting the ground state
wavefunction by such a term will increase the quadrupole moment $Q_{xy}$ by
$+1$. Physically, this will add/remove an integer amount of charge from the
corners of the system, while leaving the bulk of the system unchanged. An
integer ambiguity of this form was also seen when considering the polarization
in Sec. II.

Due to this integer ambiguity, we can predict the existence of a HOTI with
half-integer quadrupole moment protected by $C_{4}$ rotation symmetry. Since
the quadrupole response in Eq. 10 is odd under $C_{4}$ rotations, a $C_{4}$
invariant insulator must have an $xy$ quadrupole moment satisfying
$Q_{xy}=-Q_{xy}$. Since $Q_{xy}$ is only defined modulo an integer
$Q_{xy}=-Q_{xy}$ is satisfied by $Q_{xy}=0$ and $Q_{xy}=1/2$. The former is a
trivial $C_{4}$ symmetric insulator, while the latter is a HOTI with quantized
half-integer quadrupole moment. This HOTI will have half-integer corner
charges, similar to those found in the non-interacting quadrupole insulator of
Ref. Benalcazar  _et al._ , 2017a. This logic also predicts a similar
quantized quadrupole insulator with half-integer quadrupole moment protected
by $M_{x}$ mirror symmetry and $M_{y}$ mirror symmetry, since the $xy$
quadrupole response is odd under both of these symmetries.

Having established a linear response formalism, we can now turn our attention
to calculating the change in the quadrupole moment in a microscopic lattice
model. Based on the rank-2 quadrupole response action in Eq. 10, we are
interested in lattice models that couple to the scalar potential $A_{0}$, and
the rank-2 gauge field $A_{xy}$. In order for a lattice model to couple to
these gauge fields it must conserve both charge and dipole in the $x$ and
$y$-directions. For a system of lattice fermions, global charge conservation
corresponds to the symmetry that sends $c(\bm{r})\rightarrow
c(\bm{r})e^{i\alpha}$, where $c(\bm{r})$ is the lattice fermion annihilation
operator, and $\alpha$ is a constant. Similarly, dipole conservation in the
$x$ and $y$-directions corresponds to the symmetry that sends
$c(\bm{r})\rightarrow c(\bm{r})e^{i\bm{\beta}\cdot\bm{r}}$, where $\bm{\beta}$
is a constant two component vector.

To show why these symmetries are necessary, let us consider an arbitrary
lattice model that is composed of fermion operators $c(\bm{r}),$ and the
background gauge fields $A_{0}$ and $A_{xy}.$ Under gauge transformations
$\Lambda({\bf{r}})$ these fields transform as

$\begin{split}&A_{0}(\bm{r})\rightarrow
A_{0}(\bm{r})+\partial_{t}\Lambda(\bm{r}),\\\ &A_{xy}(\bm{r})\rightarrow
A_{xy}(\bm{r})+\Lambda(\bm{r})-\Lambda(\bm{r}+\hat{x})\\\
&\phantom{======}-\Lambda(\bm{r}+\hat{y})+\Lambda(\bm{r}+\hat{x}+\hat{y})\\\
&\phantom{A_{xy}(\bm{r})}\equiv
A_{xy}(\bm{r})+\Delta_{x}\Delta_{y}\Lambda(\bm{r}),\\\ &c(\bm{r})\rightarrow
c(\bm{r})e^{i\Lambda(\bm{r})},\end{split}$ (14)

where $\Delta_{i}$ is the lattice derivative in the $i$-direction, and we have
suppressed any dependence on $t$. Let us now consider a gauge transformation
of the form $\Lambda=\alpha$ (const.). The gauge fields $A_{0}$ and $A_{xy}$
are invariant under such transformations, while the fermions transform as
$c(\bm{r})\rightarrow c(\bm{r})e^{i\alpha}$. So in order for the system to be
gauge invariant, it must be invariant under shifting the phase of the fermions
by a constant amount. Similarly, we can also consider a gauge transformation
of the form $\Lambda(\bm{r})=\bm{\beta}\cdot\bm{r}$. Under this gauge
transformation both $A_{0}$ and $A_{xy}$ are invariant, while the fermions
transform as $c(\bm{r})\rightarrow c(\bm{r})e^{i\bm{\beta}\cdot\bm{r}}$. So in
order for the system to be gauge invariant, it must also be invariant under
shifting the phase of the fermions by an amount that depends linearly on
position.

While charge conservation is fairly common in lattice models, dipole
conservation is more unusual, and places strong constraints on the types of
terms that can appear in a lattice Hamiltonian. Importantly, because of dipole
conservation, single-particle hopping terms such as
$c^{\dagger}(\bm{r}+\hat{x})c(\bm{r})$ and
$c^{\dagger}(\bm{r}+\hat{y})c(\bm{r})$ are not allowed. However if a system
has multiple degrees of freedom within a unit cell, single-particle intracell
terms are allowed, since they do not change the dipole moment of the system.
This means one can include terms like $c_{i}^{\dagger}(\bm{r})c_{j}(\bm{r})$
where $i$ and $j$ label the different fermionic degrees of freedom within a
given unit cell.

Although single-particle hopping terms are not allowed, quartic interactions
can allow for pairs of electrons to have dynamics. A simple term of this form
is the ring exchange term:
$c^{\dagger}(\bm{r})c^{\dagger}(\bm{r}+\hat{y}+\hat{x})c(\bm{r}+\hat{x})c(\bm{r}+\hat{y})$.
A quick calculation confirms that this term is indeed invariant under linear
phase shifts, and does not change the dipole moment in the $x$ or
$y$-direction. Physically, this term can be thought of as a dipole hopping
term. To see this, we note that
$c^{\dagger}(\bm{r}+\hat{y}+\hat{x})c(\bm{r}+\hat{y})$ can be interpreted as
creating a dipole with dipole vector $\hat{x}$ centered at
$\bm{r}+\hat{y}+\tfrac{\hat{x}}{2}$ (the inclusion of $\tfrac{\hat{x}}{2}$
indicates that the dipole is defined on the link between
$\bm{r}+\hat{y}+\hat{x}$ and $\bm{r}+\hat{y}$). Similarly,
$c(\bm{r}+\hat{x})c^{\dagger}(\bm{r})$ creates a dipole with dipole vector
$-\hat{x}$ (equiv. annihilates a dipole with dipole vector $\hat{x}$) centered
at $\bm{r}+\tfrac{\hat{x}}{2}$. The ring exchange term thereby hops an
$\hat{x}$-oriented dipole from $\bm{r}+\tfrac{\hat{x}}{2}$ to
$\bm{r}+\hat{y}+\tfrac{\hat{x}}{2}$. This process can also be interpreted as
hopping a $\hat{y}$-oriented dipole one unit in the $\hat{x}$-direction, from
$\bm{r}+\frac{\hat{y}}{2}$ to $\bm{r}+\hat{x}+\frac{\hat{y}}{2}$. The ring
exchange term minimally couples to the rank-2 gauge field $A_{xy}$ via a
rank-2 Peierls factor of the form
$c^{\dagger}(\bm{r})c^{\dagger}(\bm{r}+\hat{y}+\hat{x})c(\bm{r}+\hat{x})c(\bm{r}+\hat{y})e^{iA_{xy}(\bm{r})}$You
_et al._ (2019); Dubinkin  _et al._ (2019). With the rank-2 Peierls factor
included, this term is invariant under the gauge transformations given in Eq.
14.

Based on these considerations, we introduce the following $2D$ dipole
conserving lattice model with 4 fermionic degrees of freedom per unit cell:

$\begin{split}H^{Q}&=\sum_{\bm{r}}\bm{c}^{\dagger}(\bm{r})h^{o}\bm{c}(\bm{r})-A_{0}(\bm{r})\bm{c}^{\dagger}(\bm{r})\bm{c}(\bm{r})\\\
&\phantom{=}-Vc_{1}^{\dagger}(\bm{r})c^{\dagger}_{2}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})c_{4}(\bm{r}+\hat{y})e^{iA_{xy}(\bm{r})}\\\
&\phantom{=}+h.c.,\\\
h^{o}&=\mu\Gamma_{0}+t(\Gamma_{2}+\Gamma_{4}).\end{split}$ (15)

Here, $c_{i}({\bf{r}})$ ($i=1...4$) are the four lattice fermion operators for
a unit cell $\bm{r}=(x,y)$ in the $2$D square lattice (see Fig. 4). The
$\Gamma$ matrices are defined as $\Gamma_{0}=\tau_{3}\otimes\tau_{0}$,
$\Gamma_{k}=-\tau_{2}\otimes\tau_{k}$, and
$\Gamma_{4}=\tau_{1}\otimes\tau_{0},$ for $k=1,2,3$, where $\tau_{1,2,3}$ are
the Pauli matrices. The $4\times 4$ matrix $h^{o}$ contains the single-
particle intracell terms, and $V$ is the amplitude of the ring exchange term.
Since the square lattice is bipartite, the sign of $V$ can be changed by an
appropriate unitary transformation. Because of this, we will take $V>0$
without loss of generality. Within the intracell term $h^{o}$, $t$ is the
amplitude for intracell hopping, while $\mu$ is a staggered onsite potential.
We have also included the coupling to the rank-2 gauge field $A_{xy},$ as well
as the scalar potential $A_{0}$. We will restrict our attention to the case
where the model is half filled ($2$ fermions per unit cell).

Figure 4: Schematic of the dipole conserving lattice model Eq. 15. Orange
squares correspond to the ring exchange interactions with amplitude $V$. Solid
(dashed) blue lines correspond to intracell couplings with amplitude $t$
($-t$). Solid (dashed) green circles correspond to onsite potential with
strength $\mu$ ($-\mu$).

As desired, this model is gauge invariant, and invariant under constant and
linear phase shifts of the fermion operators $c_{i}({\bf{r}})$. When $\mu=0$,
this model also has $C_{4}$ rotation symmetry. This symmetry acts on the
internal fermionic degrees of freedom as $\bm{c}\rightarrow
U_{\mathcal{C}_{4}}\bm{c}$ in Eq. 15, where

$\begin{split}U_{\mathcal{C}_{4}}&=\begin{bmatrix}0&0&1&0\\\ 0&0&0&1\\\
0&-1&0&0\\\ 1&0&0&0\\\ \end{bmatrix}.\end{split}$ (16)

When $\mu=0$, this model also has $M_{x}$ and $M_{y}$ mirror symmetries. The
$M_{x}$ and $M_{y}$ symmetries act on the internal fermionic degrees of
freedom as $\bm{c}\rightarrow U_{M_{x}}\bm{c}$ and $\bm{c}\rightarrow
U_{M_{y}}\bm{c}$ respectively, where $U_{M_{x}}=\tau_{1}\otimes\tau_{3}$ and
$U_{M_{y}}=\tau_{1}\otimes\tau_{1}$ ($\tau_{i}$ are the Pauli matrices). When
$\mu\neq 0$ The staggered onsite potential explicitly breaks the $C_{4}$ and
mirror symmetries. The model we consider here also has a subsystem symmetry
that acts on the fermionic degrees of freedom as
$\bm{c}(\bm{r})\rightarrow\bm{c}(\bm{r})e^{if_{x}(x)+if_{y}(t)}$, where
$f_{x}(x)$ is an arbitrary function of the $x$-coordinate only, and $f_{y}(y)$
is an arbitrary function of the $y$-coordinate only. This symmetry correspond
to charge being conserved along every row and column of the $2$D lattice.

As we shall show, by adiabatically changing the parameters in this model ($t$,
$V$ and $\mu$), it is possible to pump dipole and change the quadrupole
moment. In the following subsections, we will present two related
parameterizations of this pumping process. First, we shall present a piecewise
continuous periodic parameterization of Eq. 15, where the quadrupole moment
can be found exactly. Second, we shall present a fully continuous periodic
parameterization, where the quadrupole moment can be found via self-consistent
mean field theory.

###  III.1 Exactly Solvable Dipole Pumping Process

A simple and illustrative example of dipole pumping in the lattice model Eq.
15 is found by considering the parameterization

$\begin{split}&t(\theta)=\max(\cos(\theta),0),\\\
&V(\theta)=2\sqrt{2}\max(-\cos(\theta),0),\\\
&\mu(\theta)=\sqrt{2}\sin(\theta).\end{split}$ (17)

The values of $t$, $V$, and $\mu$, for this adiabatic evolution, are shown in
Fig 5. This system has $C_{4}$ symmetry when $\theta=0$ and $\theta=\pi$. At
$\theta=0$ the only non-zero terms are the intracell hopping terms $t$, and at
$\theta=\pi$, the only non-zero terms are the ring exchange terms $V$. Based
on this we can identify the system at $\theta=0$ is a trivial $C_{4}$
symmetric insulator. We also expect that at $\theta=\pi$, the system is a HOTI
protected by $C_{4}$ symmetry, with half-integer quadrupole moment (relative
to the trivial insulator at $\theta=0$).

Figure 5: The evolution of $t$, $V$, and $\mu$ in Eq. 15 as a function of
$\theta$ using the parameterization in Eq. 17

To show that this expectation is correct, we will use Eq. 13 to calculate the
quadrupole moment. The ground state of Eq. 15 with the parameterization given
in Eq. 17 can be found exactly in three steps. First, for
$0\leq\theta\leq\pi/2$ the amplitude of the ring exchange term $V$ vanishes,
and the Hamiltonian consists only of the intracell terms $t$ and $\mu$. The
ground state for this range of $\theta$ can thereby be expressed as a of
product single-particle wavefunctions defined on a single site. For
$0\leq\theta\leq\pi/2$ the ground state is

$\begin{split}\ket{0}=&\prod_{\bm{r}}\Big{[}a_{1}(\theta)c^{\dagger}_{1}(\bm{r})+a_{2}(\theta)c^{\dagger}_{2}(\bm{r})+a_{3}(\theta)c^{\dagger}_{3}(\bm{r})\Big{]}\\\
&\times\Big{[}a_{1}(\theta)c^{\dagger}_{1}(\bm{r})-a_{2}(\theta)c^{\dagger}_{2}(\bm{r})+a_{4}(\theta)c^{\dagger}_{4}(\bm{r})\Big{]}\ket{vac},\\\
&a_{1}(\theta)=a_{2}(\theta)=\sqrt{\frac{1-\sin(\theta)}{4}},\\\
&a_{3}(\theta)=a_{4}(\theta)=-\sqrt{\frac{1+\sin(\theta)}{2}}.\end{split}$
(18)

Second, for $\pi/2\leq\theta\leq 3\pi/2$, the amplitude of the intracell
hopping $t$ vanishes, and only the ring exchange term $V$ and on-site
potential term $\mu$ remain. In this case, the different plaquettes spanned by
$c_{1}(\bm{r})$, $c_{2}(\bm{r}+\hat{x}+\hat{y})$, $c_{3}(\bm{r}+\hat{x})$, and
$c_{4}(\bm{r}+\hat{y})$ decouple from each other. The resulting 4 site
interacting problem can be solved using exact diagonalization, and the ground
state of the system can be expressed as a product of 2 particle wavefunctions
defined each plaquette. For $\pi/2\leq\theta\leq 3\pi/2$, the ground state is
given by

$\begin{split}\ket{0}=\prod_{\bm{r}}\Big{[}&e^{iA_{xy}}a_{12}(\theta)c^{\dagger}_{1}(\bm{r})c^{\dagger}_{2}(\bm{r}+\hat{x}+\hat{y})\\\
&+a_{34}(\theta)c^{\dagger}_{3}(\bm{r}+\hat{x})c^{\dagger}_{4}(\bm{r}+\hat{y})\Big{]}\ket{vac},\\\
&a_{12}(\theta)=\sqrt{\frac{1-\sin(\theta)}{2}},\\\
&a_{34}(\theta)=\sqrt{\frac{1+\sin(\theta)}{2}}.\end{split}$ (19)

Third, for $3\pi/2\leq\theta\leq 2\pi$, the ring exchange terms vanish, and
the ground state can again be written as a product of single-particle
wavefunctions. Similar to Eq. 18 the ground state for $3\pi/2\leq\theta\leq
2\pi$ is given by

$\begin{split}\ket{0}=&\prod_{\bm{r}}e^{iA_{xy}}\Big{[}a_{1}(\theta)c^{\dagger}_{1}(\bm{r})+a_{2}(\theta)c^{\dagger}_{2}(\bm{r})+a_{3}(\theta)c^{\dagger}_{3}(\bm{r})\Big{]}\\\
&\times\Big{[}a_{1}(\theta)c^{\dagger}_{1}(\bm{r})-a_{2}(\theta)c^{\dagger}_{2}(\bm{r})+a_{4}(\theta)c^{\dagger}_{4}(\bm{r})\Big{]}\ket{vac},\\\
&a_{1}(\theta)=a_{2}(\theta)=\sqrt{\frac{1-\sin(\theta)}{4}},\\\
&a_{3}(\theta)=a_{4}(\theta)=-\sqrt{\frac{1+\sin(\theta)}{2}}.\end{split}$
(20)

Compared to Eq. 18, Eq. 20 differs by the addition of the phase $e^{iA_{xy}}$.
This phase is needed in order for the wavefunctions in Eq. 19 and Eq. 20 to
match at $\theta=3\pi/2$.

Figure 6: Schematic of the change in the charge distribution near the corner
of a lattice during the pumping process from Eq. 17. Here we have added a
charge $-2$ ion (in units of the electron charge) to each unit cell to make
the system charge neutral.

During this process, the quadrupole moment can be calculated using Eq. 13, and
the wavefunctions in Eq. 18-20. As a function of $\theta$, the quadrupole
moment is given by

$\Delta Q_{xy}=\begin{cases}0&\text{for }0\leq\theta\leq\pi/2\\\
\frac{1}{2}(1-\sin(\theta))&\text{for }\pi/2\leq\theta\leq 3\pi/2\\\
1&\text{for }3\pi/2\leq\theta\leq 2\pi\end{cases},$ (21)

where $\Delta Q_{xy}\equiv Q_{xy}(\theta)-Q_{xy}(0)$. These values of $\Delta
Q_{xy}$ are shown in Fig. 7. From Eq. 21 we can clearly see that after a full
period, the quadrupole increases by $1$. This corresponds to an integer amount
of charge being pumped to the corners, and agrees with our earlier claim that
the quadrupole moment of the system is only defined modulo $1$. Furthermore,
we can also confirm that $Q_{xy}(\pi)-Q_{xy}(0)=1/2$, and that when
$\theta=\pi$ the dipole conserving model is a HOTI protected by $C_{4}$
symmetry.

Since the model is exactly solvable over the full range of $\theta$ between
$0$ and $2\pi$ for the parameterization given in Eq. 17, we can also exactly
calculate the change in the corner charge of this system as $\theta$ is
varied. For a lattice of size $N_{x}\times N_{y}$ with open boundaries, the
ground state can be found by following the same steps we used in Eq. 18-20.
For this geometry, there is a net charge of $+\Delta Q_{xy}$ located at the
$(1,1)$ and $(N_{x},N_{y})$ sites, and a net charge of $-\Delta Q_{xy}$
located at the $(N_{x},1)$ and $(1,N_{y})$ sites, where $\Delta Q_{xy}$ is
defined as in Eq. 21 (the charge remains constant at all other sites). So, we
find that the quadrupole response does indeed predict the correct corner
charges. Additionally, we can confirm that the $C_{4}$ symmetric HOTI has
half-integer corner charges, as expected from the quadrupole response in Eq.
10. The change in the charge distribution as $\theta$ is increased from $0$ to
$2\pi$ is illustrated in Fig. 6.

Figure 7: The quadrupole moment $\Delta Q_{xy}\equiv Q_{xy}(\theta)-Q_{xy}(0)$
as a function of the adiabatic parameter $\theta$ from Eq. 21.

###  III.2 Continuous Dipole Pumping and Mean Field Theory

In the previous section, we established how an adiabatic process can lead to a
change in the quadrupole moment. However, despite being exactly solvable, the
pumping process in Eq. 17 is not continuous (it is only piecewise continuous).
Preparing for our eventual replacement of the adiabatic parameter $\theta$
with a new momentum $k_{z}$, we also want to analyze a parameterization of Eq.
15 that is fully continuous with respect to $\theta$, and displays the same
essential phenomenology we found using the piece-wise continuous
parameterization in Sec. III.1. Specifically, we want a periodic
parameterization where, as a function of $\theta,$ the model evolves from a
$C_{4}$ symmetric trivial insulator through a $C_{4}$ symmetric HOTI, and then
back to a trivial insulator. Based on Sec. III.1, we should find that the
quadrupole moment of the model should increase by $1/2$ during this evolution
from trivial to HOTI, and should change by $1$ after a full period.

With this in mind, we will consider the parameterization

$\begin{split}&t(\theta)=\gamma[1+\cos(\theta)],\\\
&V(\theta)=V(\text{const.}),\\\ &\mu(\theta)=\gamma\sin(\theta),\end{split}$
(22)

where $\gamma$ is a constant. The values of $t$, $V$ and $\mu$ as a function
of $\theta$ are plotted in Fig 8. As before, the model has $C_{4}$ symmetry
when $\mu=0$, which occurs when $\theta=0,\pi\mod(2\pi)$. When $\theta=\pi$,
the parameters in Eq. 17 and 22 are the same, and the ground state of the
model can be found exactly (see Eq. 19) Based on our results from Sec. III.1,
at $\theta=\pi$ the model is a HOTI with half-integer corner charges for all
values of $V$. Similarly, when $\theta=0$ and $V=0$, the parameters in Eq. 17
and 22 are the same and at this point the system is a trivial $C_{4}$
symmetric insulator. Since this system is gapped the system should remain in
this trivial phase up to some finite value of $V$. Beyond this value of $V$,
we expect that the will undergo a phase transition and become a $C_{4}$
symmetric HOTI.

Figure 8: The evolution of $t$, $V$, and $\mu$ in Eq. 15 as a function of
$\theta$ using the parameterization in Eq. 22 with $V=\gamma=1$.

For $V\neq 0$, this model cannot be exactly solved (except at $\theta=\pi$).
Because of this, we will use self-consistent mean field theory to analyze the
interacting model. The first step in this approximation is to decompose the
ring exchange interaction into terms that are quadratic in the lattice fermion
operators using a Hubbard-Stratonovich transformation. Here, we will choose to
use the following decomposition of the ring exchange term:

$\begin{split}&-Vc_{1}^{\dagger}(\bm{r})c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})c_{4}(\bm{r}+\hat{y})e^{iA_{xy}(\bm{r})}\\\
&\rightarrow\lambda_{1x}(\bm{r})c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{4}(\bm{r}+\hat{y})+\lambda_{2x}(\bm{r})c_{1}^{\dagger}(\bm{r})c_{3}(\bm{r}+\hat{x})\\\
&\phantom{\rightarrow}+\lambda_{1y}(\bm{r})c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})-\lambda_{2y}(\bm{r})c_{1}^{\dagger}(\bm{r})c_{4}(\bm{r}+\hat{y})\\\
&\phantom{\rightarrow}-\frac{2}{V}\lambda_{1x}(\bm{r})\lambda_{2x}(\bm{r})e^{-iA_{xy}(\bm{r})}-\frac{2}{V}\lambda_{1y}(\bm{r})\lambda_{2y}(\bm{r})e^{-iA_{xy}(\bm{r})}.\end{split}$
(23)

The equations of motion for the Hubbard-Stratonovich fields $\lambda$ are

$\begin{split}&\lambda_{1x}(\bm{r})=\frac{V}{2}e^{iA_{xy}(\bm{r})}c_{1}^{\dagger}(\bm{r})c_{3}(\bm{r}+\hat{x}),\\\
&\lambda_{2x}(\bm{r})=\frac{V}{2}e^{iA_{xy}(\bm{r})}c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{4}(\bm{r}+\hat{y}),\\\
&\lambda_{1y}(\bm{r})=-\frac{V}{2}e^{iA_{xy}(\bm{r})}c_{1}^{\dagger}(\bm{r})c_{4}(\bm{r}+\hat{y}),\\\
&\lambda_{1x}(\bm{r})=\frac{V}{2}e^{iA_{xy}(\bm{r})}c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x}).\end{split}$
(24)

As can be seen from the equations of motion, under a linear phase shift
$c(\bm{r})\rightarrow c(\bm{r})e^{i\bm{\beta}\cdot\bm{r}}$, the Hubbard-
Stratonovich fields transform as

$\begin{split}&\lambda_{1x}(\bm{r})\rightarrow\lambda_{1x}(\bm{r})e^{i\bm{\beta}\cdot\hat{x}},\phantom{==}\lambda_{2x}(\bm{r})\rightarrow\lambda_{2x}(\bm{r})e^{-i\bm{\beta}\cdot\hat{x}},\\\
&\lambda_{1y}(\bm{r})\rightarrow\lambda_{1y}(\bm{r})e^{i\bm{\beta}\cdot\hat{y}},\phantom{==}\lambda_{2y}(\bm{r})\rightarrow\lambda_{2y}(\bm{r})e^{-i\bm{\beta}\cdot\hat{y}}.\end{split}$
(25)

In the self-consistent mean field theory approximation, we assume that the
Hubbard-Stratonovich fields acquire an expectation value that satisfies the
equations of motion Eq. 24. Equivalently, this approximation can be
interpreted as increasing the number of flavors of lattice fermions that
couple to the Hubbard-Stratonovich fields from $1$ to $N$, and taking the
$N\rightarrow\infty$ limitZinn-Justin (2002). As we show in Appendix A, the
self-consistent values of $\lambda_{ai}$ ($a=1,2$, $i=x,y$) can be written as

$\begin{split}&\lambda_{1x}(\bm{r})=\lambda
e^{i\phi_{x}(\bm{r})+iA_{xy}(\bm{r})},\\\ &\lambda_{2x}(\bm{r})=\lambda
e^{-i\phi_{x}(\bm{r})},\\\ &\lambda_{1y}(\bm{r})=\lambda
e^{i\phi_{y}(\bm{r})+iA_{xy}(\bm{r})},\\\ &\lambda_{2y}(\bm{r})=\lambda
e^{-i\phi_{y}(\bm{r})},\end{split}$ (26)

where the phase fields $\phi_{x}$ and $\phi_{x}$ satisfy the relationships

$\begin{split}\Delta_{y}\phi_{x}(\bm{r})=\Delta_{x}\phi_{y}(\bm{r})=A_{xy}(\bm{r}).\end{split}$
(27)

The self-consistent values of $\lambda$ can be found numerically as a function
of $\theta$ and $V/\gamma,$ and are shown in Fig. 9. Due to the equations of
motion in Eq. 24, the phase fields $\phi_{x}$ and $\phi_{y}$ transform under a
gauge transformation $\Lambda$ as

$\begin{split}&\phi_{x}(\bm{r})\rightarrow\phi_{x}(\bm{r})+\Delta_{x}\Lambda(\bm{r}),\\\
&\phi_{y}(\bm{r})\rightarrow\phi_{y}(\bm{r})+\Delta_{y}\Lambda(\bm{r}),\end{split}$
(28)

and Eq. 27 is consistent with the rank-2 gauge symmetry.

Figure 9: The self-consistent values of $\lambda/\gamma$ as a function of
$\theta$ for $V/\gamma$ between $0.4$ and $4.0$.

If we fix the external gauge fields to vanish ($A_{0}=A_{xy}=0$), the
quadratic mean field Hamiltonian can be written in Fourier space as

$\begin{split}&H^{Q}_{\text{mf}}=\sum_{\vec{k}}\Big{[}\bm{c}^{\dagger}(\vec{k})h^{Q}_{\text{mf}}(\vec{k})\bm{c}(\vec{k})\Big{]},\\\
&h^{Q}_{\text{mf}}(\vec{k})=\mu\Gamma_{0}+t(\Gamma_{2}+\Gamma_{4})\\\
&\phantom{===}+\lambda\big{[}\cos(k_{x}+\phi_{x})\Gamma_{4}+\sin(k_{x}+\phi_{x})\Gamma_{3}\\\
&\phantom{===}+\cos(k_{y}+\phi_{y})\Gamma_{2}+\sin(k_{y}+\phi_{y})\Gamma_{1}],\\\
\end{split}$ (29)

where $\vec{k}=(k_{x},k_{y})$ is the $2$D lattice momentum. Here, we have left
the $\theta$-dependence of $t$, $\mu$ and $\lambda$ implicit. Up to a shift in
momentum due to $\phi_{x}$ and $\phi_{y}$, this mean field Hamiltonian is
equivalent to that of the non-interacting quadrupole insulator that was
analyzed in Ref. Benalcazar  _et al._ , 2017b. We note that this mean-field
Hamiltonian can be coupled to $A_{xy}$ fields using the values for the
$\lambda_{ai}$ in Eq. 26.

The spectrum of the mean field Hamiltonian is given by

$\epsilon(\bm{k})=\pm\sqrt{\epsilon_{x}^{2}(k_{x})+\epsilon_{y}^{2}(k_{y})+\mu^{2}},$
(30)

where $\epsilon_{i}(k_{i})=\sqrt{\lambda^{2}+t^{2}+2\lambda
t\cos(k_{i}+\phi_{i})},$ for $i=x,y$. The band structure of Eq. 29 consists of
two upper bands and two lower bands. There is a gap between these bands when
$\lambda/t\neq 1$ and/or $\mu\neq 0$. When $\lambda/t=1$ and $\mu=0,$ a gap
closing occurs at $(k_{x},k_{y})=(\pi-\phi_{x},\pi-\phi_{y})$. In Ref.
Benalcazar  _et al._ , 2017b it was shown that when $\mu=0$, this Hamiltonian
describes a trivial $C_{4}$ symmetric insulator for $\lambda/t<1,$ and a
$C_{4}$ symmetric HOTI with half-integer corner charges for $\lambda/t>1$.

With this in mind, we shall now consider the quadrupole response of the mean
field Hamiltonian, using our linear response formalism. To do this, we will
couple the mean field Hamiltonian to the gauge fields according to Eq. 26.
Here, we are primarily interested in finding the change in the quadrupole
moment as $\theta$ is increased from $0$ to $\pi$, as well as the total change
in the quadrupole moment as $\theta$ is increased from $0$ to $2\pi$. At these
values of $\theta$, the lattice model has $C_{4}$ symmetry, and the change in
quadrupole moment is quantized as an integer or half-integer. Indeed within
the mean field framework, we find that

$\begin{split}Q_{xy}(\theta=\pi)-Q_{xy}(\theta=0)=\begin{cases}\frac{1}{2}&\mbox{if
}V<V_{c}\\\ 0&\mbox{if }V>V_{c}\end{cases},\\\
Q_{xy}(\theta=2\pi)-Q_{xy}(\theta=0)=\begin{cases}1&\mbox{if }V<V_{c}\\\
0&\mbox{if }V>V_{c}\end{cases},\end{split}$ (31)

where $V_{c}\sim 2.5\gamma$. The details of this calculation are presented in
Appendix B.

For $V<V_{c}$ these results are consistent with an adiabatic process where a
trivial $C_{4}$ symmetric insulator adiabatically evolves into a $C_{4}$
symmetric HOTI as $\theta$ is increased from $0$ to $\pi$. Indeed, we see in
Fig. 9 that $\lambda(\theta=0)=0$ for $V<V_{C}$, and the mean field
Hamiltonian at $\theta=0$ is a trivial insulator. As noted before, at
$\theta=\pi$, the model can be solved exactly, and we can confirm that it is
indeed a $C_{4}$ symmetric HOTI. Additionally, we find that for this range of
$V$ the quadrupole moment of the model increases by $1$ after a full period.
This agrees with the results of Sec. III.1.

For $V>V_{c}$ the quadrupole moment of the model at $\theta=0$ is equal to the
quadrupole moment of the model at $\theta=\pi$. As noted before, at
$\theta=\pi$, the model is a HOTI for all values of $V$. This means that at
$\theta=0$ the model transitions from being a trivial insulator to being a
HOTI at $V=V_{c}$. In Fig. 9 we can see that at this transition the value of
$\lambda(\theta=0)$ jumps from $0$ to $\sim 2\gamma$. Indeed, for
$\lambda(\theta=0)>2\gamma$, the mean field Hamiltonian describes a non-
interacting HOTI with half-integer corner chargesBenalcazar  _et al._ (2017a).
In terms of the parameters of the lattice model in Eq. 15, this means that the
dipole conserving model is a $C_{4}$ symmetric trivial insulator for $\mu=0$,
$V/t\lesssim 1.25,$ and a $C_{4}$ symmetric HOTI for $\mu=0$, $V/t\gtrsim
1.25$. This phase diagram is shown in Fig. 10.

Figure 10: Top: The solutions to the self-consistent solutions for $\lambda/t$
as a function of $V/t$ at $\mu=0$. When $V/t\simeq 1.25$ (grey line), the
quadrupole moment of the model changes by $1/2$. Bottom: The mean field phase
diagram of Eq. 15 at $\mu=0$. Note that the relation between $t$ and $\gamma$
is given in Eq. 22.

As a final point, we would like to address the role of the phase fields
$\phi_{x}$ and $\phi_{y}$. Based on Eq. 28, the phase fields can be considered
as dynamic gauge fields. Since these fields are dynamic, they must be
integrated over. This integration will project out any states that transform
non-trivially under shifts in $\phi_{i}$. To see what states transform non-
trivially, we note that in Eq. 29 a shift in $\phi_{i}$ is equivalent to
shifting the lattice momentum $k_{i}$. Integrating over $\phi_{i}$ therefore
projects out any many-body states that transform non-trivial under a shift in
momentum. Momentum is dual to position, i.e., it is the position operator that
acts to shift the momentum. From our earlier discussion of polarization in
Sec. II, we can conclude that the many-body states that transform non-
trivially under a shift in momentum are the states with non-vanishing
polarization. So, the integration over the phase fields projects out any many-
body states with a non-vanishing polarization, and will ensure that the many-
body polarization is a good quantum number in the ground state of the mean
field Hamiltonian, as it is in the full interacting Hamiltonian. This is
crucial for applying the rank-2 formalism which relies on dipole conservation.
This can be shown explicitly, by considering the mean field Hamiltonian when
$\mu=t=0$ (this corresponds to $\theta=\pi$ in Eq. 22). For these parameters,
the mean field ground state is given by

$\begin{split}\ket{0}=&\prod_{\bm{r}}\Big{[}\frac{e^{-i\phi_{y}}c^{\dagger}_{1}(\bm{r})}{2}+\frac{e^{i\phi_{x}}c^{\dagger}_{2}(\bm{r}+\hat{x}+\hat{y})}{2}-\frac{c^{\dagger}_{4}(\bm{r}+\hat{y})}{\sqrt{2}}\Big{]}\\\
&\times\Big{[}\frac{e^{-i\phi_{x}}c^{\dagger}_{1}(\bm{r})}{2}-\frac{e^{i\phi_{y}}c^{\dagger}_{2}(\bm{r}+\hat{x}+\hat{y})}{2}-\frac{c^{\dagger}_{3}(\bm{r}+\hat{x})}{\sqrt{2}}]\ket{vac}.\end{split}$
(32)

After integrating over the phase field $\phi_{x}$ and $\phi_{y}$, Eq. 32
reduces to

$\ket{0}=\prod_{\bm{r}}\frac{1}{\sqrt{2}}\Big{[}c^{\dagger}_{1}(\bm{r})c^{\dagger}_{2}(\bm{r}-\hat{x}+\hat{y})+c^{\dagger}_{3}(\bm{r}+\hat{x})c^{\dagger}_{4}(\bm{r}+\hat{y})\Big{]}\ket{vac}.$
(33)

This is exactly the ground state of the full interacting model when $\mu=t=0$
(see Eq. 19). Here we can directly see that the integration over the phase
fields is necessary in order for the ground state of the non-dipole conserving
mean field Hamiltonian to match that of the dipole conserving interacting
Hamiltonian. We address this topic in more detail in Appendix C.

We note that this is a very interesting outcome, i.e., if we take the ground
state of the free fermion quadrupole model which has vanishing polarization in
the ground state because of symmetry quantization, we can couple it to phase
fields and after integrating them out we recover the ground state for the
ring-exchange quadrupole model which has microscopic dipole conservation at
the Hamiltonian level. Ref. Dubinkin and Hughes, 2020 showed these two ground
states were adiabatically connected in the presence of $C_{4}$ symmetry if the
dipole conservation was relaxed, and so our results forge a new connection
between these models. These results open the possibility to use a rank-2
quadrupole calculation to determine the quadrupole moment of free-fermion
systems if the polarization-free projection can be carried out as we did
above. We leave such a program to future work.

##  IV Dipolar Chern-Simons Response Action

It is well known that the polarization response of a $1$D system is related to
the Chern-Simons response of a $2$D system via dimensional reductionQi  _et
al._ (2008). Based on this, one would expect that the $2$D quadrupole response
in Eq. 10 is related to a topological response in $3$D. Such a $3$D
topological response was proposed in Ref. You  _et al._ , 2019. This response
describes a system with anomalous chiral hinge modes, and can be written in
terms of three background gauge fields, two rank-1 gauge fields $A_{0}$,
$A_{z}$, and 1 rank-2 gauge field $A_{xy}$,

$\begin{split}S_{\text{dCS}}=\frac{1}{4\pi}\int
d^{4}x[&A_{xy}\partial_{z}A_{0}+A_{0}\partial_{z}A_{xy}-A_{xy}\partial_{t}A_{z}\\\
&-A_{z}\partial_{t}A_{xy}+A_{z}\partial_{x}\partial_{y}A_{0}-A_{0}\partial_{x}\partial_{y}A_{z}].\end{split}$
(34)

This action breaks both $C_{4}$ rotation symmetry around the $z$-axis, and
time reversal symmetry $\mathcal{T}$, but is invariant under their product,
which we will refer to as $C_{4}\mathcal{T}$ symmetry. Eq. 34 is gauge
invariant up to surface terms, which shall be discussed later. Provided that
the gauge fields are non-singular, Eq. 34 can be written as a sum of total
derivatives. Because of this, the response action can be reduced to a sum of
boundary terms when defined on a manifold with boundary. The Lagrangians for
these boundary terms are

$\begin{split}&\mathcal{L}_{dCS,\pm
t}=\mp\frac{1}{4\pi}A_{z}A_{xy},\phantom{==}\mathcal{L}_{dCS,\pm
z}=\pm\frac{1}{4\pi}A_{0}A_{xy},\\\ &\mathcal{L}_{dCS,\pm
x}=\pm\frac{1}{4\pi}A_{z}\partial_{y}A_{0},\phantom{==}\mathcal{L}_{dCS,\pm
y}=\pm\frac{1}{4\pi}A_{z}\partial_{x}A_{0},\\\ &\mathcal{L}_{dCS,\pm x,\pm
y}=-\frac{1}{4\pi}A_{z}A_{0},\phantom{==}\mathcal{L}_{dCS,\pm x,\mp
y}=\frac{1}{4\pi}A_{z}A_{0}.\end{split}$ (35)

Here, $\mathcal{L}_{dCS,\pm\mu}$ ($\mu=x,y,z,t$) is the Lagrangian for the
boundaries oriented normal to the $\pm\mu$-direction. Similarly,
$\mathcal{L}_{dCS,\pm x,\pm y}$ is the Lagrangian for the hinges between
boundaries normal to the $\pm x$ and $\pm y$-direction, while
$\mathcal{L}_{dCS,\pm x,\mp y}$ is the Lagrangian for the hinges between
boundaries normal to the $\pm x$ and $\mp y$-direction.

Since the bulk action is a total derivative, all bulk currents vanish, and
there are only non-vanishing currents at the boundaries. These boundary
currents are anomalous, and their anomalous conservation laws are given by

$\begin{split}&\partial_{t}j^{\pm x}_{0}+\partial_{z}j^{\pm
x}_{z}=\pm\frac{1}{4\pi}[\partial_{y}\partial_{z}A_{0}-\partial_{y}\partial_{t}A_{z}]=\pm\frac{1}{4\pi}\partial_{y}E_{z},\\\
&\partial_{t}j^{\pm y}_{0}+\partial_{z}j^{\pm
y}_{z}=\pm\frac{1}{4\pi}[\partial_{x}\partial_{z}A_{0}-\partial_{x}\partial_{t}A_{z}]=\pm\frac{1}{4\pi}\partial_{x}E_{z},\\\
&\partial_{t}j^{\pm z}_{t}+\partial_{x}\partial_{y}j^{\pm
z}_{xy}=\mp\frac{1}{4\pi}[\partial_{x}\partial_{y}A_{0}-\partial_{t}A_{xy}]=\mp\frac{1}{4\pi}E_{xy},\\\
&\partial_{z}j^{\pm t}_{z}+\partial_{x}\partial_{y}j^{\pm
t}_{xy}=\pm\frac{1}{4\pi}[\partial_{x}\partial_{y}A_{z}-\partial_{z}A_{xy}]=\pm\frac{1}{4\pi}B,\\\
&\partial_{t}j^{\pm x\pm y}_{0}+\partial_{z}j^{\pm x\pm
y}_{z}=\pm\frac{1}{4\pi}[\partial_{z}A_{0}-\partial_{t}A_{z}]=\pm\frac{1}{4\pi}E_{z},\\\
&\partial_{t}j^{\pm x\mp y}_{0}+\partial_{z}j^{\pm x\mp
y}_{z}=\mp\frac{1}{4\pi}[\partial_{z}A_{0}-\partial_{t}A_{z}]=\mp\frac{1}{4\pi}E_{z},\end{split}$
(36)

where we have introduced the fields
$E_{xy}\equiv\partial_{x}\partial_{y}A_{0}-\partial_{t}A_{xy}$ and
$B\equiv\partial_{x}\partial_{y}A_{z}-\partial_{z}A_{xy}$, which can be
thought of as rank-2 electric and magnetic fields respectively. Here, the
superscript $\pm\mu$ ($\mu=x,y,z,t$) indicates that the current is defined on
boundaries oriented normal to the $\pm\mu$-direction. Similarly, the
superscript $\pm x\pm y$ indicates that the current is defined on hinges
between boundaries normal to the $\pm x$ and $\pm y$-direction, while $\pm
x\mp y$ indicates that the current is defined on hinges between boundaries
normal to the $\pm x$ and $\mp y$-direction. From Eq. 36, we see that a rank-2
electric field $E_{xy}$ produces an anomalous dipole current $j_{xy}$ on the
boundaries normal to the $\pm z$-direction, and a gradient of the electric
field $\partial_{i}E_{z}$ ($i=y,x$) produces an anomalous current in the
$z$-direction on boundaries normal to the $\pm x$ and $\pm y$-directions
respectively. There is also a chiral anomaly-like response on the hinges.

As noted before, Eq. 34 is only gauge invariant up to boundary terms. After a
gauge transformation $\Lambda$, the gauge variation at the boundaries is given
by

$\begin{split}&\delta\mathcal{L}_{dCS,\pm
t}=\mp\frac{1}{4\pi}\Lambda[\partial_{x}\partial_{y}A_{z}-\partial_{z}A_{xy}],\\\
&\delta\mathcal{L}_{dCS,\pm
z}=\pm\frac{1}{4\pi}\Lambda[\partial_{x}\partial_{y}A_{0}-\partial_{t}A_{xy}],\\\
&\delta\mathcal{L}_{dCS,\pm
x}=\pm\frac{1}{4\pi}\Lambda[\partial_{y}\partial_{t}A_{z}-\partial_{y}\partial_{z}A_{0}],\\\
&\delta\mathcal{L}_{dCS,\pm
y}=\pm\frac{1}{4\pi}\Lambda[\partial_{x}\partial_{t}A_{z}-\partial_{x}\partial_{z}A_{0}],\\\
&\delta\mathcal{L}_{dCS,\pm x,\pm
y}=-\frac{1}{4\pi}\Lambda[\partial_{t}A_{z}-\partial_{z}A_{0}],\\\
&\delta\mathcal{L}_{dCS,\pm x,\mp
y}=\frac{1}{4\pi}\Lambda[\partial_{t}A_{z}-\partial_{z}A_{0}].\end{split}$
(37)

In order for the full theory to be gauge invariant, there must be additional
degrees of freedom located at the boundaries and hinges in order to restore
gauge invariance. This is similar to what occurs when the $2$D Chern-Simons
action is defined on a manifold with boundary. Specifically, due to the hinge
terms $\delta\mathcal{L}_{dCD,\pm x\pm y}$, there must be chiral modes that
propagate along hinges between boundaries normal to the $\pm x$ and $\pm
y$-direction in order for the theory to be gauge invariant. Similarly, due to
$\delta\mathcal{L}_{dCD,\pm x\mp y}$ there must be anti-chiral modes that
propagate along hinges between boundaries normal to the $\pm x$ and $\mp
y$-direction. We can therefore conclude that a consistent theory described by
Eq. 34 must have chiral hinge modes. Additionally, due to the boundary terms,
$\mathcal{L}_{\pm x}$ and $\mathcal{L}_{\pm y}$, there must also be additional
modes at boundaries normal to the $\pm x$ and $\pm y$-directions in order for
the theory to be gauge invariant.

The dipolar Chern-Simons response also predicts a quadrupole analog of the
Laughlin pump. The typical quantum Hall Laughlin pumping process can be
observed by considering a $2$D Chern-Simons theory defined on cylinder. When a
unit of flux is inserted through this cylinder, the Hall current pumps charge
from one end of the cylinder to the other, and changes the dipole moment by
$1$. A similar process occurs when the dipolar Chern-Simons term is defined on
an annulus with periodic boundary conditions in the $z$-direction. In this
case, when a unit of flux is inserted in the $z$-direction, Eq. 34 predicts
that the $xy$ quadrupole moment of the system will increase by $1$. This
process can analogously be thought of as dipole pumping, and will change the
amount of charge located at the hinges of the system.

Finally, we would like to confirm that the dipolar Chern-Simons action is in
fact related to the $2$D quadrupole response by dimensional reduction. To do
this, we will consider an arbitrarily thin annulus with periodic boundary
conditions in the $z$-direction. In this limit, we can treat the flux $\theta$
passing through the $z$-direction as an adiabatic parameter of a $2$D theory,
and the dipole analog of the Laughlin pump becomes a dipole analog of a
Thouless pump. As $\theta$ is increased, the quadrupole moment of the $2$D
system is shifted due to an adiabatic pumping of dipole moment across the
system. This is exactly the phenomenology we saw before when considering the
quadrupole moment in Sec. III.

To show this formally, we can dimensionally reduce Eq. 34 by setting
$A_{z}=\Theta/L_{z}$, and taking the limit $L_{z}\rightarrow 0$. If we take
$\Theta$ to be constant, the dimensionally reduced action is given by

$\begin{split}S_{dCS,2D}=\int
d^{3}x\frac{\Theta}{4\pi}[\partial_{x}\partial_{y}A_{0}-\partial_{t}A_{xy}].\end{split}$
(38)

However, as discussed in Ref. You  _et al._ , 2019, the term in Eq. 38 is only
half the quadrupole response of the dimensionally reduced system. The other
half comes from the boundary degrees of freedom we discussed earlier. Taking
both of these contributions into account, the full response of the
dimensionally reduced system is given by

$\begin{split}S_{dCS,2D}=\int
d^{3}x\frac{\Theta}{2\pi}[\partial_{x}\partial_{y}A_{0}-\partial_{t}A_{xy}],\end{split}$
(39)

which is the quadrupole response from Eq. 10 with $Q_{xy}=\Theta/2\pi$.

##  V $C_{4}\mathcal{T}$ Symmetric Lattice model

Since the $2$D quadrupole response is related to the $3$D dipolar Chern-Simons
response via dimensional reduction, the $2$D dipole conserving insulator
presented in Sec. III should be similarly related to a chiral hinge insulator
in $3$D that realizes the dipolar Chern-Simons response. To go from the $2$D
model to the $3$D model, we will proceed in the usual fashion of identifying
the adiabatic parameter $\theta$ (which controls the dipole pumping of Eq. 15)
with the momentum along the $z$-direction, $k_{z}$. A local Hamiltonian should
depend smoothly on the momentum $k_{z}$, and so we will consider the
parameterization given in Eq. 22, which is fully continuous with respect to
$\theta$. If we substitute $\theta\rightarrow k_{z}$, the resulting $3$D
Hamiltonian is

$\begin{split}H^{dCS}&=\sum_{k_{z},\bm{r}}\bm{c}^{\dagger}(\bm{r})h^{z}(k_{z})\bm{c}(\bm{r})\\\
&\phantom{=}-Vc_{1}^{\dagger}(\bm{r})c^{\dagger}_{2}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})c_{4}(\bm{r}+\hat{y})\\\
&\phantom{=}+h.c.,\\\
h^{z}(k_{z})&=\gamma\sin(k_{z})\Gamma_{0}+\gamma[1+\cos(k_{z})](\Gamma_{2}+\Gamma_{4}).\end{split}$
(40)

An illustration of this Hamiltonian is shown in Fig. 11. Eq. 40 is invariant
under phase shifts that depend linearly on the $x$ and $y$ coordinates,
$\bm{c}(\bm{R})\rightarrow\bm{c}(\bm{R})e^{i(\beta_{1}x+\beta_{2}y)}$, and
conserves dipole in both the $x$ and $y$-directions. However, due to the
explicit dependence of $h^{z}$ on $k_{z}$, the Hamiltonian does not conserve
dipole in the $z$-direction. Eq. 40 also has $C_{4}\mathcal{T}$ symmetry and
$C_{2}$ spatial rotations in the $xy$ plane.

As expected from our discussion of the dipolar Chern-Simons response, Eq. 40
can be coupled to the background gauge fields $A_{0}$, $A_{z}$, and $A_{xy}$.
In real space, the minimally coupled Hamiltonian is given by

$\begin{split}H^{dCS}=\sum_{\bm{R}}&[\gamma\bm{c}^{\dagger}(\bm{R})T\bm{c}(\bm{R})+\gamma\bm{c}^{\dagger}(\bm{R}+\hat{z})T^{z}\bm{c}(\bm{R})e^{iA_{z}(\bm{R})}\\\
-&Vc_{1}^{\dagger}(\bm{R})c_{2}^{\dagger}(\bm{R}+\hat{x}+\hat{y})c_{3}(\bm{R}+\hat{x})c_{4}(\bm{R}+\hat{y})e^{iA_{xy}(\bm{R})}\\\
-&A_{0}(\bm{R})\bm{c}^{\dagger}(\bm{R})\bm{c}(\bm{R})+h.c.],\\\
T=&\frac{1}{2}[\Gamma_{2}+\Gamma_{4}],\phantom{==}T^{z}=\frac{1}{2}[i\Gamma_{0}+\Gamma_{2}+\Gamma_{4}],\end{split}$
(41)

where $\bm{R}\equiv(x,y,z)$ labels the points of a $3$D cubic lattice. The
matrix $T$ is made up of the intracell terms, and the matrix $T^{z}$ is made
up of intercell hopping terms in the $z$-direction. In this form, it is clear
that $T^{z}$ violates dipole conservation in the $z$-direction.

Figure 11: A cross section of the $3$D lattice model Eq. 41. The Hamiltonian
consists of intracell and $z$-direction hopping terms $h^{z}$ (blue) and a
ring exchange term $V$ in the $xy$ plane (orange).

###  V.1 Mean Field Analysis

To analyze the physics of the interacting $C_{4}\mathcal{T}$ symmetric lattice
model in Eq. 40, we will employ a self-consistent mean field theory approach.
This will largely mirror the mean field analysis of the $2$D dipole conserving
model in Sec. III.2. As we shall show, within the self-consistent mean field
framework the 4 particle ring exchange interaction $V$ causes a (Mott) gap to
form for the lattice fermions. Additionally, the resulting mean field
Hamiltonian exactly maps on to a known non-interacting HOTI with chiral hinge
modes.

Similar to Sec. III.2, we will decompose the quartic ring exchange terms via a
Hubbard-Stratonovich transformation,

$\begin{split}&-Vc_{1}^{\dagger}(\bm{R})c_{2}^{\dagger}(\bm{R}+\hat{x}+\hat{y})c_{3}(\bm{R}+\hat{x})c_{4}(\bm{R}+\hat{y})e^{iA_{xy}(\bm{R})}\\\
&\rightarrow\lambda_{1x}(\bm{R})c_{2}^{\dagger}(\bm{R}+\hat{x}+\hat{y})c_{4}(\bm{R}+\hat{y})\\\
&\phantom{\rightarrow}+\lambda_{2x}(\bm{R})c_{1}^{\dagger}(\bm{R})c_{3}(\bm{R}+\hat{x})\\\
&\phantom{\rightarrow}+\lambda_{1y}(\bm{R})c_{2}^{\dagger}(\bm{R}+\hat{x}+\hat{y})c_{3}(\bm{R}+\hat{x})\\\
&\phantom{\rightarrow}-\lambda_{2y}(\bm{R})c_{1}^{\dagger}(\bm{R})c_{4}(\bm{R}+\hat{y})\\\
&\phantom{\rightarrow}-\frac{2}{V}\lambda_{1x}(\bm{R})\lambda_{2x}(\bm{R})e^{-iA_{xy}(\bm{R})}\\\
&\phantom{\rightarrow}-\frac{2}{V}\lambda_{1y}(\bm{R})\lambda_{2y}(\bm{R})e^{-iA_{xy}(\bm{R})}.\end{split}$
(42)

The equations of motion are exactly those in Eq. 24 (upon exchange the $2$D
lattice coordinate $\bm{r}$ with the $3$D lattice coordinate $\bm{R}$). Under
a linear phase shift
$\bm{c}(\bm{R})\rightarrow\bm{c}(\bm{R})e^{i(\beta_{1}x+\beta_{2}y)}$, the
Hubbard-Stratonovich fields transform as

$\begin{split}&\lambda_{1x}(\bm{R})\rightarrow\lambda_{1x}(\bm{R})e^{i\beta_{1}},\phantom{==}\lambda_{2x}(\bm{R})\rightarrow\lambda_{2x}(\bm{R})e^{-i\beta_{1}},\\\
&\lambda_{1y}(\bm{R})\rightarrow\lambda_{1y}(\bm{R})e^{i\beta_{2}},\phantom{==}\lambda_{2y}(\bm{R})\rightarrow\lambda_{2y}(\bm{R})e^{-i\beta_{2}}.\end{split}$
(43)

As shown in Appendix A, in the self-consistent mean field theory approximation
the Hubbard-Stratonovich fields acquire expectation values of the form

$\begin{split}&\lambda_{1x}(\bm{R})=\lambda
e^{i\phi_{x}(\bm{R})+iA_{xy}(\bm{R})},\\\ &\lambda_{2x}(\bm{R})=\lambda
e^{-i\phi_{x}(\bm{R})},\\\ &\lambda_{1y}(\bm{R})=\lambda
e^{i\phi_{y}(\bm{R})+iA_{xy}(\bm{R})},\\\ &\lambda_{2y}(\bm{R})=\lambda
e^{-i\phi_{y}(\bm{R})},\end{split}$ (44)

where the phase fields $\phi_{x}$ and $\phi_{x}$ satisfy the relationship

$\begin{split}\Delta_{y}\phi_{x}(\bm{R})=\Delta_{x}\phi_{y}(\bm{R})=A_{xy}(\bm{R}).\end{split}$
(45)

The self-consistent values of $\lambda$ depend on $V/\gamma$, and are shown in
Fig. 12. As before, under a gauge transformation $\Lambda$, the phase fields
gauge transform as $\phi_{i}\rightarrow\phi_{i}+\Delta_{i}\Lambda$.

Figure 12: Top: The solutions to the self-consistent solutions for
$\lambda/\gamma$ as a function of $V/\gamma$. When $V/\gamma\simeq 2.1$ (grey
line), $\lambda/\gamma=2$ (dashed line) and the mean field model undergoes a
phase transition from a topological chiral hinge insulator to a weak HOTI.
Bottom: The mean field phase diagram of Eq. 41.

For vanishing gauge fields ($A_{0}=A_{z}=A_{xy}=0$) the mean field Hamiltonian
can be written in Fourier space as

$\begin{split}&H^{dCS}_{\text{mf}}=\sum_{\bm{k}}\Big{[}\bm{c}^{\dagger}(\bm{k})h^{dCS}_{\text{mf}}(\bm{k})\bm{c}(\bm{k})\Big{]},\\\
&h^{dCS}_{\text{mf}}(\bm{k})=\gamma\sin(k_{z})\Gamma_{0}+\gamma[1+\cos(k_{z})](\Gamma_{2}+\Gamma_{4})\\\
&\phantom{====}+\lambda\big{[}\cos(k_{x}+\phi_{x})\Gamma_{4}+\sin(k_{x}+\phi_{x})\Gamma_{3}\\\
&\phantom{====}+\cos(k_{y}+\phi_{y})\Gamma_{2}+\sin(k_{y}+\phi_{y})\Gamma_{1}\big{]},\end{split}$
(46)

where $\bm{k}=(k_{x},k_{y},k_{z})$ is the momentum of the $3$D lattice model.
As a consistency check, we note that Eq. 46 is related to the $2$D mean field
model in Eq. 29 by dimensional reduction. The single-particle energy spectrum
of the mean field Hamiltonian is given by

$\begin{split}\epsilon(\bm{k})=&\pm\big{[}7\gamma^{2}+4\lambda^{2}+8\gamma^{2}\cos(k_{z})+\gamma^{2}\cos(2k_{z})\phantom{\big{)}}\\\
&\phantom{\big{(}}+4\lambda\gamma\cos(k_{x}+\phi_{x})(\cos(k_{z})+1)\phantom{\big{)}}\\\
&\phantom{\big{(}}+4\lambda\gamma\cos(k_{y}+\phi_{y})(\cos(k_{z})+1)\big{]}^{1/2}.\end{split}$
(47)

When $\lambda/\gamma=2$, the system is gapless at
$(k_{x},k_{y},k_{z})=(\pi-\phi_{x},\pi-\phi_{y},0)$, and when $\lambda=0$, the
system is gapless at $k_{z}=\pi$ for all values of $k_{x}$ and $k_{y}$. As
noted in Sec. III.2, integrating over the dynamic phase fields $\phi_{x}$ and
$\phi_{y}$ projects out any states with non-vanishing dipole moment in the $x$
and $y$-directions.

In Sec. III.2, we found that in the mean field limit the $2$D dipole
conserving model is equivalent to the non-interacting quadrupole insulator.
Here, we see that the mean field Hamiltonian in Eq. 46 is equivalent to a non-
interacting chiral hinge insulator, which was also considered in Ref.
Benalcazar  _et al._ , 2017b, and is related to the non-interacting quadrupole
insulator by dimensional reduction. Based on this, we shall quote several key
results concerning the mean field model. When $0<\lambda/\gamma<2$ the mean
field Hamiltonian describes a topological chiral hinge insulator. In this
phase, the model has modes that propagate in the $z$-direction on hinges
between boundaries normal to the $x$ and $y$-directions. In the corresponding
non-interacting HOTI, these chiral hinge modes correspond to the non-trivial
Wannier-band Chern number of the system. For $2<\lambda/\gamma$, there are no
chiral hinge modes. This phase is connected to the $V\rightarrow\infty$ limit
of Eq. 40. In this limit, the ring exchange term gaps out each $xy$-layer of
the model individually. This suppresses all tunneling in the $z$-direction and
the $3D$ system becomes a stack of $2$D insulators. Based on our analysis in
Sec. III.1, each of these layers will have quadrupole moment $1/2$, and the
system can be adiabatically connected to an atomic insulator by breaking
translational symmetry along the $z$-direction, so we might refer to it as a
weak higher order topological insulator. At the phase transition connecting
these two phases ($\lambda/\gamma=2$), there is a bulk band crossing and the
system is gapless. Using the numerical solutions to the self-consistent mean
field equations, we find that this phase transition occurs at $V\equiv
V_{c}\simeq 2.1\gamma$. The mean field phase diagram is shown in Fig. 12.

Having determined the mean field structure of the lattice model, we can now
turn our attention to determining the effective response action for this
system. As noted in Sec. IV, the dipolar Chern-Simons action can be written as
a sum of total derivatives, and the only non-trivial responses occur at the
boundaries. Because of this, we will need to define our model on a lattice
with boundary in order to show that it exhibits a dipolar Chern Simons
response. As we will show, analyzing the boundaries of our model is made
tractable upon passing to the continuum in the $z$ and $t$-directions. This
analysis will be the main topic of the next subsection.

###  V.2 Continuum Analysis

In this subsection, we will study the $3$D model given in Eq. 41 in the
continuum, near $V=0$. At $V=0$, the system is quasi-1d, i.e., a decoupled
$3$D array of $1$D wires oriented in the $z$-direction. If we pass to the
continuum along the $z$-direction, each of these wires consists of 4 massless
Dirac fermions (which correspond to the fluctuations of the lattice fermions
with momentum near $k_{z}=\pi$), and 4 massive fermions (which correspond to
the fluctuations of the lattice fermions with momentum near $k_{z}=0$). We can
split lattice fermions into these two contributions using the identification

$c_{i}(\bm{R})\rightarrow\Psi_{i}(\bm{R})+(-1)^{z}\psi_{i}(\bm{R}),$ (48)

where $\Psi_{i}$ are the heavy fermions and $\psi_{i}$ are the light fermions.
At $V=0$, the continuum Lagrangian for these fields is

$\begin{split}&\mathcal{L}^{dCS}=\sum_{\bm{r}}\left[\bm{\psi}^{\dagger}(\bm{r})G^{-1}_{0}(\omega+A_{0}(\bm{r}),p_{z}+A_{z}(\bm{r}))\bm{\psi}(\bm{r})\right.\\\
&\phantom{====}+\left.\bm{\Psi}^{\dagger}(\bm{r})G^{-1}_{M}(\omega+A_{0}(\bm{r}),p_{z}+A_{z}(\bm{r}))\bm{\Psi}(\bm{r})\right],\\\
&G^{-1}_{0}(\omega,p_{z})=\omega\text{I}+p_{z}\Gamma^{0},\\\
&G^{-1}_{M}(\omega,p_{z})=\omega\text{I}-p_{z}\Gamma^{0}-\tfrac{M}{\sqrt{2}}(\Gamma^{2}+\Gamma^{4}),\end{split}$
(49)

where I is the $4\times 4$ identity matrix,
$\bm{\psi}=(\psi_{1},\psi_{2},\psi_{3},\psi_{4}),$ and
$\bm{\Psi}=(\Psi_{1},\Psi_{2},\Psi_{3},\Psi_{4})$. Here, we have only passed
to the continuum in the $z$-direction. The $x$ and $y$ coordinates are still
lattice coordinates. Based on Eq. 49, when $V=0$, $\psi_{1}$ and $\psi_{2}$
are anti-chiral fermions and $\psi_{3}$ and $\psi_{4}$ are chiral fermions
(with respect to their propagation along the $z$-direction). Additionally, the
heavy fermions $\Psi$ are massive with mass $M$. We will also take $M$ to be
the UV cutoff for this theory.

We can now consider the interactions in this theory. As before, we will use
the Hubbard-Stratonovich transformation in Eq. 42 to decompose the lattice
ring exchange interaction. Using Eq. 48 we find that the Hubbard-Stratonovich
fields $\lambda_{1/2,x}$ $\lambda_{1/2,y}$ couple to both the light fermions
$\psi$ and the heavy fermions $\Psi$. Since the heavy fermions $\Psi$ are
gapped, they can be integrated out, leaving a Lagrangian in terms of the light
fermions and the Hubbard-Stratonovich fields. After integrating out the heavy
fermions, the Lagrangian can be written as

$\mathcal{L}^{dCS}=\mathcal{L}_{\psi}+\mathcal{L}_{\lambda\lambda},$ (50)

where $\mathcal{L}_{\psi}$ contains all terms involving the light fermions
$\psi,$ and $\mathcal{L}_{\lambda\lambda}$ contains all couplings between the
Hubbard-Stratonovich fields. If we ignore any terms that osculate like
$(-1)^{z}$, $\mathcal{L}_{\psi}$ is given by

$\begin{split}\mathcal{L}_{\psi}=&\sum_{\bm{r}}\left[\bm{\psi}^{\dagger}(\bm{r})G^{-1}_{0}(\omega+A_{0}(\bm{r}),p_{z}+A_{z}(\bm{r}))\bm{\psi}(\bm{r})\right.\\\
&-[\lambda_{2x}(\bm{r})\psi_{1}^{\dagger}(\bm{r})\psi_{3}(\bm{r}+\hat{x})\\\
&\phantom{=}+\lambda_{1x}(\bm{r})\psi_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})\psi_{4}(\bm{r}+\hat{y})\\\
&\phantom{=}+\lambda_{2y}(\bm{r})\psi_{1}^{\dagger}(\bm{r})\psi_{4}(\bm{r}+\hat{y})\\\
&\phantom{=}+\left.\lambda_{1y}(\bm{r})\psi_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})\psi_{3}(\bm{r}+\hat{x})+h.c.]\right].\end{split}$
(51)

To determine $\mathcal{L}_{\lambda\lambda}$ we must integrate out the heavy
fermions $\Psi$. At one loop order, $\mathcal{L}_{\lambda\lambda}$ is given by

$\begin{split}\mathcal{L}_{\lambda\lambda}=&\sum_{\bm{r}}-u_{1}[\lambda_{1x}(\bm{r})\lambda^{*}_{1x}(\bm{r})+\lambda_{2x}(\bm{r})\lambda^{*}_{2x}(\bm{r})\\\
&+\lambda_{1y}(\bm{r})\lambda^{*}_{1y}(\bm{r})+\lambda_{2y}(\bm{r})\lambda^{*}_{2y}(\bm{r})]\\\
&+[u_{2}\lambda_{1x}(\bm{r})\lambda_{2x}(\bm{r}+\hat{y})+u_{2}\lambda_{1y}(\bm{r})\lambda_{2y}(\bm{r}+\hat{x})\\\
&\phantom{=}+\frac{2}{V}\lambda_{1x}(\bm{r})\lambda_{2x}(\bm{r})e^{-iA_{xy}(\bm{r})}\\\
&\phantom{=}+\frac{2}{V}\lambda_{1y}(\bm{r})\lambda_{2y}(\bm{r})e^{-iA_{xy}(\bm{r})}+h.c.],\end{split}$
(52)

where $u_{1}=\frac{\log(4)-1}{16\pi},$ and $u_{2}=\frac{1}{16\pi}$.

As before, we shall employ the self-consistent mean field theory
approximation. We find the Hubbard-Stratonovich fields acquire expectation
values of the form

$\begin{split}\lambda_{1x}(\bm{r})&=\lambda
e^{i\phi_{x}(\bm{r})+iA_{xy}(\bm{r})},\\\ \lambda_{2x}(\bm{r})&=\lambda
e^{-i\phi_{x}(\bm{r})},\\\ \lambda_{1y}(\bm{r})&=\lambda
e^{i\phi_{y}(\bm{r})+iA_{xy}(\bm{r})},\\\ \lambda_{2y}(\bm{r})&=\lambda
e^{-i\phi_{y}(\bm{r})}.\\\ \end{split}$ (53)

The value of $\lambda$ is determined by the effective potential

$\begin{split}H_{\lambda\lambda}=&\frac{8}{V^{\prime}}\lambda^{2}+\frac{1}{2\pi}\lambda^{2}\Big{(}\log(\frac{2\lambda^{2}}{M^{2}})-1\Big{)}.\end{split}$
(54)

where $V^{\prime}=(V^{-1}+\frac{u_{1}+u_{2}}{2})^{-1}$. The effective
potential is minimized by
$\lambda=\frac{1}{\sqrt{2}}Me^{-\frac{8\pi}{V^{\prime}}}$. In agreement with
the numerical results, we find that $\lambda$ vanishes when $V\rightarrow 0,$
and it increases monotonically with increasing $V$. Additionally, due to the
$u_{2}$ term in Eq. 52, at low energies the phase fields $\phi_{i}$ satisfy

$\begin{split}\Delta_{y}\phi_{x}(\bm{r})=\Delta_{x}\phi_{y}(\bm{r})=A_{xy}(\bm{r}).\\\
\end{split}$ (55)

These results agree with the lattice mean field results from Sec. V.1. Here,
as before, under a gauge transformation $\Lambda$, the phase fields $\phi_{i}$
transform as
$\phi_{i}(\bm{r})\rightarrow\phi_{i}(\bm{r})+\Delta_{i}\Lambda(\bm{r})$.

Based on Eq. 51 we see that in the mean field limit, the quartet fermion
cluster $\psi_{1}(\bm{r})$, $\psi_{2}(\bm{r}+\hat{x}+\hat{y})$,
$\psi_{3}(\bm{r}+\hat{x})$, and $\psi_{4}(\bm{r}+\hat{y})$ couples to one
another for each value of in plane lattice coordinate ${\bm{r}}.$ Hence, the
continuum Lagrangian decouples into 4-fermion clusters, which are defined on
each plaquette of the $xy$-plane. This feature will make it possible for us to
analytically consider boundaries normal to the $x$ and $y$-directions, which
we do in the next section. For $\lambda\neq 0$, these fermions all become
massive (see Eq. 51). From this, we can explicitly confirm that this model has
chiral hinge modes. As we show in Fig 13, at the top right hinge there is a
net chiral mode. There are also similar hinge modes at the other hinges of the
system. In addition, there is both a chiral mode and an anti-chiral mode
located at each at each lattice site along the boundaries normal to the $x$
and $y$-direction (see Fig 13). Since these modes all come in pairs they can
be gapped out with local symmetry-persevering perturbations.

Figure 13: The continuum model Eq. 49. In each unit cell there are 2 chiral
($\bullet$) and 2 anti-chiral ($\times$) modes. The ring exchange interaction
(orange), gaps out these 4 modes (grey circles) at each plaquette. The modes
along the edges remain gapless (white circles), but can be gapped when
symmetry-preserving perturbations are added. An odd number of chiral hinge
modes is stabilized on each corner of the $xy$ plane.

###  V.3 Effective Response Theory

We will now turn our attention to finding the bulk and boundary responses for
the $3$D $C_{4}\mathcal{T}$ symmetric chiral hinge insulator. The response
action for this system is composed of terms that depend on the three
background gauge fields, $A_{0}$, $A_{z}$ and $A_{xy},$ as well as the two
phase fields $\phi_{x}$ and $\phi_{y}$. The phase fields must be included
since they depend on the background gauge field $A_{xy}$ (see Eq. 55). The
coefficients for the various terms in the response action are determined by
the current-current correlation functionsFradkin (2013). Here the current-
current correlation functions are determined using the continuum mean field
Lagrangian from Sec. V.2.

First, we shall consider the effective response action for the bulk of the
system. In the limit of low frequency and momentum, the effective Lagrangian
for the bulk is

$\begin{split}\mathcal{L}^{dCS}_{eff,\text{bulk}}=&(2A_{xy}-\partial_{y}\phi_{x}-\partial_{x}\phi_{y})\partial_{t}A_{z}\\\
-&(2A_{xy}-\partial_{y}\phi_{x}-\partial_{x}\phi_{y})\partial_{z}A_{0}\\\
=&0,\end{split}$ (56)

where we have passed to the continuum in the $x$ and $y$-directions, and used
Eq. 55 in the last line. This result is consistent with the dipolar Chern-
Simons response action given in Eq. 34, since that response action is a total
derivative.

To probe the non-vanishing the boundary effects, we can consider the mean
field Lagrangian defined on a lattice with boundary. As noted in Sec. V.2,
there are gapless fermions at the boundaries, but they can be gapped out with
symmetry preserving perturbations. Upon doing so, we can find that the
effective response Lagrangian for the background gauge fields and the phase
fields at low frequency and momentum is given by

$\begin{split}L^{dCS}_{eff,\pm
y}=&\pm\frac{1}{4\pi}[A_{0}\partial_{x}A_{z}+\phi_{x}\partial_{z}A_{0}-\phi_{x}\partial_{t}A_{z}],\\\
L^{dCS}_{eff,\pm
x}=&\pm\frac{1}{4\pi}[A_{0}\partial_{y}A_{z}+\phi_{y}\partial_{z}A_{0}-\phi_{y}\partial_{t}A_{z}].\end{split}$
(57)

The first term on the right hand sides of each line of Eq. 57 is exactly the
boundary term of the dipolar Chern-Simons response action. The remaining two
terms are the couplings between the gauge fields $A_{0}$ and $A_{z}$ and the
phase fields $\phi_{x}$ and $\phi_{y}$. With the addition of the phase fields,
the boundary Lagrangian is gauge invariant up to hinge terms.

Finally, we can consider the action on the hinges. As noted before, there is a
gapless chiral mode at the $\pm x,\pm y$ hinges and a gapless anti-chiral mode
at the $\pm x,\mp y$ hinges. If we ignore these gapless fermions, the
effective hinge Lagrangian is given by

$\begin{split}L^{dCS}_{eff,\pm x,\pm y}=&-\frac{1}{4\pi}A_{0}A_{z},\\\
L^{dCS}_{eff,\pm x,\mp y}=&+\frac{1}{4\pi}A_{0}A_{z}.\\\ \end{split}$ (58)

This is exactly the hinge term from the dipolar Chern-Simons response. From
Eq. 56-58 we can confirm that the effective action of Eq. 40 matches the
dipolar Chern-Simons action in the bulk as well as boundaries normal to the
$x$ and $y$-directions, and the hinges that separate them.

It will also be useful to consider how the effective Lagrangian transforms
under a gauge transformation. Using Eq. 57 and 58, we find that under a gauge
transformation $\Lambda$, the effective Lagrangian is shifted by the hinge
term

$\begin{split}&\delta\mathcal{L}^{dCS}_{eff,\pm x,\pm
y}=-\frac{1}{4\pi}\Lambda[\partial_{t}A_{z}-\partial_{z}A_{0}],\\\
&\delta\mathcal{L}^{dCS}_{eff,\pm x,\mp
y}=\frac{1}{4\pi}\Lambda[\partial_{t}A_{z}-\partial_{z}A_{0}],\end{split}$
(59)

which is the same hinge term given in Eq. 37. Gauge invariance is restored by
the aforementioned chiral hinge modes of the system (which were ignored in our
derivation of the effective Lagrangian). It is well known that a single chiral
fermion is not gauge invariant, and the gauge variation of a single chiral
fermion exactly cancels out the gauge variation $\delta L_{eff,\pm x,\pm
y}$Naculich (1988); Chandrasekharan (1994). Similarly, the gauge variation of
a single anti-chiral fermion exactly cancels out the gauge variation of
$L_{eff,\pm x,\mp y}$. From this, we can conclude that the lattice degrees of
freedom (the phase fields and the chiral hinge modes) cancel out the gauge
anomalies of the dipolar Chern-Simons action (Eq. 37) and make the full theory
gauge invariant, as desired.

As a final point, we note that in this analysis we have relied on
translational invariance in the $t$ and $z$-directions. Because of this we
cannot make any statements about the boundaries normal the $t$ or
$z$-directions. Determining how to analyze these boundaries is an interesting
topic for further research. In particular, due to the coupling to the rank-2
gauge field $A_{xy}$, boundaries normal to the $z$-direction may host
interesting and exotic physics.

##  VI $M\mathcal{T}$ Symmetric Lattice model

In Sec. V we considered a dipole conserving $C_{4}\mathcal{T}$ symmetric $3$D
model with chiral hinge modes. In this section, we will consider a related
dipole conserving $3$D model, which breaks $C_{4}\mathcal{T}$ symmetry, but is
instead invariant under $M_{x}\mathcal{T}$ and $M_{y}\mathcal{T}$. For
simplicity we shall use the shorthand $M\mathcal{T}$ to refer to both of the
symmetries. As we shall show, the $M\mathcal{T}$ model also displays chiral
hinge modes, similar to those of the $C_{4}\mathcal{T}$ symmetric model, but
has other interesting features that distinguish it from the previous case.

Here, we shall consider the Hamiltonian

$\begin{split}H^{M\mathcal{T}}&=\sum_{k_{z},\bm{r}}\left[\bm{c}^{\dagger}(\bm{r})h^{M\mathcal{T}}(k_{z})\bm{c}(\bm{r})\right.\\\
&\phantom{=}-Vc_{1}^{\dagger}(\bm{r})c^{\dagger}_{2}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})c_{4}(\bm{r}+\hat{y})\\\
&\phantom{=}+\left.h.c.\right],\\\
h^{M\mathcal{T}}(k_{z})&=\gamma\sin(k_{z})\Gamma_{0}+[1+\cos(k_{z})](\gamma^{\prime}\Gamma_{2}+\gamma\Gamma_{4}).\end{split}$
(60)

As desired, this model is invariant under phase shifts that depend linearly on
the $x$ and $y$ coordinates,
$\bm{c}(\bm{R})\rightarrow\bm{c}(\bm{R})e^{i(\beta_{1}x+\beta_{2}y)}$. For
$\gamma^{\prime}=\gamma,$ Eq. 62 has $C_{4}\mathcal{T}$ and is the same as Eq.
40. When $\gamma^{\prime}\neq\gamma$, the $C_{4}\mathcal{T}$ symmetry is
explicitly broken, and Eq. 60 only has $M\mathcal{T}$ symmetry. The
$M_{x}\mathcal{T}$ and $M_{y}\mathcal{T}$ symmetries act on the fermionic
degrees of freedom in each unit cell as $\bm{c}\rightarrow U_{M_{x}}\bm{c}$
and $\bm{c}\rightarrow U_{M_{y}}\bm{c}$ respectively, where
$U_{M_{x}}=\tau_{1}\otimes\tau_{3},$ and $U_{M_{y}}=\tau_{1}\otimes\tau_{1}$.
Here we are interested in the case where $C_{4}\mathcal{T}$ symmetry is
explicitly broken, and without loss of generality we will take
$\gamma^{\prime}<\gamma$.

###  VI.1 Mean Field Analysis

As before, we shall use self-consistent mean field theory to analyze the
interacting $M\mathcal{T}$ symmetric lattice model. To do this, we shall
decompose the ring exchange interaction using the decomposition from Eq. 42.
In the self-consistent mean field limit, the Hubbard-Stratonovich fields
acquire the expectation values of the form

$\begin{split}&\lambda_{1x}(\bm{R})=\lambda_{x}e^{i\phi_{x}(\bm{R})+iA_{xy}(\bm{R})},\\\
&\lambda_{2x}(\bm{R})=\lambda_{x}e^{-i\phi_{x}(\bm{R})},\\\
&\lambda_{1y}(\bm{R})=\lambda_{2y}(\bm{R})=0,\\\
&\Delta_{y}\phi_{x}(\bm{R})=A_{xy}(\bm{R}).\end{split}$ (61)

The self-consistent values of $\lambda_{x}$ depend on both $V/\gamma$ and
$\gamma^{\prime}/\gamma,$ and can be found numerically. The self-consistent
values of $\lambda_{x}/\gamma$ as a function of $V/\gamma$ for
$\gamma^{\prime}/\gamma=1/2$ are shown in Fig. 14 (details of this calculation
are given in Appendix A). Unlike the self-consistent solutions for the
$C_{4}\mathcal{T}$ model in Eq. 44, these solutions break $C_{4}\mathcal{T}$
symmetry, but remain invariant under $M\mathcal{T}$ symmetry. Since
$\lambda_{1y}(\bm{R})=\lambda_{2y}(\bm{R})=0$ in the mean field limit, there
are no terms that hop fermions along the $y$-direction, although there are
terms that hop fermions along the $x$-direction. Because of this, in the mean-
field limit, this model will be quasi-$2$D, i.e., a decoupled array of $2$D
planes stacked along the $y$-direction.Here, the planes are stacked along the
$y$-direction because we chose to set $\gamma^{\prime}<\gamma$ in Eq. 60. If
instead we had chosen to set $\gamma^{\prime}>\gamma$, the planes would be
stacked along the $x$-direction.

Figure 14: The self-consistent solutions of $\lambda_{x}/\gamma$ as a function
of $V/\gamma$ for $\gamma^{\prime}/\gamma=1/2$. At the mean field level, the
bulk of the system remains gapped for all values of $\lambda_{x}$.

For periodic geometries, the quasi-2D nature of the mean field Hamiltonian can
be made manifest by writing it in the following form:

$\begin{split}H^{M\mathcal{T}}_{\text{mf}}=&\sum_{\bm{k}^{\prime},y}\bm{c}^{\dagger}(\bm{k}^{\prime},y)h^{M\mathcal{T}}_{\text{mf}}(\bm{k}^{\prime},y)\bm{c}(\bm{k}^{\prime},y),\\\
h^{M\mathcal{T}}_{\text{mf}}(\bm{k}^{\prime},y)=&\gamma\sin(k_{z})\Gamma_{0}+[1+\cos(k_{z})](\gamma^{\prime}\Gamma_{2}+\gamma\Gamma_{4})\big{)}\\\
&+\lambda_{x}\big{(}\cos(k_{x}+\phi_{x})\Gamma_{4}+\sin(k_{x}+\phi_{x})\Gamma_{3}\big{)},\\\
\end{split}$ (62)

where $\bm{k}^{\prime}=(k_{x},k_{z})$, and we have omitted the gauge fields
for simplicity. Here, we can directly see that the mean field Hamiltonian
describes a layered system where each layer is an $xz$-plane, and the
$y$-coordinate labels the different layers. Because the different layers are
fully decoupled, the mean field Hamiltonian explicitly conserves dipole in the
$y$-direction. As before, the phase field $\phi_{x}$ projects out states with
non-vanishing polarization in the $x$-direction.

When $\gamma^{\prime}=0$, each layer in Eq. 62 is composed of two decoupled
$2$D insulators, one of which depends only on the $c_{1}$ and $c_{3}$
fermions, and one of which depends only on the $c_{2}$ and $c_{4}$ fermions.
When $0<\lambda_{x}/\gamma<2,$ the first of these insulators has Chern number
$-1,$ and the second has Chern number $+1$. Since the two insulators have
opposite Chern number, there are no net chiral modes associated with such a
layer, and any surface modes can be gapped out with symmetry preserving
perturbations. The $\gamma^{\prime}$ term couples these two insulators and
turns them into a single insulator with vanishing Chern-number. When
$\lambda_{x}/\gamma>2$ both the insulators have Chern number 0.

The energy spectrum of the mean field Hamiltonian is given by

$\begin{split}\epsilon(\bm{k}^{\prime})=&\Big{(}\frac{3}{2}\gamma^{\prime
2}+2\gamma^{2}+\lambda_{x}^{2}+4[\gamma^{\prime 2}+\gamma^{2}]\cos(k_{z})\\\
&+2\gamma\lambda_{x}[1+\cos(k_{z})]\cos(k_{x}+\phi_{x})+\gamma^{\prime}\cos(2k_{z})\Big{)}^{1/2}.\end{split}$
(63)

For $\gamma^{\prime}\neq 0$, the spectrum is gapped for all values of
$\lambda_{x}/\gamma>0$. When $\gamma^{\prime}=0$, and $\lambda_{x}/\gamma=2$
there is a gap closing at $(k_{x},k_{z})=(\pi-\phi_{x},0)$. As noted before,
at this point the two decoupled insulators in each layer transition from
having Chern number $\pm 1$ to having Chern number 0.

To determine the existence of any hinge modes in this model, we will consider
a lattice with open boundaries along the $y$-direction. Such a geometry is
straightforward to analyze given the quasi-$2$D nature of the mean field
Hamiltonian. Specifically, let us consider a lattice of length $N_{y}$ in the
$y$-direction. In the mean field limit, this corresponds to $N_{y}$ decoupled
$2$D systems, which are indexed by their $y$-coordinate. For the layers away
from the boundaries ($y\neq 1,N_{y}$), the mean field Hamiltonian is the same
as $h(\bm{k}^{\prime},y)$ given in Eq. 62. Based on Eq. 61, in the $y=1$ layer
(i.e., the boundary normal to the $-\hat{y}$-direction), there are no Hubbard-
Stratonovich fields coupling the $c_{2}$ and $c_{4}$ lattice fermions, and the
mean field Hamiltonian is given by

$\begin{split}h^{M\mathcal{T}}_{\text{mf}}(\bm{k}^{\prime},1)=&\gamma\big{(}\sin(k_{z})\Gamma_{0}+[1+\cos(k_{z})](\gamma^{\prime}\Gamma_{2}+\gamma\Gamma_{4})\big{)}\\\
&+\lambda\big{(}\cos(k_{x}+\phi_{x})\Gamma^{+}_{4}+\sin(k_{x}+\phi_{x})\Gamma^{+}_{3}\big{)},\\\
\end{split}$ (64)

where $\Gamma^{\pm}_{i}=(\Gamma_{i}\pm\eta\Gamma_{i}\eta)/2$, and
$\eta=\text{diag}(-1,1,1,1)$. For $\lambda_{x}<\lambda^{c}_{x}\equiv
2(\gamma+2\gamma^{\prime 2}/\gamma)$, this layer has Chern number $-1$. At
$\lambda_{x}=\lambda^{c}_{x}$ this boundary layer is gapless, and for
$\lambda_{x}>\lambda^{c}_{x}$ the boundary layer is an insulator with Chern
number $0$. At the $y=L_{y}$ layer (i.e., the boundary normal to the
$+\hat{y}$-direction), there are no Hubbard-Stratonovich fields coupling the
$c_{1}$ and $c_{3}$ lattice fermions, and the mean field Hamiltonian is given
by

$\begin{split}h^{M\mathcal{T}}_{\text{mf}}(\bm{k}^{\prime},N_{y})=&\gamma\big{(}\sin(k_{z})\Gamma_{0}+[1+\cos(k_{z})](\gamma^{\prime}\Gamma_{2}+\gamma\Gamma_{4})\big{)}\\\
&+\lambda\big{(}\cos(k_{x}+\phi_{x})\Gamma^{-}_{4}+\sin(k_{x}+\phi_{x})\Gamma^{-}_{3}\big{)}.\\\
\end{split}$ (65)

Similar to before, for $0<\lambda_{x}<\lambda_{x}^{c}$, this layer has Chern
number $+1$. At $\lambda_{x}=\lambda_{x}^{c}$ this boundary layer is also
gapless, and for $\lambda_{x}>\lambda_{x}^{c}$ this boundary layer is an
insulator with Chern number $0$. It is worth noting that the boundary
Hamiltonians in Eq. 64 and 65 can be modified by the addition of additional of
symmetry preserving terms that are localized at the boundaries of the system.
Because of this, the value of value of $\lambda^{c}_{x}$ is not uniquely
determined by the bulk of the system.

As we have seen, the mean field model can be interpreted as a layered system,
where for $0<\lambda_{x}<\lambda_{x}^{c}$ each layer in the bulk consists of
two $2$D insulators with opposite Chern numbers, which are coupled via
$\gamma$. The layers at boundaries normal to the $\pm\hat{y}$-direction, have
Chern number $\pm 1$. This is shown schematically in Fig. 15. Based on this,
we can confirm the existence of chiral hinge modes in the $M\mathcal{T}$
symmetric model for $0<\lambda_{x}<\lambda_{x}^{c}$. In terms of the lattice
parameters, this means that the chiral hinge modes persist up to a critical
value of $V$ that depends on the value of $\gamma^{\prime}/\gamma$, as well as
any additional boundary terms. At $\lambda_{x}=\lambda_{x}^{c}$ the mean field
model transitions from a chiral hinge insulator to a weak HOTI. Provided that
$\gamma^{\prime}\neq 0$, the energy-gap of the bulk of the system remains open
during this transition at the mean field level. Only the single-particle
energy gap of the boundary closes. In contrast, we found that when the
$C_{4}\mathcal{T}$ model transitions from being a chiral hinge insulator to a
weak HOTI, the bulk energy-gap closes. Phase transitions where only the
boundary energy-gap closes have been previously studied in the context of
boundary obstructed topological phasesBenalcazar  _et al._ (2017a); Khalaf
_et al._ (2019). For these systems, it has been shown that while the bulk
energy gap does not close during such a phase transition, the Wannier gap in
the bulk does close, and signals a transition between the topologically
distinct states.

While the mean-field Hamiltonian we have derived is only a boundary obstructed
topological phase, let us comment on the initial interacting model. It is
worth remarking that in non-dipole conserving models, e.g., free-fermion band
theories, chiral hinge modes cannot be protected by $M\mathcal{T}$ symmetry
alone in $3$D. This can be concluded from the fact that if we consider a $3$D
model with chiral hinge modes, it is possible to add $2$D insulators with
Chern number $\pm 1$, to the boundaries normal to the $\mp\hat{y}$ direction
without breaking $M\mathcal{T}$ symmetry. This will cause there to be both a
chiral and anti-chiral fermion at each hinge, which can be gapped out with
symmetry preserving perturbations. However, adding a $2$D insulator with Chern
number $\pm 1$ to the boundaries of a system necessarily violates dipole
conservation. This can be concluded from the fact that adiabatically shifting
the momentum of an insulator with non-vanishing Chern number causes the
insulator to polarize, via the Laughlin pumping process. This clearly violates
dipole conservation (lattice models of Chern insulators also have single-
particle tunneling terms and tunable orbital magnetization, both of which
violate dipole conservation). Because of this, the chiral hinge modes in the
$M\mathcal{T}$ model we have considered here are protected, and cannot be
gapped out without breaking dipole conservation in either the $x$ or
$y$-direction.

Figure 15: Schematic of the interacting $M\mathcal{T}$ symmetric model in Eq.
60 with open boundary conditions (top), and the corresponding mean field
Hamiltonian (bottom). The mean field Hamiltonian describes a system composed
of $2$D layers stacked along the $y$-direction. In the bulk, these layers have
Chern number $0$, while the layers at the boundaries normal to the $\pm
y$-direction have Chern number $\pm 1$.

Having established the existence of the chiral hinge modes, would now like to
use the mean field Hamiltonian to find the response action for the chiral
hinge insulator phase of the $M\mathcal{T}$ model. This can be done
efficiently due to the quasi-2D nature of the mean field model. If we
reintroduce the gauge fields to the mean field Hamiltonian, we find that in
the bulk (Eq. 61) the response action vanishes. At the boundaries normal to
the $x$ and $y$-directions, we find that in the topological phase, the
boundary and hinge responses are given by

$\begin{split}L_{eff,\pm x}=0,\\\ L_{eff,\pm
y}=\pm&\frac{1}{2\pi}\big{[}A_{0}\partial_{x}A_{z}+\phi_{x}\partial_{z}A_{0}-\phi_{x}\partial_{t}A_{y}\big{]},\\\
L_{eff,\pm x,\pm y}=-&\frac{1}{4\pi}A_{0}A_{z}.\\\ L_{eff,\pm x,\mp
y}=+&\frac{1}{4\pi}A_{0}A_{z}.\end{split}$ (66)

In Appendix D, we show that this response can also be derived from a continuum
analysis, similar to what was done in Sec. V.2. Here, we can see that the
effective response vanishes for boundaries normal to the $x$-direction, while
for boundaries normal to the $y$-direction, the response action is equivalent
to a $2$D Chern-Simons action, if we identify the phase field $\phi_{x}$ with
the $x$-direction component of a rank-1 gauge field, $A_{x}$ (recall that the
phase field gauge transforms as
$\phi_{x}\rightarrow\phi_{x}+\partial_{x}\Lambda$). This $2$D Chern-Simons
term can be understood from the fact that in the mean field limit, the layers
normal to the $\pm\hat{y}$-direction have Chern number $\pm 1$. As expected,
the anomalies at the hinges are canceled by the inclusion of chiral fermionic
modes.

Comparing the mean field analysis of the $M\mathcal{T}$ model to that of the
$C_{4}\mathcal{T}$ model in Sec. V, we find that both of these models have a
topological phase with chiral hinge modes that persists up to a finite value
of ring exchange amplitude $V$, after which they transition into a weak HOTI
phase with no chiral hinge modes (though it does have non-chiral hinge modes
if we preserve translation symmetry along the hinge). Despite this, we find
that the effective surface Lagrangian of the $M\mathcal{T}$ model (Eq. 66)
differs significantly from that of the $C_{4}\mathcal{T}$ model (Eq. 57). For
the $C_{4}\mathcal{T}$ model, the boundary Lagrangian is non-vanishing for all
four boundaries, and all boundary terms have a prefactor of $\frac{1}{4\pi}$.
However, for the $M\mathcal{T}$ model, the boundary Lagrangian is non-
vanishing for the boundaries normal only to the $\pm\hat{y}$-directions, and
these terms have a prefactor of $\frac{1}{2\pi}$. In terms of the bulk
response action, this difference corresponds to the addition of a boundary
term

$\delta
S_{eff}=\int\frac{d^{4}x}{4\pi}[\partial_{y}(A_{0}\partial_{x}A_{z})-\partial_{x}(A_{0}\partial_{y}A_{z})].$
(67)

This term is invariant under $M\mathcal{T}$, but not $C_{4}\mathcal{T}$, as
expected. Importantly, this term does not affect the hinge terms, which are
the same for the $M\mathcal{T}$ model and the $C_{4}\mathcal{T}$ model. So, it
appears that the $C_{4}\mathcal{T}$ and $M\mathcal{T}$ lattice models share
the same hinge physics, but have different boundary physics.

###  VI.2 Dimensional Reduction

We will now consider dimensionally reducing the $M\mathcal{T}$ model in Eq. 60
from $3$D to $2$D. This should result in a $2$D model parameterized for dipole
pumping, analogous to the model given in Sec. III. If we use the standard
procedure of identifying the $z$-direction momentum $k_{z}$ with an adiabatic
parameter $\theta$, the dimensionally reduced Hamiltonian is given by

$\begin{split}H^{QM}&=\sum_{\bm{r}}\left[\bm{c}^{\dagger}(\bm{r})h^{oM}\bm{c}(\bm{r})\right.\\\
&\phantom{=}-Vc_{1}^{\dagger}(\bm{r})c^{\dagger}_{2}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})c_{4}(\bm{r}+\hat{y})\\\
&\phantom{=}+\left.h.c.\right],\\\
h^{oM}&=\gamma\sin(\theta)\Gamma_{0}+[1+\cos(\theta)](\gamma^{\prime}\Gamma_{2}+\gamma\Gamma_{4}).\end{split}$
(68)

For $\gamma^{\prime}=\gamma$, this is exactly the model from Eq. 15 with the
parameterization given in Eq. 22. For $\gamma^{\prime}\neq\gamma$ the model
has $M_{x}$ and $M_{y}$ symmetry when $\theta=0$, but still has $C_{4}$
symmetry when $\theta=\pi$. As before we will take $\gamma^{\prime}<\gamma$.

Following our analysis in Sec. III.2, we will use mean field theory to
determine the behavior of Eq. 68 at various values of $\theta$. Using the same
decomposition of the ring exchange term given in Eq. 23, we find that in the
mean field limit, the Hubbard-Stratonovich fields acquire expectation values
of the form

$\begin{split}&\lambda_{1x}(\bm{r})=\lambda_{x}e^{i\phi_{x}(\bm{r})+iA_{xy}(\bm{r})},\\\
&\lambda_{2x}(\bm{r})=\lambda_{x}e^{-i\phi_{x}(\bm{r})},\\\
&\lambda_{1y}(\bm{r})=\lambda_{2y}(\bm{r})=0,\\\
&\Delta_{y}\phi_{x}(\bm{r})=A_{xy}(\bm{r}).\end{split}$ (69)

Here, the self-consistent values of $\lambda_{x}$ depend on $\theta$,
$V/\gamma$, and $\gamma^{\prime}/\gamma$. The self-consistent values of
$\lambda_{x}$ are shown in Fig. 16.

Figure 16: The self-consistent values of $\lambda_{x}/\gamma$ as a function of
$\theta$ for $\gamma^{\prime}/\gamma=1/2$ and values of $V/\gamma$ between
$0.4$ and $4.0.$

As can be seen from Eq. 69, in the mean field limit there are no terms that
hop a single fermion along the $y$-direction. Because of this, the mean field
fermionic Hamiltonian is quasi-$1$D. Using the self-consistent values of
$\lambda_{x},$ the mean field Hamiltonian on periodic geometries can be
written in a manifestly quasi-$1$D form as

$\begin{split}H^{QM}_{\text{mf}}=&\sum_{k_{x},y}\bm{c}^{\dagger}(k_{x},y)h^{QM}_{\text{mf}}(k_{x},y)\bm{c}(k_{x},y),\\\
h^{QM}_{\text{mf}}(k_{x},y)=&\gamma\big{(}\sin(\theta)\Gamma_{0}+[1+\cos(\theta)](\gamma^{\prime}\Gamma_{2}+\gamma\Gamma_{4})\big{)}\\\
&+\lambda_{x}\big{(}\cos(k_{x}+\phi_{x})\Gamma_{4}+\sin(k_{x}+\phi_{x})\Gamma_{3}\big{)}.\\\
\end{split}$ (70)

The energy spectrum of the mean field Hamiltonian is the same as in Eq. 63,
upon identifying $k_{z}\rightarrow\theta$ (and treating $\lambda_{x}$ a
function of $\theta$). From this, we see that for $\gamma^{\prime}\neq 0$, the
mean field Hamiltonian remains gapped for all values of $\lambda_{x}/\gamma$.

To simplify the description, we note that the mean field Hamiltonian in Eq. 70
can be treated as an array of wires which are aligned along the $x$-direction,
and indexed by their $y$ coordinate. When $\gamma^{\prime}=0$, each wire in
Eq. 70 is composed of two decoupled $1$D insulators, one which only depends on
the $c_{1}$ and $c_{3}$ fermion operators, and one which only depends on the
$c_{2}$ and $c_{4}$ fermion operators. The first of these insulators is
equivalent to the SSH Hamiltonian given in Eq. 6, with parameterization

$(u,v,\mu)=(\gamma[1+\cos(\theta)],\lambda_{x},\gamma\sin(\theta)),$ (71)

while the second is equivalent to the SSH Hamiltonian with parameterization

$(u,v,\mu)=(\gamma[1+\cos(\theta)],\lambda_{x},-\gamma\sin(\theta)).$ (72)

Using Eq. 5 we find that these two insulators have opposite polarization
responses. In particular, as $\theta$ is varied, the two insulators pump an
opposite amount of charge, such that the net boundary charge remains
unchanged. This is to be expected from dimensional reduction since, in the
mean field limit of the $3$D model (Eq. 62), each layer consists of two
decoupled insulators with opposite Chern-numbers when $\gamma^{\prime}=0$. The
$\gamma^{\prime}$ term in Eq. 70 couples the SSH models at each layer,
resulting in a single insulator with vanishing polarization for all values of
$\theta$.

To determine the existence of any corner modes, we will consider a system with
open boundaries in the $y$-direction of length $N_{y}$. For the layers away
from the boundaries ($y\neq 1,N_{y}$), the mean field Hamiltonian is the same
as $h(k_{x},y)$ given in Eq. 70. At the $y=1$ layer (i.e., the boundary normal
to the $-\hat{y}$-direction), there are no Hubbard-Stratonovich fields
coupling the $c_{2}$ and $c_{4}$ lattice fermions, and the mean field
Hamiltonian is equivalent to the boundary Hamiltonian from Eq. 64 upon
identifying $k_{z}\rightarrow\theta$. Similarly, the mean field Hamiltonian at
the $y=N_{y}$ layer (i.e., the boundary normal to the $+\hat{y}$-direction) is
equivalent to the boundary Hamiltonian from Eq. 65 upon identifying
$k_{z}\rightarrow\theta$.

In the quasi-$1$D limit this model takes, the corner charges are contributed
only by the polarization of the boundary layers. The boundary polarization can
be calculated using Eq. 5. We find that, as a function of $\theta$, the change
in polarization of the $y=1$ layer is equal and opposite to the change in
polarization of the $y=L_{y}$ layer. Provided that $V$ exceeds a finite
critical value that depends on $\gamma^{\prime}/\gamma$, the polarization of
the $y=L_{y}$ ($y=1$ layer) increases (decreases) as $\theta$ is increased.
Over a full period, we find that the polarization of the boundary normal to
the $\pm\hat{y}$-direction changes by $\pm 1$. Additionally, when the mean
field model has $M_{x}$ and $M_{y}$ symmetry at $\theta=\pi$, the boundary
polarization is $\pm 1/2$. This value of the boundary polarization indicates
that there will be half-integer charges localized at the corners of the
system. As discussed in Sec. III, for systems coupled to the rank-2 gauge
fields, half-integer corner charges correspond to a quadrupole moment
$Q_{xy}=1/2$ (relative to a trivial insulator). To confirm that the quadrupole
moment at $\theta=\pi$ is in-fact $1/2$, we can couple the mean field
Hamiltonian to the rank-2 gauge field and use linear response to determine the
quadrupole moment (see Appendix B for further details). Upon doing this, we
indeed find that $Q_{xy}(\theta=\pi)-Q_{xy}(\theta=0)=1/2$ for $V$ less than
the aforementioned critical value.

As a final point, we would like to compare the mean field analysis of the two
$2$D dipole conserving Hamiltonians we have considered so far. In Sec. III.2
we considered a mean field decomposition of Eq. 15 that is manifestly
invariant under $C_{4}$ symmetry, and showed that in the mean field limit, the
model remains fully $2$D (see Eq. 29). In this section, we considered a mean
field decomposition of Eq. 70 that is invariant under only $M_{x}$ and $M_{y}$
symmetry, and showed that in the mean field limit, the model is quasi-$1$D
(equivalent to an array of decoupled $1$D wires). However, if we compare the
two interacting models, we see that at $\theta=\pi$ Eq. 15 and 70 are the
same, and both Hamiltonians consist only of the ring exchange term $V$.
Because of this, at $\theta=\pi$ the two mean field models should be
equivalent. While this is not immediately obvious, it is true that the two
mean field models are equivalent at $\theta=\pi$ provided we integrate over
the phase fields $\phi_{x}$ and $\phi_{y}$. As discussed in Sec. III.2, when
using the $C_{4}$ symmetric mean field decomposition, the ground state of the
mean field Hamiltonian at $\theta=\pi$ is given by Eq. 32, and upon
integrating over the phase fields, it is equivalent to the ground state of the
full interacting model Eq. 33. Similarly, when using the $M_{x}$ and $M_{y}$
symmetric mean field decomposition, the ground state at $\theta=\pi$ can be
written in real space as

$\begin{split}\ket{0}=&\prod_{\bm{r}}\Big{[}\frac{e^{i\phi_{x}}c^{\dagger}_{2}(\bm{r}+\hat{x}+\hat{y})}{\sqrt{2}}-\frac{c^{\dagger}_{4}(\bm{r}+\hat{y})}{\sqrt{2}}\Big{]}\\\
&\times\Big{[}\frac{e^{-i\phi_{x}}c^{\dagger}_{1}(\bm{r})}{\sqrt{2}}-\frac{c^{\dagger}_{3}(\bm{r}+\hat{x})}{\sqrt{2}}\Big{]}\ket{vac}.\end{split}$
(73)

Upon integrating over the phase fields $\phi_{x}$ we see that this ground
state is also equal to the ground state of the full interacting model in Eq.
33. We can therefore conclude that although the mean field Hamiltonians appear
very different at $\theta=\pi$, they do in fact lead to the same physics,
provided that the dynamic phase fields are properly accounted for.

##  VII Conclusion and Outlook

In this work we considered how various topological multipolar theories can be
realized in dipole conserving lattice models. In $2$D we considered the
quadrupole response and used a linear response formalism to determine the
quadrupole moment of a dipole conserving many-body system. We showed that
during a periodic adiabatic process, the quadrupole moment can change by an
integer, and that there is a $C_{4}$ symmetric HOTI with quadrupole moment
$1/2$ relative to a trivial insulator. We verified these results by
considering an exactly solvable interacting lattice model. We also showed that
for a properly chosen parameterization, this $2$D model can be related to a
$3$D $C_{4}\mathcal{T}$ symmetric model with chiral hinge modes. The boundary
responses of this $3$D model match those of the dipolar Chern-Simons action,
which is related to the rank-2 quadrupole response via dimensional reduction.
We also considered a variation of the $3$D model, that is instead invariant
under $M_{x}\mathcal{T}$ and $M_{y}\mathcal{T}$ symmetries, and showed that
they had similar hinge modes and topological responses. It is worth restating
that some of the analyses of these models has relied on self-consistent mean
field theory (which is equivalent to taking increasing the number of flavors
of fermions from $1$ to $N$ and taking the limit $N\rightarrow\infty$). While
we believe that the results of this approximation are valid, it would be
useful to find a way to verify these results that do not rely on the mean
field approximation.

An interesting facet of our analysis is that the interacting models we
considered here map exactly onto known non-interacting HOTIs in the mean field
limit. The $2$D dipole conserving $C_{4}$ symmetric HOTI and the $3$D dipoles
conserving $C_{4}\mathcal{T}$ symmetric HOTI we presented here, correspond to
the non-interacting $2$D quadrupole insulator, and $3$D chiral hinge insulator
of Ref. Benalcazar  _et al._ , 2017b respectively. The dipole conserving
models and non-interacting models also display the same topological phenomena
(quantized corner charges in $2$D and chiral hinge modes in $3$D). This
connection is interesting, since in Ref. Benalcazar  _et al._ , 2017b, the
topological nature of the non-interacting models is manifest in the non-
trivial topology of the Wannier bands. In the interacting models presented
here, there are no Wannier bands (outside the mean field limit) and the
topological nature of the models is manifest in non-trivial rank-2 topological
responses. Based on this, it appears that there may be deeper connections
between rank-2 gauge theories and HOTIs than initially expected.

It should be possible to simulate the models we have presented here in cold
atom systems. Due to the exceptional degree of control in cold atom systems,
it is possible to tune a system such that any single-particle hopping terms
are negligible, and the dominant terms are interactions, such as the ring-
exchange interaction we considered hereDai  _et al._ (2016). This should allow
for an explicit construction of the $2$D model we have considered here, and a
simulation of the pumping process we discussed. In principle it should also be
possible to construct a cold atom analog $C_{4}\mathcal{T}$ chiral hinge
insulator as well. Beyond cold atoms, other metamaterial systems exhibiting
tunable non-linear couplings may be able to realize the models we have
considered here. In known solid state materials there are no rank-2 gauge
fields that couple to the charge degrees of freedom. However, provided we are
treating the rank-2 gauge fields as background probes, it is possible to
identify the rank-2 gauge fields with the derivatives of physical rank-1 gauge
fields in certain situationsDubinkin  _et al._ (2019). It would be an
interesting question for further research to consider how the rank-1 analogs
of the phenomena we discussed here could occur in physical materials, and how
they could be measured in experiments.

## Acknowledgments

We thank O. Dubinkin for helpful discussions. JMM is supported by the National
Science Foundation Graduate Research Fellowship Program under Grant No. DGE –
1746047. TLH thanks the US Office of Naval Research (ONR) Multidisciplinary
University Research Initiative (MURI) grant N00014-20- 1-2325 on Robust
Photonic Materials with High-Order Topological Protection for support.

## References

  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Reviews of modern physics 82, 3045 (2010).
  * Senthil (2015) T. Senthil, Annu. Rev. Condens. Matter Phys. 6, 299 (2015).
  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Physical review letters 95, 226801 (2005).
  * Bernevig  _et al._ (2006) B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, science 314, 1757 (2006).
  * Kitaev (2009) A. Kitaev, in  _AIP conference proceedings_ , Vol. 1134 (American Institute of Physics, 2009) pp. 22–30.
  * Ryu  _et al._ (2010) S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. Ludwig, New Journal of Physics 12, 065010 (2010).
  * Hasan and Moore (2011) M. Z. Hasan and J. E. Moore, Annu. Rev. Condens. Matter Phys. 2, 55 (2011).
  * Chen  _et al._ (2013) X. Chen, Z.-C. Gu, Z.-X. Liu, and X.-G. Wen, Physical Review B 87, 155114 (2013).
  * Kapustin  _et al._ (2015) A. Kapustin, R. Thorngren, A. Turzillo, and Z. Wang, Journal of High Energy Physics 2015, 1 (2015).
  * Su  _et al._ (1979) W. Su, J. Schrieffer, and A. J. Heeger, Physical review letters 42, 1698 (1979).
  * Jackiw and Schrieffer (1981) R. Jackiw and J. R. Schrieffer, Nuclear Physics B 190, 253 (1981).
  * Fu  _et al._ (2007) L. Fu, C. L. Kane, and E. J. Mele, Physical review letters 98, 106803 (2007).
  * Moore and Balents (2007) J. E. Moore and L. Balents, Physical Review B 75, 121306 (2007).
  * Klitzing  _et al._ (1980) K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
  * King-Smith and Vanderbilt (1993a) R. King-Smith and D. Vanderbilt, Physical Review B 47, 1651 (1993a).
  * Qi  _et al._ (2008) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Physical Review B 78, 195424 (2008).
  * Qi  _et al._ (2013) X.-L. Qi, E. Witten, and S.-C. Zhang, Physical Review B 87, 134519 (2013).
  * Zhang  _et al._ (1989) S. C. Zhang, T. H. Hansson, and S. Kivelson, Physical review letters 62, 82 (1989).
  * Laughlin (1981) R. B. Laughlin, Physical Review B 23, 5632 (1981).
  * Thouless (1983) D. Thouless, Phys. Rev. B 27, 6083 (1983).
  * Fu (2011) L. Fu, Physical Review Letters 106, 106802 (2011).
  * Hughes  _et al._ (2011) T. L. Hughes, E. Prodan, and B. A. Bernevig, Physical Review B 83, 245132 (2011).
  * Hsieh  _et al._ (2012) T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, and L. Fu, Nature communications 3, 1 (2012).
  * Isobe and Fu (2015) H. Isobe and L. Fu, Physical Review B 92, 081304 (2015).
  * Cheng  _et al._ (2016) M. Cheng, M. Zaletel, M. Barkeshli, A. Vishwanath, and P. Bonderson, Physical Review X 6, 041068 (2016).
  * Bradlyn  _et al._ (2017) B. Bradlyn, L. Elcoro, J. Cano, M. Vergniory, Z. Wang, C. Felser, M. Aroyo, and B. A. Bernevig, Nature 547, 298 (2017).
  * Benalcazar  _et al._ (2017a) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017a).
  * Benalcazar  _et al._ (2017b) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Physical Review B 96, 245115 (2017b).
  * Langbehn  _et al._ (2017) J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Physical review letters 119, 246401 (2017).
  * Song  _et al._ (2017a) Z. Song, Z. Fang, and C. Fang, Physical review letters 119, 246402 (2017a).
  * Serra-Garcia  _et al._ (2018) M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Nature 555, 342 (2018).
  * Peterson  _et al._ (2018) C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. Bahl, Nature 555, 346 (2018).
  * Imhof  _et al._ (2018) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, _et al._ , Nature Physics 14, 925 (2018).
  * Khalaf  _et al._ (2019) E. Khalaf, W. A. Benalcazar, T. L. Hughes, and R. Queiroz, arXiv preprint arXiv:1908.00011 (2019).
  * Song  _et al._ (2017b) H. Song, S.-J. Huang, L. Fu, and M. Hermele, Physical Review X 7, 011020 (2017b).
  * Song and Schnyder (2017) X.-Y. Song and A. P. Schnyder, Physical Review B 95, 195108 (2017).
  * Rasmussen and Lu (2018) A. Rasmussen and Y.-M. Lu, arXiv preprint arXiv:1810.12317 (2018).
  * Thorngren and Else (2018) R. Thorngren and D. V. Else, Physical Review X 8, 011040 (2018).
  * You  _et al._ (2019) Y. You, F. Burnell, and T. L. Hughes, arXiv preprint arXiv:1909.05868 (2019).
  * Nandkishore and Hermele (2019a) R. M. Nandkishore and M. Hermele, Ann. Rev. of Cond. Matt. Phys. 10, 295 (2019a).
  * Pretko  _et al._ (2020) M. Pretko, X. Chen, and Y. You, Int. J of Mod. Phys. A 35, 2030003 (2020).
  * Pretko (2017a) M. Pretko, Physical Review B 95, 115139 (2017a).
  * Pretko (2017b) M. Pretko, Physical Review B 96, 035119 (2017b).
  * Pretko (2018) M. Pretko, Physical Review B 98, 115134 (2018).
  * Nandkishore and Hermele (2019b) R. M. Nandkishore and M. Hermele, Annual Review of Condensed Matter Physics 10, 295 (2019b).
  * Haah (2011) J. Haah, Physical Review A 83, 042330 (2011).
  * Vijay  _et al._ (2016) S. Vijay, J. Haah, and L. Fu, Physical Review B 94, 235157 (2016).
  * Hsieh and Halász (2017) T. H. Hsieh and G. B. Halász, Physical Review B 96, 165105 (2017).
  * Pretko (2017c) M. Pretko, Physical Review D 96, 024051 (2017c).
  * Ma  _et al._ (2018) H. Ma, M. Hermele, and X. Chen, Physical Review B 98, 035111 (2018).
  * Gromov (2019) A. Gromov, Physical review letters 122, 076403 (2019).
  * Xu (2006) C. Xu, arXiv preprint cond-mat/0602443 (2006).
  * Xu and Hořava (2010) C. Xu and P. Hořava, Physical Review D 81, 104033 (2010).
  * Pretko and Radzihovsky (2018) M. Pretko and L. Radzihovsky, Physical review letters 120, 195301 (2018).
  * Seiberg (2019) N. Seiberg, arXiv preprint arXiv:1909.10544 (2019).
  * Seiberg and Shao (2020) N. Seiberg and S.-H. Shao, arXiv preprint arXiv:2004.00015 (2020).
  * Dubinkin  _et al._ (2020) O. Dubinkin, A. Rasmussen, and T. L. Hughes, Annals of Physics , 168297 (2020).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Reviews of Modern Physics 83, 1057 (2011).
  * Resta and Vanderbilt (2007) R. Resta and D. Vanderbilt, in  _Physics of Ferroelectrics_ (Springer, 2007) pp. 31–68.
  * Resta (1992) R. Resta, Ferroelectrics 136, 51 (1992).
  * King-Smith and Vanderbilt (1993b) R. King-Smith and D. Vanderbilt, Physical Review B 47, 1651 (1993b).
  * Resta (1994) R. Resta, Reviews of modern physics 66, 899 (1994).
  * Dubinkin  _et al._ (2019) O. Dubinkin, J. May-Mann, and T. L. Hughes, arXiv preprint arXiv:1909.07403 (2019).
  * Zinn-Justin (2002) J. Zinn-Justin,  _Quantum field theory and critical phenomena_ , Vol. 113 (Clarendon Press, Oxford, 2002).
  * Dubinkin and Hughes (2020) O. Dubinkin and T. L. Hughes, arXiv preprint arXiv:2002.08385 (2020).
  * Fradkin (2013) E. Fradkin,  _Field theories of condensed matter physics_ (Cambridge University Press, 2013).
  * Naculich (1988) S. G. Naculich, Nuclear Physics B 296, 837 (1988).
  * Chandrasekharan (1994) S. Chandrasekharan, Physical Review D 49, 1980 (1994).
  * Dai  _et al._ (2016) H.-N. Dai, B. Yang, A. Reingruber, H. Sun, X.-F. Xu, Y.-A. Chen, Z.-S. Yuan, and J.-W. Pan, arXiv preprint arXiv:1602.05709 (2016).

##  Appendix A Details of the Self-Consistent Mean Field Theory Calculations

In this appendix we shall go over the details of the self-consistent mean
field theory approximations we used in the main body of the text. The method
we shall discuss is general, although the details of the calculation will