\section{Introduction}
Twisted bilayer graphene (TBG)~\cite{Cao2018,Cao2018b,Balents2020,Andrei2020}
forms an incommensurate moir\'e lattice. The description of its electronic
degrees of freedom, to a very good approximation, can be done using
concepts of periodic crystals, due to the weak scattering
between opposite valleys~\cite{LopesdosSantos2017,Trambly2010,Morell2010,Bistritzer2011}.
However, structurally, there is no periodicity in the system if the
twist angle takes a generic value. As was demonstrated
by Ochoa~\cite{Ochoa2019}, this has profound implications for the elasticity theory of TBG.
Phason modes $\vb{w}\mleft(\vb{x},t\mright)$, which correspond
to acoustic branches of the incommensurate lattice,
dominate lattice vibrations on the scale of the
moir\'e period and give rise to an additional twist stiffness $\kappa$ in the elastic energy
\begin{equation}
E_{\text{el}} \rightarrow E_{\text{el}} + \frac{\kappa}{8} \int \dd[2]{x} \, \mleft[\partial_{x} w_{y} - \partial_{y} w_{x}\mright]^{2}\,. \label{eq:stiffness}
\end{equation}
While in standard elasticity theory such a term is not allowed by rotational
symmetry~\cite{Chaikin2000}, in the case of TBG the adhesive potential between the two
layers gives rise to $\kappa > 0$.
Notice that the elasticity theory of planar quasicrystals~\cite{Levine1985,De1987,Ding1993}
can also be formulated in terms of Eq.~\eqref{eq:stiffness}.
The influence of the twist term and phason
excitations on electron-lattice couplings was discussed
in Refs.~\cite{Ochoa2019,Lian2019}, while their role
in the context of electronic nematicity was analyzed in Ref.~\cite{Fernandes2020}.

Phasons in TBG have a very straightforward interpretation.
Ignoring the structural relaxation of the lattice, they
can be identified with the relative displacement between layers~\cite{Ochoa2019}.
To illustrate a twist of the phason mode, we therefore
consider a moir\'e superlattice with sixfold rotational symmetry,
described in terms of the density profile
\begin{equation}
\varrho_{\theta}(\vb{r}) = \sum_{\vb{G}} \abs{\rho_{\vb{G}}} \Elr^{\iu \mleft(\vb{G} \vdot \vb{r} - \phi_{\vb{G}}\mright)}\,. \label{density}
\end{equation}
$\vb{G}$ are reciprocal lattice vectors of the moir\'e superlattice
and $\phi_{\vb{G}} = \vb{w} \vdot \vb{G}$.
In Fig.~\ref{fig:phason-density}, we show the density
profile with $\vb{w} = \theta \vu{z} \vcross \vb{r}$
for $\theta = 0$ and with an excited phason, i.e., for finite rotation of the twist angle $\theta$.
For simplicity, we consider only the leading harmonics, $\vb{G}_1$, $\vb{G}_2$,
and $\vb{G}_3 \equiv - \vb{G}_1 - \vb{G}_2$, with $\abs{\rho_{\vb{G}_i}} \equiv \rho$
and $\vb{G}_{1,2}$ denoting the standard primitive vectors of the triangular lattice.

The dynamics of phasons in incommensurate crystals also differs
from that of standard acoustic phonons by the fact that it is strongly
damped at long times~\cite{Lubensky1985,Baggioli2000}.
This is due to the friction that opposes the
relative motion of the incommensurate mass-density wave.
In usual elasticity, this damping does not occur
because the displacement couples to the generator of translations,
i.e.\ the momentum density. In TBG, however, anharmonic
terms of the adhesion potential between the layers give rise to a finite damping~\cite{Ochoa2021}.

\begin{figure}[t]
\includegraphics[scale=0.4]{Figure_density_v3.pdf}
\caption{Density profile for the triangular moir\'e superlattice, as given by Eq.~\eqref{density},
in the cases without rotation [panel (a), $\varrho_{\theta=0}(\vb{r})$]
and with a rotation of $\theta/(2 \pi) = 0.03$ [panel (b), $\varrho_{\theta}(\vb{r})$].
The phason excitation is illustrated by the difference in density profiles,
$\var{\varrho_{\theta}(\vb{r})} = \varrho_{\theta}(\vb{r}) - \varrho_{\theta=0}(\vb{r})$,
which is plotted in panel (c), and whose absolute value is plotted in panel (d).
The underlying triangular lattice is a guide to the eye.}
\label{fig:phason-density}
\end{figure}

In this paper, we discuss the impact of the twist
stiffness and phason damping on inevitable
defects of the moir\'e lattice such as dislocations and
disclinations. To this end,
we extend the recently formulated fracton formulation of standard
elasticity~\cite{Pretko2018,Pai2018,Pretko2019} to
incommensurate crystals with overdamped dynamics by developing an appropriate tensor gauge field theory.
The analogy to a tensor version of electromagnetism allows
us to draw general conclusions about the lattice defect dynamics of TBG.
The finite twist stiffness of Eq.~\eqref{eq:stiffness}
increases the degrees of freedom compared
to standard elasticity theory, giving rise to
a new type of defect that we call \emph{discompression}.
As a result, defects of the moir\'e lattice,
which can be described in terms of fractons, are governed by an additional conservation law
that dramatically affects the mobility of lattice defects.
In particular, we find that in an incommensurate crystal with phason excitations,
both glide and climb of dislocations are forbidden and
the dynamics of defects is always sub-diffusive. This
implies that dislocation networks
in twisted bilayer graphene are expected to be exceptionally stable.
We also show that our theory directly applies to other
incommensurate systems such as planar quasicrystals.

The impact of topological lattice defects on the mechanical and electronic properties of graphene and bilayer graphene has been widely discussed~\cite{Louie2010,Kirkland2012,Butz2013}. Lattice disorder has also been recognized as a key ingredient in TBG, with experiments reporting sharp variations of the twist angle and of uniaxial heterostrain over relatively short length scales, which are manifested in local elecronic properties~\cite{Zeldov2020,Bediako2021}. Here, we focus on the topological defects of the emergent moir\'e incommensurate lattice in TBG. While its elastic properties derive from those of the coupled graphene layers, we adopt a coarse-grained low-energy description, as in Ref.~\cite{Ochoa2019}, that does not require mapping the defects of the moir\'e lattice onto individual carbon atoms.

Defects of two-dimensional lattices can be efficiently captured in
terms of a dual description of elasticity theory, where disclinations
emerge as effective charges, dislocations as dipoles of these charges,
and elastic forces between them are transmitted by gauge
fields~\cite{Nelson1979,Young1979,Kleinert1989,Zaanen2004,Cvetkovic2007,Kleinert2008}.
This approach readily demonstrates why dislocations in commensurate lattices can
glide in the direction of the Burgers
vector $\vb{b}$, but not climb perpendicular
to it~\cite{Cvetkovic2007}.
While defects in aperiodic crystals were
discussed in Ref.~\cite{Bohsunng1987,Kleman2003,Kleman2013},
their dynamics remains an open problem.

A very elegant formulation of this duality
was recently achieved in Refs.~\cite{Pretko2018,Pai2018,Pretko2019,Radzihovsky2020a,Gromov2019,Gromov2020b,Nguyen2020}.
It was recognized that the dual formulation
of the usual elasticity theory is, in fact,
a fracton field theory. Fracton phases describe quantum phases of matter with
excitations of restricted mobility~\cite{Nandkishore_review}.
If considered in isolation, a
fracton is immobile, either along certain directions or as a whole.
It can only move collectively, via interactions with other degrees
of freedom. Fractons were initially discussed in the context of non-ergodic
quantum glass models~\cite{Chamon2005} and stabilizer codes for self-correcting
quantum memory~\cite{Haah2011}. More recently, it was recognized that
fractons can be efficiently described in terms of tensor gauge theories~\cite{Pretko2017}.
The restricted mobility of fractons emerges in terms of dipole conservation
laws and gives rise to anomalous, sub-diffusive
hydrodynamics~\cite{Gromov2020,Feldmeier2020,Moudgalya2021,Radzihovsky2020}.

\section{Elasticity theory of twisted bilayer graphene crystals}
We start with the action of elasticity
theory~\cite{Kleinert1989,Chaikin2000},
supplemented by the twist stiffness of Eq.~\eqref{eq:stiffness}.
To include dissipative dynamics, we use
the Schwinger-Keldysh formalism~\cite{Kamenev2011}:
\begin{align}
\begin{aligned}
S &= \frac{1}{2} \int_{\vb{x},t,t' \in \mathcal{C}} w_{i}(\vb{x},t) D^{-1}\mleft(t-t'\mright) w_{i}(\vb{x},t')  \\
&- \frac{1}{2} \int_{\vb{x},t \in \mathcal{C}} \mleft[C_{ij,k\ell} w_{ij}(\vb{x},t) w_{k\ell}(\vb{x},t) + \kappa \vartheta_{w}^2(\vb{x},t)\mright]\,. \label{eq:action}
\end{aligned}
\end{align}
Here, $\int_{\vb{x},t \in \mathcal{C}} = \int_{\mathcal{C}} \dd{t} \int_{V} \dd[2]{x}$
indicates that the time integration has to be performed
on the round-trip Keldysh contour, while the spatial
integration goes over the crystal volume $V$.
$w_{ij} \equiv \frac{1}{2} \mleft(\partial_i w_{j} + \partial_j w_{i}\mright)$
is the symmetric phason strain, the elastic constants $C_{ij,k\ell}$ are
summarized in Ref.~\cite{Ochoa2019}, and the bond angle variation
is $\vartheta_{w} \equiv \frac{1}{2} \LCs_{ij} \partial_i w_{j}$.
Damping enters through the nonlocal-in-time contribution
\begin{equation}
 D^{-1}\mleft(t-t'\mright) = - \Dd\mleft(t-t'\mright) \partial_{t'}^{2} + \gamma\mleft(t-t'\mright)\,,
\end{equation}
where we assume Ohmic damping $\Im \gamma^{R}(\omega) = \gamma_{0} \omega$
for the retarded self energy. Without damping, the equation of motion
become elastic wave equations, while the dynamics
becomes diffusive at large $\gamma_0$. For further details on the Keldysh formlism, see Appendix~\ref{app:A}.

We shall use a notation that will prove efficient in 2D.
In it, we always start with lower index vectors or tensors $A_{i}$.
The raising of their indices we define as a contraction
with the Levi-Civita symbol:
\begin{equation}
A^{i} \equiv \LCs_{ij} A_{j}\,.
\end{equation}
From $A_{i} = \mleft(A_{x}, A_{y}\mright)$ and $A^{i} \equiv \LCs_{ij} A_{j}$
we see that the upper index vector $A^{i} = \mleft(A_{y}, - A_{x}\mright)$
is the lower index vector rotated in the clockwise direction by $\uppi/2$. This
ensures that $\vb{A} \vdot \vb{B} = A_{i} B_{i} = A^{i} B^{i}$.
On the other hand,
$\mleft(\vb{A} \vcross \vb{B}\mright)_{z} = A_{i} B^{i} = - A^{i} B_{i}$
and in particular $A_i A^i = 0$. The divergence and two-dimensional curl
may likewise be expressed in this notation as
\begin{align}
\grad \vdot \vb{A} &= \partial_{i} A_{i}\,, \nonumber \\
\mleft(\grad \vcross \vb{A}\mright)_{z} &= \partial_{i} A^{i} = - \partial^{i} A_{i}\,.
\end{align}
For single-valued functions $f$, $\partial_{i} \partial^{i} f = 0$. For
multivalued functions, the partial derivative do not commute at the branch cuts.
The most important benefit of this notation, reminiscent of the
van der Waerden spinor notation~\cite{vdWaerden1928},
is that the correspondence between elasticity and tensor
electromagnetism~\cite{Pretko2018,Pai2018,Pretko2019,Radzihovsky2020a}
amounts to simply raising all the indices.

\subsection{Dual gauge theory}

Building on the approach of Refs.~\cite{Pretko2018,Pai2018,Pretko2019},
the first step in formulating a dual description is to
introduce the fields
$\partial_t w_{i} \to \pi_i$, $C_{ij,k\ell} w_{k\ell} \to \sigma_{ij}$,
and $\kappa \vartheta_{w} \to M$
through a series of Hubbard-Stratonovich transformations.
They can be associated with momentum, stress, and torque, respectively.
This yields the real-time action
\begin{align}
S &= \frac{1}{2} \int_{\vb{x}, t \in \mathcal{C}} \mleft[C_{ij,k\ell}^{-1} \sigma_{ij}(\vb{x},t) \sigma_{k\ell}(\vb{x},t) + \kappa^{-1} M^{2}(\vb{x},t)\mright] \nonumber \\
&\quad - \frac{1}{2} \int_{\vb{x},t,t' \in \mathcal{C}} \pi_{i}(\vb{x},t) \Gamma\mleft(t-t'\mright)\pi_{i}(\vb{x},t') + S_{w}\,, \label{afterHStransf}
\end{align}
where the retarded  version of $\Gamma$ is
$\Gamma^{R}(\omega) = \omega^{2} / \mleft(\omega^{2} + \iu \gamma_{0} \omega\mright)$.
The last term
\begin{equation}
S_{w} = \int_{\vb{x},t \in \mathcal{C}} \mleft[\pi_{i} \mleft(\partial_{t} w_{i}\mright) - \sigma_{ij} w_{ij} - M \vartheta_{w}\mright] \label{eq:stress-strain-term}
\end{equation}
still contains the original phason field.

The next step
is to decompose $w_{i} = \tilde{w}_{i} + w_{i}^{(s)}$
into a smooth part $\tilde{w}_{i}$
and a singular part $w_{i}^{(s)}$ due to topological defects.
After integrating out the smooth part,
we obtain the constraint
\begin{equation}
\partial_t \pi_j = \partial_i \Sigma_{ij}\,, \label{constraint}
\end{equation}
which resembles Newton's second law with a non-symmetric stress
\begin{equation}
\Sigma_{ij}=\sigma_{ij} + \frac{1}{2} \LCs_{ij} M\,,
\end{equation}
consisting of
the symmetric stress tensor $\sigma_{ij}$ and the torque $M$.
Whereas $\sigma_{ij}$ also appears in ordinary crystals,
the torque only arises due
to the finite twist stiffness.

The last step is to introduce (tensorial) gauge potentials that enforce
the constraint~\eqref{constraint}:
\begin{align}
\begin{aligned}
\Sigma_{ij} &= - \partial_t A_{ij} - \partial^i \phi_j\,, \\
\pi_j &= - \partial_i A_{ij}\,.
\end{aligned}
\end{align}
These gauge potentials are invariant under the gauge transformations
\begin{align}
\begin{aligned}
A_{ij} &\mapsto A_{ij} + \partial^{i} \Lambda_{j}\,, \\
\phi_{j} &\mapsto \phi_{j} - \partial_t \Lambda_{j}\,.
\end{aligned} \label{gauge_transformation}
\end{align}
The stress-strain coupling Eq.~\eqref{eq:stress-strain-term} now takes
the familiar form of a charge/current-potential coupling:
\begin{equation}
S_{w} = \int_{\vb{x},t \in \mathcal{C}} \mleft(A_{ij} J_{ij} - \phi_{j} Q_{j}\mright)\,. \label{lagrangian_gauge}
\end{equation}
Due to the tensorial nature of the gauge fields,
the charge density
\begin{equation}
Q_{j} = \partial^{i} \partial_{i} w_{j}^{(s)}
\end{equation}
is a vector
and the corresponding current density,
\begin{equation}
J_{ij} = \mleft(\partial_i \partial_t - \partial_t \partial_i\mright) w_{j}^{(s)}\,,
\end{equation}
a tensor.
Demanding the invariance of~\eqref{lagrangian_gauge} under
the gauge transformations
yields the continuity equation:
\begin{equation}
\partial_{t} Q_{j} + \partial^{i} J_{ij} = 0\,. \label{continuity_p}
\end{equation}

\subsection{Electromagnetic analogy}
Let us develop, in analogy to Ref.~\cite{Pretko2018,Pai2018,Pretko2019,Radzihovsky2020a} for standard elasticity,
a physical illuminating electromagnetic analogy of
the conserved vector charge $Q_j$. Consider a single
charge
\begin{equation}
Q_{j}(\vb{x},t) = q_{j} \Dd\mleft(\vb{x}-\vb{r}(t)\mright)
\end{equation}
at
position $\vb{r}(t)$ that moves with velocity $\vb{v}(t) = \dot{\vb{r}}(t)$.
Eq.~\eqref{continuity_p} is then fulfilled by
\begin{equation}
J_{ij}(\vb{x},t) = v^{i} q_{j} \Dd\mleft(\vb{x}-\vb{r}(t)\mright)\,.
\end{equation}
Inserting this into the action~\eqref{lagrangian_gauge}
yields the Lorentz force
\begin{equation}
F_{i} = \var{S_{w}} / \var{r_{i}} = - \mleft(\Sigma^{i}_{~j} + v^{i} \pi_{j}\mright) q_{j}\,.
\end{equation}
The problem behaves like the electrodynamics
of non-symmetric tensor electric fields $ \Sigma^{i}_{~j}$
and vector magnetic fields $\pi_{i}$.
This  analogy also helps us
understand why the overdamped dynamics of the phasons does not spoil
the entire gauge description. The dissipative propagator $\Gamma(t-t')$ for
$\pi_i(t)$ in Eq.~\eqref{afterHStransf} translates,
in standard electromagnetism, to a frequency-dependent
magnetic permeability $\Gamma(\omega)$, which clearly
does not violate the gauge description of electromagnetism.

Besides the effects of dissipation, a crucial
new aspect of the incommensurate moir\'e lattice elasticity
is the vector continuity equation~\eqref{continuity_p}
that can be rephrased in terms of \emph{two} conserved scalar densities:
\begin{align}
\begin{aligned}
\rho^{(\ell)} &= \partial_{j} Q_{j}\,, \\
\rho^{(t)} &= \partial^{j} Q_{j}\,,
\end{aligned}
\end{align}
with continuity equations
\begin{align}
\begin{aligned}
\partial_{t} \rho^{(\ell)} + \partial_{j} \partial^{i} J_{ij} &= 0\,, \\
\partial_{t} \rho^{(t)} + \partial^{j} \partial^{i} J_{ij} &= 0\,.
\end{aligned} \label{continuity_div_p}
\end{align}
The transverse density $\rho^{(t)}$ also appears
in standard elasticity, where it was identified
as the density of disclinations, with $Q_i$ denoting
the Burgers vector density~\cite{Zaanen2004,Pretko2018}.
The longitudinal density $\rho^{(\ell)}$
is the new conserved charge present only in incommensurate crystals.

\subsection{Interpretation of the charges}
The charges of the elastic gauge theory
are topological lattice defects that make
the displacement field $w_{i}$ multivalued~\cite{Kleinert2008}.
We can write the singular part of a dislocation displacement field as
\begin{equation}
\vb{w}^{(s,Q)}(\vb{x}) = \frac{\Im \Log(z)}{2 \uppi} \, \vb{b}\,, \label{us_disloc}
\end{equation}
where $z = x_1 + \iu x_2$, $\Log(z)$ is the principal branch of the logarithm
with a branch cut along the negative $x_1$ axis,
and $\vb{b} = b_1 \vu{e}_1 + b_2 \vu{e}_2$ is the Burgers vector.
At the origin the partial derivatives do not commute, yielding
$Q_j(\vb{x}) = b_j \Dd(\vb{x})$, which confirms that the
vector charge is the Burgers vector density.

Next, we construct a disclination from a Dirac string of dislocations.
By integrating $\int_{-\infty}^{0} \dd{x_1'} \, \vb{w}^{(s,Q)}(x_1-x_1', x_2)$,
ignoring all except the singular
parts, and setting $b_1=0$ and
$b_2 \neq 0$, we obtain
\begin{equation}
\vb{w}^{(s,t)}(\vb{x}) = - b_2 \frac{\Im\mleft[z \Log(z)\mright]}{2 \uppi} \, \vu{e}_2\,, \label{us_discli}
\end{equation}
with $Q_j = b_2 \Kd_{2j} \Theta(-x_1) \Dd(x_2)$,
$\rho^{(\ell)} = b_2 \Theta(-x_1) \Dd'(x_2)$, and
$\rho^{(t)} = b_2 \Dd(\vb{x})$.
This confirms that $\rho^{(t)}$ is the density of disclinations.
In addition, the result for $\rho^{(l)}$ implies that a disclination
also corresponds to a Dirac string along the negative $x_1$-axis
made of dipoles ($\sim \Dd'(x_2)$) of the longitudinal charge.

Similarly, we may consider a Dirac string of dislocations
whose Burgers vectors point along the branch cut,
instead of perpendicular to it
($b_1 \neq 0$ and $b_2 = 0$),
\begin{equation}
\vb{w}^{(s,\ell)}(\vb{x}) = - b_1 \frac{\Im\mleft[z \Log(z)\mright]}{2 \uppi} \, \vu{e}_1\,. \label{us_discompress}
\end{equation}
The charges are $Q_j = b_1 \Kd_{1j} \Theta(-x_1) \Dd(x_2)$,
$\rho^{(\ell)} = - b_1 \Dd(\vb{x})$, and
$\rho^{(t)} = b_1 \Theta(-x_1) \Dd'(x_2)$.
Since $\vb{w}^{(s,\ell)}$ points parallel to the Dirac string,
and its jump in value is proportional to the distance from
the origin, it represents an abrupt change in the strain.
We shall thus call this defect a \emph{discompression}.
Below Eq.~\eqref{us_discli} we showed that a disclination
is equivalent to a Dirac string of discompression dipoles.
The result for $\rho^{(t)}$ implies, in turn, that a discompression
equals a Dirac string of disclination dipoles. These ``fusion rules'' for the
various defects are sketched in the lower part of Fig.~\ref{fig:defects}.

In Fig.~\ref{fig:defects}, we show these defects for a triangular lattice.
Disclinations and discompressions are energetically
very expensive as both come along with macroscopic
regions of compression or bond-angle mismatch.
More formally, both are forbidden as single
charges by an appropriate generalization of Weingarten's theorem~\cite{Kleinert2008}.
Hence, a moir\'e lattice is charge neutral
with regards to $\rho^{(t,\ell)}$, while their
dipoles and higher moments remain perfectly valid excitations.
Fig.~\ref{fig:defects}(c) illustrates the known fact that a dislocation
corresponds to a dipole of disclinations (panels (a) and (b))
with dipole moment perpendicular to the Burgers vector $\vb{b}$.
Our analysis shows that a dislocation [Fig.~\ref{fig:defects}(f)] can
also emerge from two discompressions of opposite charge (panels (d) and (e)).
The dipole moment is now oriented along the Burgers vector.
Single discompressions of Fig.~\ref{fig:defects}(d) and (e) can be
generated in terms of a Volterra process
where we cut along the negative $x_1$-axis and displace matter
above and below the cut in opposite directions parallel to the cut,
before glueing back together. This is distinct from disclinations
of Fig.~\ref{fig:defects}(a) and (b) where the displacement
is perpendicular to the cut and opens a wedge.

\begin{figure}[t]
\centering
\includegraphics[scale=0.8]{Figure_defects_columns_v1.pdf}
\caption{Various defects for a triangular lattice and the relations between them.
(a) Positive disclination ($b_2 = 2 / \sqrt{3}$).
(b) Negative disclination ($b_2 = - 2 / \sqrt{3}$).
(c) Dislocation as a dipole of disclinations.
(d) Positive discompression ($b_1 = 1/2$).
(e) Negative discompression ($b_1 = -1/2$).
(f) Dislocation as a dipole of discompressions.
The table summarizes the fusion rules that explain how to represent charges in terms of Dirac strings of dipoles of other charges.}
\label{fig:defects}
\end{figure}

\subsection{Consequences of the continuity equations}
Eq.~\eqref{continuity_p} implies the
conservation of the total Burgers vector.
Let us analyze this vector conservation law in terms of
the continuity equations~\eqref{continuity_div_p} of the scalar densities
that both contain two spatial derivatives. This implies, in addition to the conservation
of the total disclination and discompression charges,
that both of their total dipole moments are conserved as well:
\begin{equation}
\frac{\dd{}}{\dd{t}} \int \dd[2]{x} \  x_k \, \rho^{(\ell,t)} = 0\,. \label{conservation_disclin_di}
\end{equation}
Both charges are therefore immobile fractons.
Only correlated movements at fixed dipole density are allowed.

\subsection{Energetics and equations of motion}
Thus far we analyzed conservation laws that follow from
the gauge invariance of the dual theory. To understand the
actual dynamics of  defects requires, however,  an analysis of
the energetics of vector charges. To this end we determine
the equations of motion of the gauge fields that
mediate the stresses. The Euler-Lagrange equation of $A_{ij}$ is:
\begin{equation}
\partial_{t} \tilde{\Sigma}_{ij}(t) = \int \dd{t'} \, \Gamma^{R}\mleft(t-t'\mright) \partial_{i}\pi_{j}(t') - J_{ij}(t)\,, \label{eq:A_ij-EL}
\end{equation}
with $\tilde{\Sigma}_{ij} = C_{ij,k\ell}^{-1} \sigma_{k\ell} + \LCs_{ij} \kappa^{-1} M$.
It is the analog of Ampere's law in electrodynamics.
The analog of Gauss's law $\partial^{i} \tilde{\Sigma}_{ij} = Q_{j}$
follows from the variation of the potential $\phi_j$ and demonstrates
that $Q_j$ is indeed the source of stress.

We can use Ampere's law
to determine the dynamics of the quadrupolar moments.
To achieve this, we first contract the
Euler-Lagrange equation~\eqref{eq:A_ij-EL} with $\Kd_{ij}$
and $\LCs_{ij}$. This gives
\begin{align}
\partial_t C_{ii,k\ell}^{-1} \sigma_{k\ell} &= \partial_i (\Gamma^{R} \pi_i) - J_{ii}\,, \label{contr_delta} \\
\frac{2}{\kappa} \partial_t M &= \partial_i (\Gamma^{R} \pi^i) - J^{~i}_i\,, \label{contr_epsilon}
\end{align}
where
\begin{equation}
\mleft(\Gamma^{R} \pi_j\mright)(\vb{x},t) = \int \dd{t'} \, \Gamma^{R}\mleft(t - t'\mright) \pi_j(\vb{x},t')\,.
\end{equation}

Now we use the conservation laws~\eqref{continuity_div_p},
partially integrate the quadrupolar moments twice, and then
use the equations of motion. This yields
\begin{align}
\frac{\dd{}}{\dd{t}} \int_V \abs{\vb{x}}^2 \rho^{(t)} ~ &\overset{\mathclap{\eqref{continuity_div_p}}}{=} ~ - \int_V x_k x_k \partial^{j} \partial^{i} J_{ij} \label{sup:quad-der1} \\
&= - 2 \int_V J_{ii} + \int_V \partial_{j} \mleft(2 x_i J^{ji} - \abs{\vb{x}}^2 \partial_{i} J^{ij}\mright) \notag \\
&\overset{\mathclap{\eqref{contr_delta}}}{=} ~ 2 \int_V \partial_t C_{ii,k\ell}^{-1} \sigma_{k\ell} \notag \\
&+ \int_V \partial_{j} \mleft(2 x_i J^{ji} - \abs{\vb{x}}^2 \partial_{i} J^{ij} - 2 \Gamma^R \pi_{j}\mright) \notag
\end{align}
and
\begin{align}
\frac{\dd{}}{\dd{t}} \int_V \abs{\vb{x}}^2 \rho^{(\ell)} \label{sup:quad-der2} ~ &\overset{\mathclap{\eqref{continuity_div_p}}}{=} ~ - \int_V x_k x_k \partial_{j} \partial^{i} J_{ij} \\
&= - 2 \int_V J^{~i}_{i} + \int_V \partial_{j} \mleft(\abs{\vb{x}}^2 \partial_{i} J^{i}_{~\,j} - 2 x_i J^{j}_{~\,i}\mright) \notag \\
&\overset{\mathclap{\eqref{contr_epsilon}}}{=} ~ \frac{4}{\kappa} \int_V \partial_t M \notag \\
&+ \int_V \partial_{j} \mleft(\abs{\vb{x}}^2 \partial_{i} J^{i}_{~\,j} - 2 x_i J^{j}_{~\,i} - 2 \Gamma^R \pi^{j}\mright)\,. \notag
\end{align}
Finally, we rewrite the sources by exploiting the
Euler-Lagrange equations of the action~\eqref{afterHStransf},
$C_{ij,k\ell}^{-1} \sigma_{k\ell} = w_{ij}$ and $M = \kappa \vartheta_{w}$,
and thus obtain
\begin{equation}
\frac{\dd{}}{\dd{t}} \int \dd[2]{x} \, \abs{\vb{x}}^2 \rho^{(\ell,t)} = \frac{\dd{}}{\dd{t}} S^{(\ell,t)}\,. \label{conservation_quad}
\end{equation}
The source terms for disclination quadrupoles is
\begin{equation}
    S^{(t)} = 2 \int \dd[2]{x} \, \partial_i w_{i}
\end{equation}
and corresponds to a costly volume change of the system. In turn,
\begin{equation}
    S^{(\ell)} = 4 \int \dd[2]{x} \, \vartheta_{w}
\end{equation}
for discompression
quadrupoles corresponds to a global twist.
These quadrupoles can easily be
interpreted if one considers an isolated dislocation $Q_j(\vb{x}) = b_j \Dd\mleft(\vb{x} - \vb{r}(t)\mright)$ whose $\int \dd[2]{x} \, \abs{\vb{x}}^2 \rho^{(t)} = 2 b^j r_j(t)$
and $\int \dd[2]{x} \, \abs{\vb{x}}^2 \rho^{(\ell)} = - 2 b_j r_j(t)$.
Thus the quadrupolar moments are the components of the dislocation position
parallel and perpendicular to the Burgers vector.
Eq.~\eqref{conservation_quad} was obtained up to boundary terms. The boundary terms vanish or, more physically, cannot be
changed by local operations that act in the bulk of the medium.
Thus~\eqref{conservation_quad} expresses how
in the bulk the sources of changing
quadrupolar moments are the energetically costly
compression and disorientation of the crystal.

In our analysis we freely used expressions
that are valid only on-shell, the various conservation laws (\ref{continuity_p}, \ref{continuity_div_p}, \ref{conservation_disclin_di}, \ref{conservation_quad}),
equations of motion (\ref{eq:eomdisplfull}, \ref{eq:A_ij-EL}, $C_{ij,k\ell}^{-1} \sigma_{k\ell} = \partial_i w_{j}$, $M = \kappa \vartheta_{w}$),
expressions for the sources $S^{(t)}$ and  $S^{(\ell)}$,
etc., can all be formulated as rigorous statements
that apply to the respective path-integral averages.
This is proved by using the Schwinger-Dyson equations.

From Eq.~\eqref{conservation_quad} we may therefore conclude that a climb, i.e.\ the
motion perpendicular to $b_j$ of a dislocation,
is accompanied by an energetically costly change in the volume of the crystal.
In addition, a glide, i.e.\ the motion parallel to $b_j$ of a dislocation,
is accompanied by an energetically costly change of the orientation of the crystal.
While the first result holds in any crystal,
the second is due to the non-zero twist stiffness
of incommensurate ones. In its presence the motion of dislocations is forbidden
entirely instead of just being restricted to one direction.

\subsection{Sub-diffusive hydrodynamics}

Thus far we have analyzed the motion of isolated dipoles and quadrupoles. The situation
becomes more complicated for higher-order multipoles. However, using the hydrodynamic
description of fractons~\cite{Gromov2020,Feldmeier2020,Moudgalya2021,Radzihovsky2020} one readily finds that there
cannot be a Fick's law where gradients of the vector
charges yield currents (i.e.\ $J_{ij} \propto \partial_i Q_j$ or $\propto \partial^i Q_j$),
leading to ordinary diffusion of all charges. If this were true, the mere existence
of scalar charges would induce a current $J_{ii} \sim \rho^{(\ell,t)}$, which is not allowed.
Expanding in gradients, the leading
symmetry-allowed term connects $J_{ij}$ with
a third spatial derivative of $Q_{j}$.
This gives rise to a coupled dynamics of the two
scalar densities that form two sub-diffusive
modes with dispersions $\omega = \iu B_{\pm} k^4$.

\section{Generalization to generic planar quasicrystals and incommensurate crystals}
\label{sec:III}
Finally, we demonstrate that our results are not restricted to TBG and
apply directly to planar quasicrystals and incommensurate crystals with two length scales.
To this end  we relate the elastic gauge theory with twisting term
considered in the main text
to the elasticity theory of two-dimensional quasicrystals.
In these systems, the low-energy elastic variables are
doubled compared to periodic crystals~\cite{Levine1985,De1987,Ding1993}.
In addition to the phonons $\vb{u}$, there are phasons $\vb{w}$ that can
be shown to have a twist stiffness like in Eq.~\eqref{eq:stiffness}.
The defects of the system are characterized
by two Burgers vectors $\vb{b}_{u} = \oint \dd{\vb{u}}$ and $\vb{b}_{w} = \oint \dd{\vb{w}}$~\cite{De1987}.
Our results imply that pure displacement dislocations can glide,
whereas those that involve $\vb{b}_{w}$ cannot.
This remains true even in the presence of phonon-phason coupling.

A symmetry analysis of the problem yields the following result for
the elastic energy~\cite{Levine1985,De1987,Ding1993}:
\begin{equation}
\mathcal{E}_{\text{el}} = \frac{1}{2} C_{ij,k\ell} u_{ij} u_{k\ell} + R_{ij,k\ell} u_{ij} w_{k\ell} + \frac{1}{2} K_{ij,k\ell} w_{ij} w_{k\ell}\,,
\end{equation}
where $u_{ij} = \frac{1}{2} \mleft(\partial_{i} u_{j} + \partial_{j} u_{i}\mright)$
is the usual symmetric strain tensor describing phonons and $w_{ij} = \partial_{i} w_{j}$
is the non-symmetric phason tensor.

For the moment, let us focus on the case of a twelvefold symmetric planar quasicrystals
that has the special property of no phonon-phason mixing, i.e., $R_{ij,k\ell} = 0$~\cite{Ding1993}.
Because of the high rotational symmetry, the elastic constants $C_{ij,k\ell}$
have the same form as those of isotropic media:
\begin{equation}
C_{ij,k\ell} = \lambda \, \Kd_{ij} \Kd_{k\ell} + \mu \, \mleft(\Kd_{ik} \Kd_{j\ell} + \Kd_{i\ell} \Kd_{jk}\mright)\,. \label{sup:eq:C_ijkl}
\end{equation}
On the other hand,
\begin{equation}
K_{ij,k\ell} = K_{1} \Kd_{ik} \Kd_{j\ell} + K_{2} \mleft(\Kd_{ij} \Kd_{k\ell} - \Kd_{i\ell} \Kd_{jk}\mright) + \tilde{K}_{3} \LCs^{S}_{ij} \LCs^{S}_{k\ell}\,, \label{sup:eq:K_ijkl}
\end{equation}
where $\LCs^{S}_{11} = \LCs^{S}_{22} = 0$ and $\LCs^{S}_{12} = \LCs^{S}_{21} = 1$.
We note that $K_{ij,k\ell}$ satisfies only the major index symmetry
$K_{ij,k\ell} = K_{k\ell,ij}$, but not the minor ones
$K_{ij,k\ell} \neq K_{ji,k\ell}$ and $K_{ij,k\ell} \neq K_{ij,\ell k}$.

To make contact with our system~\eqref{eq:action},
we merely have to recognize that $K_{ij,k\ell}$
can be rewritten as
\begin{equation}
K_{ij,k\ell} = \tilde{C}_{ij,k\ell} + \frac{\kappa}{4} \LCs_{ij} \LCs_{k\ell}\,, \label{sup:decomposition}
\end{equation}
where
\begin{equation}
\tilde{C}_{ij,k\ell} = \tilde{\lambda} \, \Kd_{ij} \Kd_{k\ell} + \tilde{\mu} \, \mleft(\Kd_{ik} \Kd_{j\ell} + \Kd_{i\ell} \Kd_{jk}\mright) + \tilde{K}_{3} \LCs^{S}_{ij} \LCs^{S}_{k\ell}
\end{equation}
is the minor-index-symmetric part that
satisfies $\tilde{C}_{ij,k\ell} = \tilde{C}_{ji,k\ell} = \tilde{C}_{ij,\ell k}$
and has $\tilde{\lambda} = K_{2}$ and $\tilde{\mu} = \frac{1}{2} \mleft(K_{1} - K_{2}\mright)$.
The minor-index-antisymmetric part plays the role of a twisting term with stiffness
\begin{equation}
\kappa = 2 \mleft(K_{1} + K_{2}\mright)\,. \label{sup:eff-kappa}
\end{equation}
Thus the phonon displacement field $u_{i}$ behaves like in the elasticity theory
of a conventional crystal, while the phason $w_{i}$ is governed by an elastic energy identical
to the one of twisted bilayer graphene, given in Eq.~\eqref{eq:stiffness} or~\eqref{eq:action}.
Hence all our conclusions about the mobility of defects of the phason modes carry over.

The situation is more complicated for quasicrystals with fivefold or eightfold symmetry.
In these cases, $C_{ij,k\ell}$ is still isotropic, as in~\eqref{sup:eq:C_ijkl}, and
$K_{ij,k\ell}$ still has the form~\eqref{sup:eq:K_ijkl},
although with $\tilde{K}_3 = 0$ in the case of fivefold symmetry,
and can therefore be decomposed
according to~\eqref{sup:decomposition}. Most significantly,
a phonon-phason coupling of the form
\begin{equation}
R_{ij,k\ell} = R \mleft(\Kd_{i1} - \Kd_{i2}\mright) \mleft[\Kd_{ij} \Kd_{k\ell} + \Kd_{ik} \Kd_{j\ell} - \Kd_{i\ell} \Kd_{jk}\mright]
\end{equation}
is allowed. Consequently, one has to redo the whole gauge formulation.

The end results, however, turn out to be insensitive to phonon-phason coupling.
Specifically, for the phonon field one finds
the usual results~\cite{Pretko2018,Pai2018,Pretko2019}:
\begin{align}
\partial_{t} \rho^{(t)}_{u} + \partial^{j} \partial^{i} J_{uij}^{S} &= 0\,, \\
\frac{\dd{}}{\dd{t}} \int \dd[2]{x} \  x_k \, \rho^{(t)}_{u} &= 0\,, \\
\frac{\dd{}}{\dd{t}} \int \dd[2]{x} \, \abs{\vb{x}}^2 \rho^{(t)}_{u} &= \frac{\dd{}}{\dd{t}} S^{(t)}_{u}\,,
\end{align}
where $Q_{uj} = \partial^{i} \partial_{i} u_{j}^{(s)}$,
$J_{uij} = \mleft(\partial_i \partial_t - \partial_t \partial_i\mright) u_{j}^{(s)}$,
$\rho^{(t)}_{u} = \partial^{j} Q_{uj}$,
$J_{uij}^{S} = \frac{1}{2} \mleft(J_{uij} + J_{uji}\mright)$, and
$S^{(t)}_{u} = 2 \int \dd[2]{x} \, \partial_i u_i$. For the
phason fields one similarly finds that \eqref{continuity_p},
\eqref{continuity_div_p}, \eqref{conservation_disclin_di},
and \eqref{conservation_quad} continue to hold, but with the
appropriate $Q_{wj} = \partial^{i} \partial_{i} w_{j}^{(s)}$,
$J_{wij} = \mleft(\partial_i \partial_t - \partial_t \partial_i\mright) w_j^{(s)}$,
$\rho^{(\ell)}_{w} = \partial_{j} Q_{wj}$, $\rho^{(t)}_{w} = \partial^{j} Q_{wj}$,
$S^{(t)}_{w} = 2 \int \dd[2]{x} \, \partial_i w_i$, and
$S^{(\ell)}_{w} = 4 \int \dd[2]{x} \, \vartheta_w$.
Therefore, the glide of defects that have a finite
phason Burgers vectors $\vb{b}_{w} = \oint \dd{\vb{w}}$
is suppressed, and the parameter that controls the supression
is the effective twisting stiffness~\eqref{sup:eff-kappa}.

\emph{Details of the phonon-phason coupled gauge theory:}
It is convenient to introduce an additional index $\mu,\nu \in \{u,w\}$
that differentiates the phonon and phason displacement fields. That way,
$u_{ui} = u_{i}$, $u_{wi} = w_{i}$, $u_{uij} = u_{ij}$, $u_{wij} = w_{ij}$, and
\begin{gather}
C_{uij,uk\ell} = C_{ij,k\ell}\,, \quad C_{wij,wk\ell} = K_{ij,k\ell}\,, \\
C_{uij,wk\ell} = C_{wk\ell,uij} = R_{ij,k\ell}\,.
\end{gather}
The appropriate generalization of the action~\eqref{eq:action} we may now write as
\begin{align}
\begin{aligned}
S &= \frac{1}{2} \int_{\vb{x},t,t' \in \mathcal{C}} u_{\mu i}(t) D^{-1}_{\mu}\mleft(t-t'\mright) u_{\mu i}(t') \\
&\qquad - \frac{1}{2} \int_{\vb{x},t \in \mathcal{C}} C_{\mu ij,\nu k\ell} u_{\mu ij} u_{\nu k\ell}\,.
\end{aligned}
\end{align}
Next, we introduce the Hubbard-Stratonovich fields
$\partial_t u_{\mu i} \to \pi_{\mu i}$ and $C_{\mu ij,\nu k\ell} u_{\nu k\ell} \to \Sigma_{\mu ij}$,
integrate out the regular parts of $u_{\mu i}$, and then enforce the
constraints $\partial_t \pi_{\mu j} = \partial_i \Sigma_{\mu ij}$
through gauge potentials:
\begin{align}
\begin{aligned}
\Sigma_{uij} &= - \partial_t A_{uij}^{S} - \partial^{i} \partial^{j} \phi_{u}\,, \\
\pi_{uj} &= - \partial_{i} A_{uij}^{S}\,,
\end{aligned} \\
\begin{aligned}
\Sigma_{wij} &= - \partial_t A_{wij} - \partial^{i} \phi_{wj}\,, \\
\pi_{wj} &= - \partial_{i} A_{wij}\,,
\end{aligned}
\end{align}
where $A_{uij}^{S} = A_{uji}^{S}$ is symmetric. This yields the action
\begin{align}
\begin{aligned}
S &= - \frac{1}{2} \int_{\vb{x},t,t' \in \mathcal{C}} \pi_{\mu i}(t) \Gamma_{\mu}\mleft(t-t'\mright)\pi_{\mu i}(t') \\
&\qquad + \frac{1}{2} \int_{\vb{x}, t \in \mathcal{C}} C_{\mu ij,\nu k\ell}^{-1} \Sigma_{\mu ij} \Sigma_{\nu k\ell} + S_{uw}\,,
\end{aligned} \label{sup:afterHStransf-2}
\end{align}
where $\Gamma_{\mu}\mleft(t-t'\mright) = \lvec{\partial_{t}} D_{\mu}\mleft(t-t'\mright) \rvec{\partial_{t'}}$ and
\begin{align}
\begin{aligned}
S_{uw} = \int_{\vb{x},t \in \mathcal{C}} \big[&A_{uij}^{S} J_{uij}^{S} + \phi_{u} \rho^{(t)}_{u} \\
&+ A_{vij} J_{vij} - \phi_{vj} Q_{vj}\big]\,,
\end{aligned}
\label{sup:charge-potential-coupling}
\end{align}
with the previously defined charges and current densities.

In this action we see why the phonon-phason coupling does
not affect the argument leading to the quadrupolar
conservation laws. On the one hand, the charge conservation
laws are a consequence of the local gauge symmetry and
from the form of the source term~\eqref{sup:charge-potential-coupling}
we immediately see that they are unaffected by $C_{\mu ij,\nu k\ell}^{-1}$
or $R_{ij,k\ell}$. On the other hand, let us
take a look at the Euler-Lagrange
equations of $A_{\mu ij}$ that allowed for an additional
partial integration in~\eqref{sup:quad-der1} and~\eqref{sup:quad-der2}:
\begin{equation}
\partial_{t} \tilde{\Sigma}_{\mu ij}(t) = \int \dd{t'} \, \Gamma^{R}_{\mu}\mleft(t-t'\mright) \partial_{i}\pi_{\mu j}(t') - J_{\mu ij}(t)\,.
\end{equation}
Once one expresses $\tilde{\Sigma}_{\mu ij} = C_{\mu ij,\nu k\ell}^{-1} \Sigma_{\nu k\ell}$
in terms of the original fields $u_{\mu ij}$,
one finds that $\tilde{\Sigma}_{\mu ij} = u_{\mu ij}$
according to the Euler-Lagrange equations of $\Sigma_{\mu ij}$.
Thus $C_{\mu ij,\nu k\ell}$ again drop out of the argument.
Therefore, $C_{\mu ij,\nu k\ell}$ formally influence the argument
only through the symmetry properties that they
entail for $\Sigma_{\mu ij}$. Physically, however,
$C_{\mu ij,\nu k\ell}$ also set the energy
scales of glide and climb suppressions.

\section{Conclusions}
The dual formulation of elasticity theory that
exploits the concept of fractons allows for rather general insights
into the mobility of dissipative incommensurate crystals, such as twisted bilayer graphene.
In particular, we find that these systems have, in distinction to usual crystals, completely
immobile dislocations in the low defect density limit
that become sub-diffusive in the high-density limit.
Hence, while the electronic properties of graphene can, to a very good approximation, be treated using concepts of periodic crystals~\cite{LopesdosSantos2017,Trambly2010,Morell2010,Bistritzer2011}, the incommensurate nature of the material is much more visible in its mechanical properties.
In order to estimate whether the effects discussed here are quantitatively relevant,  we follow~\cite{Ochoa2019}. Here, $\kappa$ is shown to be of the order of $10 \text{eV}/\text{nm}^2  \sin\mleft(\Theta/2\mright)$. For twist angles $\Theta \sim 1^{\circ}$ this yields $\kappa \sim 0.1 \text{eV}/\text{nm}^2$. Hence, on the length scale of the moir\'e crystal $\sim 10 \text{nm}$, this stiffness is clearly relevant in a wide temperature regime.

Our results are not unique to twisted bilayer graphene. As we demonstrate in section~\ref{sec:III}, they can be generalized to generic planar quasicrystals and incommensurate crystals. The lack of mobility of defects stabilizes networks of defects, such as the soliton network discussed in Ref.~\cite{Ochoa2019}. Another implication is that these incommensurate systems should mechanically be rather brittle.

The power of the formalism that led to the results of this work
is drawn primarily from the deep intuition that
follows from the analogy to electromagnetism,
 including multipole expansions or the electromagnetism of dissipative media.
An interesting open question is how the properties of these defects,
and particularly the discompressions, impact the global mechanical
properties of TBG, as well as its local electronic spectrum, where stable defect configurations on the electronic spectrum are expected to lead to localized bound state formation. For isolated disclinations and discompressions one might even expect topologically protected bound states in the electronic spectrum~\cite{Peterson2021}.

\emph{Note added:} After this work was completed, we became aware of interesting related works concerned with a fracton description~\cite{Surowka2021} and topologically protected defect motion~\cite{Else2021} of quasicrystals.

\begin{acknowledgments}
We thank H. Ochoa for fruitful discussions. J.~G.\ was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) project SCHM 1031/12-1.
G.~P.\ and J.~S.\ were supported by the DFG - TRR 288 - 422213477 Elasto-Q-Mat (project A07).
R.M.F. was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under award no. DE-SC0020045.
R.M.F. also acknowledges a Mercator Fellowship (TRR 288 - 422213477) from the German Research Foundation (DFG).
\end{acknowledgments}


\section{Introduction}

Surveillance cameras, also known as closed-circuit television (CCTV) systems, have proliferated in the last several decades as the costs to record and store video have fallen dramatically.
As of 2016, there were an estimated 350 million surveillance cameras worldwide~\cite{ihs2016}.
The United States, with an estimated 50 million CCTV cameras installed, is believed to have the highest per capita number of surveillance cameras (15.3 CCTV cameras per 100 people) in the world~\cite{precise2019}.

Past work has found that surveillance cameras may play an important role in crime prevention and investigation, but there is also growing concern about the dangers cameras pose to privacy and equity.
Further, recent advances in facial recognition technology significantly amplify both the potential costs and the potential benefits of widespread surveillance,
as it is now possible to identify and track specific individuals across space and time.
While these technical advances promise to aid law enforcement efforts, they may also unjustly concentrate policing on more heavily monitored communities.
This surveillance may also hinder longstanding freedoms of speech and association, as it becomes easier to identify those participating in public gatherings, potentially dissuading dissent.

Despite the wide-ranging implications of surveillance cameras on public safety, police enforcement, and democratic governance, relatively little is known about the precise number and placement of cameras, hampering efforts to assess their impacts.
Past work to gauge the prevalence and spatial distribution of surveillance cameras has either examined aggregate production or shipping numbers, or relied on public disclosures in select jurisdictions---approaches that suffer from limitations of scale and scope.

To address these limitations,
\citet{turtiainen2020cctv} note that researchers could, in theory, map surveillance cameras by applying computer vision algorithms to street view data, which provide nearly complete visual coverage of many cities.
Building on that insight, here we describe and implement a scalable method for measuring the distribution of outdoor surveillance cameras across the United States, and, more generally, across the world.
Specifically, we couple computer vision algorithms with a verification pipeline by expert human annotators, together with statistical adjustment, to analyze a large-scale corpus of street view images.
In this manner, we leverage the proliferation of cameras and image data themselves to quantify the prevalence of surveillance technology.

To carry out this analysis, we use the public repository of images collected as part of Google's Street View service, launched in 2007. Since its inception, millions of 360-degree panoramas have been collected by cameras mounted on the roof rack of Google Street View cars, covering more than 10 million miles across 83 countries~\cite{raman2017}.
The rich archive of historical street view images provides opportunities to understand the evolution of the built environment, particularly the adoption of surveillance cameras on a global scale.
However, it is still extremely challenging---if not impossible---for humans to eyeball millions of images and spot cameras from the diverse street view context: a camera usually consists of 30--50 pixels out of over 400,000 pixels in a standard 640 $\times$ 640 street view image.

To scour this collection of images, we train and apply a computer vision algorithm to first filter street view data to those candidate images likely to contain a surveillance camera.
We specifically start with a random selection of 1.6 million street view images from 10 large U.S. cities and 6 other major cities, which contains approximately 6,000 positive model detections.
This curated set of candidate images is then examined by
human experts for verification.
To go from verified camera detections in our sample to city-wide estimates, we further
estimate both the recall of our model (which we find is 0.63), and
the proportion of the city covered by our sample.
This latter quantity is computed
based on the recorded camera position and angle, coupled with high-precision data on the road network and building footprints.

We find substantial variation in the density of visible surveillance cameras across the 16 cities we consider, ranging from 0.07 cameras per linear kilometer along the road network in Seattle to 0.95 cameras per kilometer in Seoul.
Examining the 10 U.S. cities in greater detail, we find that surveillance cameras are concentrated in
commercial, industrial, and mixed city zones,
and also in areas with higher shares of non-white residents.
This concentration of cameras in majority-minority neighborhoods persists even after adjusting for zone, pointing to the potential disparate impacts of surveillance technology on communities of color.

\section{Related work}
Our work connects to several interrelated strands of research in computer vision, urban computing, and privacy, which we briefly summarize below.

\subsection{Street View Understanding}
Visual scene understanding~\cite{hoiem2015guest} is one of the most fundamental and challenging goals in computer vision.
In part because of its potential to support self-driving vehicles, both the industrial and scientific communities have put considerable effort and investment into designing and creating labeled street view datasets for training and evaluating deep learning models, such as CamVid~\cite{brostow2008segmentation}, the KITTI Vision Benchmark Suite~\cite{geiger2013vision}, Cityscapes~\cite{cordts2016cityscapes}, and Mapillary Vistas~\cite{neuhold2017mapillary}.
Based on these datasets, several studies exploit the characteristics of urban-scene images and propose object segmentation~\cite{choi2020cars,liu2015layered} and change detection~\cite{alcantarilla2018street} algorithms for general street view understanding.

Related research has focused on detecting specific elements from street images, including greenery~\cite{li2015assessing}, buildings~\cite{kang2018building}, and city infrastructure such as utility poles~\cite{zhang2018using}.
Of particular relevance to our work,
\citet{neuhold2017mapillary} built an image segmentation model to identify---among other objects---CCTVs in street view data.
The publicly available \citeauthor{neuhold2017mapillary} Mapillary Vistas Dataset contains over 20,000 labeled images but fewer than 100 labeled cameras, leading to relatively poor performance on the specific task of detecting cameras.
More recently,
\citet{turtiainen2020cctv} developed a state-of-the-art object detection model tailored specifically to CCTVs,
based on nearly 10,000 images of cameras that they collected and labeled.
That dataset, however, has not been publicly released at the time of writing.
As a result, we constructed (and have released) our own labeled camera dataset and built a camera detection model using standard computer vision techniques.

\subsection{Urban Computing}
Urban computing aims to tackle major issues in cities---such as traffic control, public health, and economic development---by modeling and analyzing urban data.
A large body of research has shown that it is possible to infer socioeconomic information from satellite images~\cite{jean2016combining, Sheng_2020_CVPR_Workshops}, monitor human mobility~\cite{xu2018human}, and identify geo-tagged social network activities~\cite{schwartz2014social}.
Recent studies using street view images have dramatically increased the accuracy of processed data, as well as the geographic resolution analyzed.
By manually scoring street view images from 2,709 city blocks, \citet{hwang2014divergent} find that gentrification in Chicago from 2007 to 2009 was negatively associated with the concentration of minority groups.
\citet{mooney2016use} labeled the characteristics of 532 intersections in New York City, such as curb cuts and crosswalks, to assess environmental contributions to pedestrian injury.

As an alternative to relying on human experts to annotate street view images, modern computer vision algorithms have a much higher throughput with close to zero cost, enabling researchers to scale to multiple cities.
For example, \citet{gebru2017using} enumerated 22 million automobiles (8\% of all vehicles in the U.S.) in 50 million street view images to accurately
estimate local income, race, education, and voting patterns.
In our work, we draw on the merits of both approaches, combining high-recall computer vision algorithms with high-precision human verification in a unified estimation pipeline.
\input{figure/flowchart}
\subsection{Surveillance and Privacy}
 While past work has found that surveillance cameras play an important role in crime investigation~\cite{king2008citris} and deterrence~\cite{welsh2015effectiveness}, cameras also pose significant challenges to privacy.
Legal scholars have long considered the ramifications of cameras on First Amendment freedoms and the constitutional right to privacy~\cite{robb1980police}.
More recently, scholars have been concerned with the role of surveillance cameras in predictive policing \cite{joh2016discretion},
in enabling the adverse effects of facial recognition and computer vision~\cite{stanley2019robot,calo2010fake,buolamwini2018shades,nkonde2020black},
and with the threat of surveillance hacking~\cite{hermann2018hack,quintin2015license}.
These concerns have led to bans on facial-recognition technology by law enforcement in San Francisco, Boston, and Portland~\cite{banfacialrecognition},
as well as the drafting of federal legislation~\cite{frtbill}.

Despite these concerns, there has been limited success in identifying the number and geospatial distribution of cameras.
The Electronic Frontier Foundation (EFF) recently acquired the locations of cameras accessible by prosecutors in San Francisco~\cite{maass2019camera}.
Other private-market researchers have estimated the prevalence of installed cameras at a national level through unit shipments~\cite{jenkins2019surveil}.
However, neither of these approaches are able to estimate the prevalence and specific locations of public and private cameras at scale, hindering downstream analysis on the impacts of surveillance.

\section{Data and Methods}

For 16 major cities, we estimate the total number and spatial distribution
of surveillance cameras visible from the street.
We specifically consider the 10 cities with the highest urban density in the U.S., among those with at least 500,000 residents, and 6 other major cities in Asia and Europe.
Our statistical estimation procedure involves three key steps.
First, we compile a dataset of street view images both with and without cameras and label these images with segmentation masks.
We then train a camera segmentation model on this dataset,
and, importantly, estimate the
recall of our detection algorithm on a held-out validation dataset.
Second, we run our camera detection algorithm on a random sample of street view images.
All positive camera detections are then reviewed by human experts to remove false positives.
Finally, by combining the geometry of the camera angle,
the road network, and building footprints,
we calculate our sample's coverage of the road network.
These three steps are outlined in
Figure~\ref{fig:flowchart}.
In the following sections, we describe the data used in our analysis and more fully detail each step in our estimation pipeline.

\subsection{Data}
\input{table/city-stats}
\input{figure/example}
We analyze the 16 cities listed in Table~\ref{tab:citystats}.
For each city, we obtained the road network and building footprints from OpenStreetMap~\cite{OpenStreetMap,boeing2017osmnx}.
U.S.\ Census maps were used to restrict the geospatial data to the city's administrative borders.
All street view images used for model training and camera detection were accessed through the Google Static Streetview API.\footnote{\url{https://developers.google.com/maps/documentation/streetview}}
We further used San Francisco camera location data from the EFF~\cite{maass2019camera} to construct training and evaluation datasets for our model.

\subsection{Step 1: Model Training and Evaluation}

We start by creating training and evaluation datasets for our camera detection model.
For each of the 2,660 geo-tagged cameras in San Francisco identified by the EFF,
we pulled the closest street view images from 2012--2019 (if there is a scene available within 30 meters).
Manually labeling the resulting 13,240 images yielded 861 positive instances containing 977 cameras.
We note that many of the cameras listed in the EFF dataset appear to be indoors or otherwise are not visible from the street.
In Figure~\ref{fig:example}, we show several labeled examples.

We frame our camera detection problem as a binary image segmentation problem
to maximize learning from a limited number of samples.
We split the positive images by location into 70\%/15\%/15\% training/validation/test sets, making sure images from the same site always belong to the same split. We further include all camera instances from Mapillary Vista into our training dataset. By mixing with the negative images, we end up with 5,298 images for training, 1,040 for validation, and 1,040 for testing.

For ease, we use off-the-shelf methods to train our computer vision model (for state-of-the-art camera detection, see \citet{turtiainen2020cctv}).
In particular, our segmentation model follows the architecture of DeepLab V3+~\cite{chen2017rethinking,chen2018encoder} with an EfficientNet-b3~\cite{tan2019efficientnet} backbone.
We apply a random horizontal flip and randomly crop the original image (640 x 640 pixels) to 320 x 320 before feeding it into the model during training.
In the inference phase, we first crop the input image into four patches and then merge the output segmentation maps back to the original size.
The segmentation model's performance is shown in Table~\ref{tab:performance}.
To aggregate the pixel level prediction to the instance level, we first apply a morphological dilation with a 3x3 kernel to merge detected areas and filter false detections by size.  After validating several combinations of pixel-level probability thresholds and size filters,
we decided to use a probability threshold of 0.75 and a size threshold of 50 pixels, which yields precision and recall equal to 0.58 and 0.63, respectively (see Figure~\ref{fig:performance}).

In Figure~\ref{fig:failure}, we present several illustrative failures of our detection model. The model is occasionally confused by objects that share some of the visual features of cameras, such as building structures, parking meters, and street lamps. In some instances, our model also merges multiple cameras into one detection. These problems are mitigated by the human verification step, as described in the next section.
\input{table/performance}
\input{figure/performance}
\input{figure/failure}

\subsection{Step 2: Camera Detection and Verification}
\input{figure/spatial}
\input{figure/geometry}
\input{figure/verification}
For each city, we sampled street view images at $N=100,000$ points chosen uniformly at random from the road network.\footnote{For reference, there are more than 400,000 points covered with distinct street view panoramas in San Francisco.}
For approximately 3\% of the selected points, there was
no street view coverage within 10 meters, in which case we discarded and then re-sampled the location.
Figure~\ref{fig:spatial} shows the spatial distribution of the sampled points for three example cities: San Francisco, New York, and Chicago.

For each location, we then selected a 360-degree street view panorama.
For London, Paris, and the 10 American cities, we selected the oldest available image taken between 2015 and 2021;
for the remaining cities, we selected the oldest available image in the Google Maps corpus, which goes back to 2007.
We note that this sampling strategy is the result of a coding error; our intention was to select the \emph{newest} available image at each location.
Finally, for each location sampled, we randomly selected one out of the two 90-degree views with a midpoint perpendicular to the orientation of the road (see Figure~\ref{fig:geometry}). This approach provides the maximum view of the roadside.

We ran our camera detection model on the resulting set of 100,000 images for each of the 16 cities.
Annotators received the raw image and bounding boxes highlighting the predicted cameras, which were automatically generated from the model segmentation outputs. This process yielded 6,281 positive images with a total of 6,469 camera detections,
all of which were then verified by a human annotator.

Figure~\ref{fig:verification} illustrates the pipeline from the raw image to segmentation and bounding boxes to human verification.
In our subsequent analysis, we only consider these human-verified camera detections.

\input{figure/result_detection}
\subsection{Step 3: Road Network Coverage Estimation}
\input{figure/coverage}

The final step in our procedure is to estimate the fraction of the visible area covered by our randomly sampled images.
To estimate how much of the total length of a city's road network ($D$) has been covered by our sample ($D_{\text{covered}}$), we estimate the average length of road covered by one street view image ($\bar{d}$),
which we can then multiply by the number of images sampled ($N$).

We estimate $\bar{d}$ based on the geometry illustrated in Figure~\ref{fig:geometry}.
Each street view image comes with the exact latitude and longitude of where it was taken.
Given an image's point location $p$, we find the closest point $p'$ within the nearby buildings' footprint,
and denote the distance between $p'$ and $p$ with $\delta$.
As discussed above, we chose the street view's heading to be perpendicular to the road orientation,
and restricted it to a 90-degree view.
As a result, we estimate the length of road segment covered by the image to be $d  = 2\delta$.

We repeat this procedure for each sampled street view image.
We remove the relatively small number of images taken at locations more than 30 meters from a building---corresponding to 60 meters of street coverage---since at further distances, cameras become too small to be reliably detected by either humans or computer vision algorithms.
This results in images that cover about 25--30 meters of street.
For example, the mean road segment covered by an image $\bar{d}$ is 24, 29, and 28 meters in San Francisco, Chicago, and New York City, respectively, as shown in Figure~\ref{fig:coverage}.

Now, we estimate the proportion of a city's road network covered by our sample as $c = (N \bar{d}) / (2D)$,
where $N$ is the total number of samples for a city (within a given time period) and $D$ is the total length of the city's road network.
Note that the factor of 2 is to account for the fact that
our street view images only cover one of the two sides of a street at any sampled point.

Finally, putting all the above pieces together, we can estimate the number of cameras $K_{i}$ in city $i$:
\begin{equation}
 \hat{K}_{i} = \frac{n_{i}}{c_{i}r} \label{eq:estimation},
\end{equation}
where
$r$ is the recall of our model,
$n_{i}$ is the number of verified camera detections,
and $c_{i}$ is
the proportion of the road network of a city covered by our sample.

Similarly, we model variance by treating each sampled instance as a draw from a Bernoulli distribution with detection probability $p_{i} = n_{i} / N_{i}$. Assuming that recall and coverage are both exact, we can estimate the standard error of the number of cameras $K_{i}$ as:
\begin{equation}
 \hat{\text{se}}(\hat{K}_{i}) = \frac{\sqrt{N_{i} \cdot p_{i} \cdot (1-p_{i})}}{c_{i}r}.
 \label{eq:variance}
\end{equation}

\section{Results}
Applying the methods described above, we now estimate the total number and spatial distribution of cameras on the road network for all 16 cities.
In addition, for the U.S. cities, we estimate the prevalence of cameras across city zones, and examine the racial composition of the neighborhoods in which cameras are concentrated.

\subsection{Camera Prevalence}
\input{table/estimation}
\input{figure/result_density_cities}
Table \ref{tab:results} shows the number of identified cameras for each city---after human verification---along with point estimates and 95\% confidence intervals for camera density and for the total number of cameras, following Eqs.~\eqref{eq:estimation} and \eqref{eq:variance}.
The same density estimates are also depicted in descending order in Figure \ref{fig:result-density}.
We find that camera density varies widely between cities: For example, Boston and New York City, the U.S. cities with the highest camera density, have almost four times as many cameras per kilometer than Seattle and Los Angeles.\footnote{%
Our computer vision model was trained on San Francisco data, and so it is possible that the camera identification rate in San Francisco is inflated due to over-fitting.
We note, however, that while model precision varies between cities, Philadelphia and Boston both have higher precision than San Francisco, which suggests our model does indeed transfer well across contexts.
}
We note that our estimates exclude indoor cameras,
as well as outdoor cameras not captured by street view images.
Perhaps due to these limitations, our estimate of 10,100 cameras in New York City is lower than the 18,000 cameras that the NYPD reportedly has access to~\cite{nypd_camera_count}.

\subsection{Camera Placement}
\input{figure/result_zone}
\input{figure/result_race}

The detection maps in Figure~\ref{fig:detections} show that cameras are not distributed uniformly across a city. Despite sampling uniformly over the road network, we find densely covered regions in each city, representing neighborhoods with a high concentration of cameras.
We examine these patterns in more detail for the 10 U.S. cities we consider, analyzing the rate of (verified) camera identifications per street image across zoning designations and neighborhood racial composition.

Figure~\ref{fig:result-zone} shows the camera identification rate for different zoning designations aggregated over the 10 U.S. cities we analyze.
We find that images from mixed, industrial, and commercial zones are more likely to contain an identified camera than images from public (such as parks and other public facilities) and residential areas.
For example, the identification rate in mixed zones (2.1\%) is more than three times the rate in residential zones (0.6\%).
This pattern holds for the majority of our chosen cities.

To compute the camera identification rate, we assigned each sampled point to the zoning designation of the closest parcel of land. To do so, we collected land use and zoning designation data for all 10 cities, and then standardized the zoning code into one of the following five categories: mixed, industrial, commercial, public, and residential. Zones with codes that represent planned development or did not clearly fit into the aforementioned categories are labeled as unknown and
omitted in the following analysis. We find that 60\% of sampled points are classified as residential, and unknown codes comprise less than 3\% of sampled points.

We next examine the relationship between camera identification rate and the share of residents in the surrounding area that identify as belonging to a minority racial or ethnic group,
aggregated over our 10 U.S. cities.
To compute this relationship, we assigned each sampled image to the minority proportion of the census block group in which it is located, as estimated by the 2018 American Community Survey.
For purposes of this analysis, we define ``minority'' as comprising those individuals who identify either as Hispanic (regardless of their race) or who do not identify as white.

Figure \ref{fig:result-race} shows the results.
The blue line is a regression (with both linear and quadratic terms) fit to the data, and indicates that an increase in the share of minority residents in a neighborhood is associated with an increase in camera identification rate.
For example, the identification rate in census blocks with a 50\% minority share (0.38\%) is twice as high as in those blocks with a 10\% minority share (0.2\%).
We see qualitatively similar results with higher-order polynomial fits.

\input{table/regression}
The observed concentration of cameras in majority-minority neighborhoods persists even after adjusting for zone category.
Specifically, Table~\ref{tab:regression}
shows the results of a linear probability model
that predicts camera detections as a function of city, zone, and racial composition---where we again use a quadratic term to account for the curvature seen in Figure~\ref{fig:result-race}.
The fitted model confirms that camera identifications increase with the minority share of residents, plateauing at approximately 60\% share and then remaining relatively flat, as in Figure~\ref{fig:result-race}.
It is unclear what is driving the apparent concentration of cameras in high minority neighborhoods.
However, regardless of the underlying mechanism, these results point to the potential impacts that video surveillance can have on communities of color.

\section{Discussion and Conclusion}
By applying computer vision, human verification, and statistical analysis to large-scale, geo-tagged image data,
we have---for the first time---estimated the number and spatial distribution of outdoor surveillance cameras
in 16 major cities around the world.
Further, the approach we have developed has the potential to scale to even more cities across the country and the world, providing a new perspective on the state of video surveillance.

In the 16 cities we analyzed,
we found considerable variation in the estimated number of surveillance cameras.
Among U.S. cities, our analysis also shows that cameras are more likely to be found in industrial, commercial, and mixed zones as compared to residential areas.
Finally, even after adjusting for zone category,
we find a greater concentration of cameras in
majority-minority neighborhoods, highlighting the need to carefully consider the potential disparate impacts of surveillance technology on communities of color.

While our computational approach is able to provide a novel quantitative perspective into the state of surveillance, it
is still subject to several important limitations, which we outline below.
First, our method relies on being able to see cameras from the street, and, more specifically, from street view images.
Indoor cameras, as well as outdoor cameras obscured from view,
are not counted by our estimation pipeline.
Further, due to the limited resolution of street view images, small cameras---such as increasingly popular doorbell cameras---are difficult to detect by either humans or algorithms.
Higher resolution and higher coverage image data could mitigate these issues in the future.
However, our current results likely underestimate the density of cameras in a city.

Second, our human annotators may not perfectly label cameras in the candidate images selected by the model, skewing our final estimates.
For example, it is possible that they rule out an actual camera (leading to an undercount) or, conversely, that they report a camera that is not in fact there (leading to an overcount).
To minimize these errors, every candidate image is independently labeled by three human annotators,
but at least some errors likely remain.

Third, errors in the estimated recall of our computer vision model---and, similarly, errors in the estimated coverage of our images---can bias our final estimates.
Estimating city-specific model recall is particularly challenging, as it requires city-specific labeled datasets.
In our analysis, we thus estimated recall for a single city, San Francisco, in which the locations of some surveillance cameras had already been compiled, which we then apply to other jurisdictions.
Further, our variance estimates treat the recall and coverage as known quantities.
Accounting for errors in their measurement would increase the variance of our final estimates.

Finally, our method does not provide any information about the cameras other than what can be inferred from their appearance.
For example, we cannot determine whether identified cameras are decoys,
are malfunctioning, or otherwise are not in use.
We likewise cannot always tell who owns the cameras (e.g., a government agency or a private citizen), who has access to the video,
and whether the camera footage is stored.
All of these factors are critical in assessing the downstream consequences of video surveillance.
Although difficult, future work may be able to answer some of these questions by conducting
a more intensive audit of a sample of the identified cameras.

Despite these limitations,
we believe our approach and results constitute an important step toward understanding the use of surveillance technology across the world.
More broadly, our
general statistical estimation pipeline can be extended and applied to characterize the prevalence and spatial distribution of a variety of other city elements detectable from street images.
Looking forward, we hope this work spurs further theoretical and empirical research at the intersection of computer vision, urban computing, and public policy.

\section*{Publication Note}
This version of the paper is updated from our original publication in two important respects.
First, we now credit
\citet{turtiainen2020cctv} both for creating a state-of-the-art camera detection model and for suggesting that computer vision could, in theory, be applied to street view data to map surveillance cameras.
We were aware of their work when initially conducting our research, but we unfortunately failed to include a citation to their paper.
We thank \citeauthor{turtiainen2020cctv} for bringing this to our attention and we apologize for the omission.
Second, we discovered a coding error in our image sampling strategy that corrupted our analysis of camera density over time. We have now removed the results of that analysis.

\balance
\bibliographystyle{ACM-Reference-Format}