# Finding hidden-feature depending laws inside a data set and classifying it
using Neural Network.

Thilo Moshagen  Nihal Acharya Adde  Ajay Navilarekal Rajgopal

###### Abstract

The logcosh loss function for neural networks has been developed to combine
the advantage of the absolute error loss function of not overweighting
outliers with the advantage of the mean square error of continuous derivative
near the mean, which makes the last phase of learning easier. It is clear, and
one experiences it soon, that in the case of clustered data, an artificial
neural network with logcosh loss learns the bigger cluster rather than the
mean of the two. Even more so, the ANN, when used for regression of a set-
valued function, will learn a value close to one of the choices, in other
words, one branch of the set-valued function, while a mean-square-error NN
will learn the value in between. This work suggests a method that uses
artificial neural networks with logcosh loss to find the branches of set-
valued mappings in parameter-outcome sample sets and classifies the samples
according to those branches.

Keywords— Neural Networks, Clustering, Classification, Model selection, Loss
function, Objective Function, ANOVA, Hypothesis testing

##  1 Introduction

Given a set of data tuple, _Clustering_ algorithms [AR] decide which elements
of the set belong together, i.e. form a subset in the sense that they have
closer mutual distance among each other. Further, there are a lot of well
established and also new methods to deal with the question of whether two or
more sets of samples belong to the same population or not. Mainly, this
considers the field of statistical hypothesis testing [SOA99] which is a
testable hypothesis based on observed data modelled as the realised values
taken by a collection of random variables. Also, in a different setting, a
data model can be defined as a set of mathematical laws that might be valid
inside a data set and describes how the data elements relate to one another.
Given that there exist some measurements and parameters that caused these
measurements, the model selection tells which model is most likely valid for
the observed measurements to happen. Variants of ANOVA (Analysis of variance)
combine the two and is used to analyse the differences among group means in a
sample [KS14].

We consider a method that answers the question whether a set of vector-valued
samples, where some components can be seen as cause and at least one as an
effect, obeys some possibly unknown rule, or if it rather splits into groups
that fulfil different rules. In other words, assuming that any input data
point consists of components
$\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$
that are presumably independent and component(s)
$\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}$
that depend on them by some generally unknown rules, the suggested methods
finds the rules in the shape of an artificial neural networks’ weights _and_
clusters the data into groups obeying each of the found rules.

The three key features of the suggested method are, first, the use of neural
networks to extract one of the unknown rules that are valid in parts of the
data. This extraction is done by supervised learning, which is a regression in
mathematical terms. The second key feature is that supervised learning is done
with a loss function that is approximately linear in the distance to zero and
thus puts less weight on far-off data than the square error loss function. For
example, the $L_{1}$-norm fulfils this. But here, the logcosh loss function
was used as it facilitates learning, while still having the desired property.
When used for regression, such a loss function leads to learning a function
that approximates well the strongest cluster of output data, while hardly
taking into account clusters with fewer members. Data lying away from the
found regression graph thus is probably obeying another law; Distance to the
regression function found by the artificial neural network is then used as a
classification criterion. This is the third key feature. The data that is
approximated well by the found regression function is considered to be
governed by that function. With the badly approximated data, a new network is
taught, and all points where its forecast matches are considered to be
governed by that second regression. This procedure is continued until no
relevant data remains unclassified. This is how in brief the suggested method
works.

##  2 Problem Setting

###  2.1 Mathematical Description

Let

$\left\\{(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}},\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}})\in\mathbb{R}^{d}\times\mathbb{R}^{d}\right\\}$
(1)

be a set of data points, where it can be assumed that
$\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}$
depends on
$\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$.
To simplify the setting and also due to the fact that artificial neural
networks do not encourage vector valued output, we restrict ourselves to

$\left\\{(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}},y)\in\Omega\subset\mathbb{R}^{d}\times\mathbb{R}\right\\}$
(2)

where $\Omega$ is the domain in which observations are defined.

The presence of clusters
$\left\\{(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}},y)\right\\}_{i,i=1,...,M}$,
where inside each cluster, the
$(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}},y)$-tuples
obey a different law is now mathematically described as follows:

Each clusters’ independent variable points are subsumed in the set
$\hat{X}_{i}\subset X$, $X$ being all
$\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$
in $\Omega$, where the law valid inside it is (first defined on the samples
only, with the hat denoting this):

$\displaystyle\hat{\phi}_{i}:\hat{X}_{i}$
$\displaystyle\longrightarrow\mathbb{R}$ (3)
$\displaystyle\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$
$\displaystyle\mapsto\hat{\phi}_{i}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})=y(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})$
(4)

where $\hat{\Phi}_{i}$ is a single-valued function, mapping each point in
$\hat{X}_{i}$ to a unique value in the range. The existence of multiple
$\hat{\Phi}_{i}$ is due to hidden features, for which nearby
$\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$
can have very distant $y$. Each $\hat{\Phi}_{i}$ induces a continuous function
$\Phi_{i}$ in some super-set $X_{i}$ of $\hat{X}_{i}$ by regression, the
continuous regression counterpart of the measurements. Those then have
reasonably bounded derivatives - which a mapping $\Phi$ that maps all $x$
would not have. There may exist a certain subset of $X$ which gives the same
output for all $\Phi_{i}$, while in the $X_{i}$ the $\Phi_{i}$ give different
values. Thus, $\Phi_{i}$ may be seen as defined only on $X_{i}$, or
alternatively on $X$, in which case the $\Phi_{i}$ coincide in parts of
$\Omega$. This can be seen as a multi-valued or set-valued function

$\displaystyle\Phi(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})=\left\\{\begin{array}[]{cc}\Phi_{1}\\\
\vdots\\\ \Phi_{M}\end{array}\right\\}.$ (8)

The set-valuedness in this nomenclature is expressed by this vector-
valuedness. It captures the property that the data input-output pairs indeed
belong to different situations or populations. The task to solve in this
nomenclature is: Given the set $X$, find the rules $\Phi_{i}$ and the subsets
$X_{i}$ where they are valid.

###  2.2 Outline of Strategy

One seeks to learn each $X_{i}$’s rule $\Phi_{i}$ by regression, which for
general $\Phi_{i}$ is done best by artificial neural network, using the
logcosh loss function: It weights the outliers less, similar to the MAE loss
while it exhibits good performance during gradient descent as MSE. The network
trained with logcosh loss will thus learn the biggest cluster $\Phi_{1}$
efficiently because it weights smaller clusters away from the biggest one only
linearly with distance, unlike the squared error losses, and thus classifies
the data as belonging to the biggest cluster or not. In our research, we train
the network with logcosh loss function in an aim to classify the clustered
data. This approach is demonstrated using a simple 1-dimensional and
2-dimensional problem.

##  3  Artificial Neural Network Regression Quality as a Classification
Criterion

Supervised learning of an Artificial Neural Network [GBC16] has the task of
learning a function that maps an input to an output based on example input-
output pairs. It is where the set of input variables
$\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$
and the output variables
$\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}$
are available and one has to use an algorithm to learn the mapping function
from the input to the output. The goal is to approximate the mapping function
so well that the new unseen input data
$\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$
can be used to predict the output variables
$\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}$
for that data. An ANN is based on a collection of connected nodes called
neurons which loosely represents the neurons in a biological brain. Each
connection transmits signals from one neuron to the other. The signal at a
connection is a real number, and the output of each neuron is computed by some
non-linear function of the sum of its inputs. These connections are called
edges. Neurons and edges typically have a weight that adjusts as learning
proceeds. Through backpropagation, the network tries to find optimal weights
and biases to represent the model. In other words, the artificial neural
network can be represented as an optimization problem which ultimately is
equivalent to minimising the loss function of the data. Therefore the choice
of the loss function becomes vital for modelling an efficient network. Our
task is to find an appropriate model that fits the regression model by one of
the rules and classifies the clustered data by it. Therefore, in the following
section, we discuss the different available loss functions and choose an
appropriate loss function in an aim to classify the clustered data.

###  3.1 Loss Functions Properties

One key feature of the suggested method is the choice of the loss function. We
will point out in the following that for a regression problem, the minimizer
of loss functions that rise linearly with the distance lies inside a cluster,
while for the quadratic loss functions, it lies between clusters. The choice
of loss function depends on a number of factors including the presence of
outliers, choice of the machine learning algorithm, time efficiency of
gradient descent, ease of finding the derivatives and confidence of
predictions. [NZL18] investigated some representative loss functions and
analysed the latent properties of them. The main goal of the investigation was
to find the reason why bilateral loss functions are more suitable for
regression task, while unilateral loss functions are more suitable for
classification task. This section covers in detail the different loss
functions which can be used for our regression problem as discussed by
[Gro18].

####  3.1.1 Mean Square Error (MSE) or L2 loss

This function originates from the theory of regression, least-squares method.
Mean Square Error (MSE) is the most commonly used regression loss function.
MSE is the sum of squared distances between our target variable $y$ and
predicted values $y_{p}$.

$\mathbf{M.S.E.}=\frac{\sum_{i=1}^{n}(y^{i}-y_{p}^{i})^{2}}{n}$ (9)

It is well known that, here, few distant points outweigh the closer points.
The MSE loss establishes that our trained model takes outliers seriously as
the contribution to loss by an outlier in input is magnified by squaring and
so learning results are biased in favor of the outliers. This can be an
advantage - predictions in zones with outliers do not produce huge errors to
the outliers since the MSE took them into account. MSE is thus good to use if
the target data conditioned on the input is normally distributed around a mean
value and in the absence of outliers. It has a continuous derivative and
therefore the minimisation with gradient methods works well. The described
property is a disadvantage for our setting, as one cluster consist of outliers
seen from the other clusters’ perspective, thus the MSE minimiser would be
right in between clusters.

Figure 1 shows the plots of mean square error loss vs. predictions, where the
target value is 0, and the predicted values range between -100 to 100. The
loss (Y-axis) reaches its minimum value at the prediction (X-axis) = 0. The
range of the loss is 0 to $\infty$.

Figure 1: Plot of Mean Square Error (MSE) Loss

####  3.1.2 Mean Absolute Error (MAE) or L1 loss

Mean Absolute Error (MAE) is just the mean of absolute errors between the
actual value $y$ and the value predicted $y_{p}$. So it measures the average
magnitude of errors in a set of predictions, without considering their
directions.

$\mathbf{M.A.E.}=\frac{\sum_{i=1}^{n}|(y^{i}-y_{p}^{i})|}{n}$ (10)

As one can see, for this loss function, both the big and small distances
contribute the same. The advantage of MAE covers the disadvantage of MSE. As
we consider the absolute value, the errors will be weighted on the same linear
scale. Therefore, unlike the previous case, MAE doesn’t put too much weight on
the outliers. However, it does not have a continuous derivative and thus
frequently oscillates around a minimum during gradient descent. The MSE does a
better job there as it has a continuous derivative and provides a stable
solution. Figure 2 shows the plot of mean absolute error loss with respect to
the prediction while the target value is 0, similar to the previous case.

Figure 2: Plot of Mean Absolute Error (MAE) Loss

####  3.1.3 Huber loss

Huber loss is just the absolute error but transforms to squared error for
small values of error. It is an attempt to overcome MAE’s disadvantage of
nonsmooth derivative. Huber loss is less sensitive to outliers in data than
the squared error loss. It is also differentiable at 0. It is basically
absolute error, which becomes quadratic when the error is small. How small
that error has to be to make it quadratic depends on a hyperparameter
$\delta$, which can be tuned. Huber loss approaches MSE when
$\delta\rightarrow 0$ and MAE when $\delta\rightarrow\infty$ (large numbers).
It is defined as

$L_{\delta}(y,y_{p})=\left\\{\begin{array}[]{ll}\frac{1}{2}(y-y_{p})^{2}&\mbox{if
}|y-y_{p}|\leq\delta\\\ \delta|y-y_{p}|-\frac{1}{2}\delta^{2}&\mbox{otherwise
}\end{array}\right\\}$ (11)

The choice of $\delta$ becomes increasingly important depending on what one
considers as an outlier. Residuals larger than delta are minimized with L1
while residuals smaller than delta are minimized with L2. Hubber loss combines
the advantages of both the loss functions. It can be really helpful in some
cases, as it curves around the minima which decreases the gradient. However,
the problem with Huber loss is that we might need to train hyperparameter
delta which is an iterative process. Figure 3 shows the plot of Huber loss vs.
predictions for different values of delta $\delta$.

Figure 3: Plot of Huber Loss

####  3.1.4 Log-Cosh loss

Log-cosh is the logarithm of the hyperbolic cosine of the prediction error.
Given the actual value $y$ and the predicted value $y_{p}$, the log-cosh is
defined as

$L(y,y_{p})=\sum_{i=1}^{n}|\operatorname{log}(\cosh((y^{i}-y_{p}^{i})))|$ (12)

$\operatorname{log}(\cosh(x))$ is approximately equal to $\frac{x^{2}}{2}$ for
small values of x and to $|x|-\operatorname{log}(2)$ for larger values. It is
twice differentiable everywhere unlike Huber loss. Therefore, the log-cosh
loss function is similar to mean absolute error with respect to its moderate
weighting of outliers, while it behaves stable during gradient descent search.
Figure 4 shows the plots of logcosh loss vs predictions, where the target
value is 0, and the predicted values range between -10 to 10.

Figure 4: Plot of Log-cosh Loss

Therefor, in our research log-cosh loss function was used and indeed showed
good results in classifying the data based on the hidden features.

Figure 5 compares the 3 different losses functions.

Figure 5: Plot of different Losses : MSE, MAE and Log-cosh

When clustered data is present, an artificial neural network with logcosh loss
function learns the bigger cluster rather than the mean of the two and hence
can be used to classify the clustered data. In the case of MSE, due to the
squaring of the error function, few faraway points are weighted more than the
nearby points. When learning clustered data, the network with MSE loss
function gets affected by these outlying clusters and tries to find the minima
between them and thereby fails to learn the bigger cluster. For linearly
growing loss functions like logcosh and MAE, just the sum of distances counts
and few far-away points do not count more than several nearby points and
therefore, a regression value near or through the heavier cluster is learnt.
Though the MAE loss function has this property of the bigger cluster, it is
non-smooth and has a non-continuous derivative resulting in oscillating
behaviour. As mentioned above, since the logcosh loss function is a
combination of MAE for larger values and MSE for the smaller values, it
successfully learns the bigger cluster and gives a stable solution. These
features of logcosh loss function are exploited in our research.

###  3.2 One-Dimensional Test Case

####  3.2.1 Test Problem

We now consider a simple 1D example based on the concept discussed in the
section 2. Two simple single-valued polynomial functions were selected and
combined in different fractions to achieve a multi-valued data set. This
section discusses the problem setting of the 1-dimensional case and thereafter
the network behaviour based on the chosen data set.

To create a multi valued Data set, 2 simple functions were selected as below.

$\Phi_{1}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})={((x-4)(x+4))}^{2},\hskip
14.22636ptx\in[-6,6]$ (13)
$\displaystyle\Phi_{2}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})=\left\\{\begin{array}[]{cc}{((x-4)(x+4))}^{2}&\hskip
14.22636ptx\in[-6,-4)\\\ 0&\hskip 14.22636ptx\in[-4,4]\\\
{((x-4)(x+4))}^{2}&\hskip 14.22636ptx\in(4,6]\\\ \end{array}\right\\}.$ (17)

where $\Phi_{1}$ and $\Phi_{2}$ are two single-valued functions which are
defined within the interval $[-6,6]$.

####  3.2.2 Training Strategy

The data set was split such that $80\%$ of the data were used for training and
the rest $20\%$ were used as test data. Initially, both the functions were
trained individually with a basic regression neural network and then tested on
the test data to validate the network.

(a) $60\%$ or more of function $\phi_{1}$

(b) $60\%$ or more of function $\phi_{2}$

Figure 6: Plot of Test and Predicted Data for the functions $\phi_{1}$ and
$\phi_{2}$ for 1-dimensional test case.

As seen in figure 6, it is clear that the neural network was able to
approximate the given functions by reducing the loss function to the minimum.

As discussed in Section 2, to set up a multi-valued data set we combine
fraction of both the sets $\Phi_{1}$ and $\Phi_{2}$ respectively, to form a
new data set $\Phi$ as per our requirement. The two data sets were combined in
different fractions, trained using our neural network and then tested on the
test data which is $20\%$ of the new combined data. The noise was added to the
data set to replicate the real-world data. The network was trained using
logcosh loss function to examine the network behaviour. To compare the
functionality of different loss functions, the network was also trained with
MSE and MAE loss function using a similar setting.

The combined function can be written as follows :

$\displaystyle\Phi(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})=\left\\{\begin{array}[]{cc}\Phi_{1}\\\
\Phi_{2}\end{array}\right\\}.$ (20)

where
$\Phi(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})$
is a combination of both multi-valued (the multi-valued region where each
$\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$
has two possible outputs $y$ as shown in the figure 7 and single valued
function (defined on $[-6,-4)\bigcup(4,6]$).

Figure 7: Plot of data set with noise for 1-dimensional case

####  3.2.3 Network behaviour

In this section, the behaviour of our network based on the chosen network
architecture is discussed. As discussed earlier the network was trained with a
different fraction of the two chosen functions and then tested on the test
data. The network was completely trained using log-cosh loss function. It can
be seen that, when using the log-cosh loss, the network predicted one of the
two chosen function with high accuracy and not the mean of the two functions.
The network predicted the function $\Phi_{1}$, when $60\%$ or more of the
function $\Phi_{1}$-generated data was chosen in the combined data set and it
predicted the function $\Phi_{2}$ otherwise, as shown in the figure 8.

(a) $60\%$ or more of function $\Phi_{1}$

(b) $60\%$ or more of function $\Phi_{2}$

Figure 8: Plot of Test and Predicted Data for the mixed data set in the
1-dimensional case. The plots illustrate the behaviour of the network when
tested on the test data. The network accurately predicts one of the 2
functions depending on the fraction of functions considered when trained using
log-cosh loss function.

As could be expected (see Section 3.1.4), the logcosh loss function learned
the bigger cluster of data, unlike the mean square error loss which learned
the mean of the two functions or the absolute error which would oscillate
between the two chosen functions as shown in Figure 9. As mentioned in section
3.1.4, the MSE loss function gets affected by the minor cluster due to
squaring and thus finds the weighted mean between the two depending on
composition of the cluster. Unlike MSE, MAE functions similar to that of
logcosh and tries to find one of the two clusters. However, since it is non
smooth and has non-continuous derivative, the prediction oscillates between
the clusters when the composition of the clusters are nearly equal.

(a) Network behaviour when the model was trained with MSE loss function. The
network predicts the weighted mean of the 2 functions depending on the
fractional composition

(b) Network behaviour when the model was trained with MAE loss function.
Though the network predicts one of the 2 clusters often, the prediction
oscillates between the two when the fractions of the clusters are nearly
equal.

Figure 9: Plot of Test and Predicted Data for the mixed data set in the
1-dimensional case when the network is trained using (a) MSE and (b) MAE loss
function.

###  3.3 Two-Dimensional Test Case

####  3.3.1 Test Problem

We now choose a 2-dimensional case based on the concept discussed in section
2. Similar to the 1D case, two 2 dimensional single-valued functions were
combined in different fractions to form the multi-valued data set $\Phi$,
learned by the neural network and finally, the behaviour of our network based
on these data sets was analysed.

The two functions

$f_{1}(x,y)=xy(2x+2y)$ (21) $f_{2}(x,y)=xy(x^{2}+y^{2})$ (22)

were used as arguments to the sigmoid function. The main reason to use the
sigmoid function was to keep the range between (0,1).

$\Phi_{1}(x,y)=\text{sigmoid}(f_{1}(x,y))=\frac{1}{1+e^{-f_{1}(x,y)}}$ (23)
$\Phi_{2}(x,y)=\text{sigmoid}(f_{2}(x,y))=\frac{1}{1+e^{-f_{2}(x,y)}}$ (24)

To set up a multi-valued data set we combined both the data sets $\Phi_{1}$
and $\Phi_{2}$ of the above functions in different fractions to form a
combined data set $\Phi$ as per our requirement. Noise was added to the data
set to replicate the real-world scenario.

####  3.3.2 Training Strategy

The neural network was trained with this data set and then predicted on the
test data which is $20\%$ of the total combined data.

Figure 10: Plot of data set without noise for a 2-dimensional data set.

Figure 10 shows the plot of the combined data set without noise, where red and
orange represent the function $\Phi_{1}$ and function $\Phi_{2}$ respectively.
As discussed earlier, in this case, for given nearby $(x,y)_{1}$ and
$(x,y)_{2}$, we have two distant values $z_{1}$ and $z_{2}$ despite being very
close to each other. The network with logcosh loss function is trained with
different fractions of the sets $\Phi_{1}$ and $\Phi_{2}$ joined into $\Phi$
in an aim to classify the two.

####  3.3.3  Network behaviour

A very noisy data set was used to train the network – the 2 populations cannot
be easily distinguished by visualization. After training the network with a
combination $\Phi$ of different fractions of $\Phi_{1}$ and $\Phi_{2}$,
similar to the 1D case, a clear rule was visible when the logcosh loss
function was used.

(a) $60\%$ or more of function $\Phi_{1}$ values in data

(b) $60\%$ or more of function $\Phi_{2}$ values in data

Figure 11: Plot of Test and Predicted Data for the mixed data set in the
2-dimensional case. The network accurately predicts one of the two functions
depending on the fractional composition of the data set.

The network predicted the function $\Phi_{1}$ when $60\%$ or more of $\Phi$
consisted values of $\Phi_{1}$ and $\Phi_{2}$ when $60\%$ or more of $\Phi$
consisted values of $\Phi_{2}$ as shown in the figure 11. In the plots, the
red scatter points represent the function $\Phi_{1}$ with noise and the red
surface plot represents the values $\Phi_{1}$ of function ${\Phi}$ without
noise. Similarly, for the function $\Phi_{2}$, orange scatter points and
orange surface plot represents the function with and without noise
respectively. Finally, the blue scatter points represent the predicted value.
The functions were plotted without noise for better visualisation. From figure
11, it is clear that the network learnt one of the 2 functions accurately
without being influenced by noise. It can be therefore confirmed that the
neural network predicts the bigger of the two clusters when logcosh loss
function is used.

##  4 Conclusion

Based on the network behaviour, we claim that a network with logcosh loss
function can be used to classify the data when clusters of data exist. It can
be concluded that in case of clustered data, an artificial neural network with
logcosh learns the bigger cluster rather than the mean of the two. Even more
so, the ANN when used for regression of a set-valued function, will learn a
value close to one of the choices, in other words, one branch of the set-
valued function, while a mean-square-error NN will learn the value in between.
Based on the above result we have a neural network that not only helps in
classifying the data based on the invisible features but also predicts the
majority cluster with high accuracy. In the real world scenario, the
unavailability of enough parameters to build the regression model is always a
major problem and therefore it becomes increasingly difficult to represent the
model based on the available limited data. Using this theory, we can classify
the clusters of data based on an invisible feature which is not available to
us beforehand. It can be also used to validate if there are enough features to
represent the model. In other words, we can confirm if a feature is essential
to represent the model.

## References

  * [AR] Charu C. Aggarwal and Chandan K. Reddy “Data Clustering” In _O’Reilly Online Learning_ ChapmanHall/CRC URL: https://www.oreilly.com/library/view/data-clustering/9781466558229/
  * [GBC16] Ian Goodfellow, Yoshua Bengio and Aaron Courville “Deep Learning” http://www.deeplearningbook.org MIT Press, 2016
  * [Gro18] Prince Grover “5 Regression Loss Functions All Machine Learners Should Know”, 2018 URL: https://heartbeat.fritz.ai/5-regression-loss-functions-all-machine-learners-should-know-4fb140e9d4b0
  * [KS14] Jörg Kaufmann and AG Schering “Analysis of Variance ANOVA” In _Wiley StatsRef: Statistics Reference Online_ American Cancer Society, 2014 DOI: https://doi.org/10.1002/9781118445112.stat06938
  * [NZL18] Feiping Nie, Hu Zhanxuan and Xuelong Li “An investigation for loss functions widely used in machine learning” In _Communications in Information and Systems_ 18, 2018, pp. 37–52 DOI: 10.4310/CIS.2018.v18.n1.a2
  * [SOA99] Alan Stuart, J. Ord and Steven Arnold “Kendall’s advanced theory of statistics. Vol.2A: Classical inference and the linear model”, 1999


# Effective Field Theory of Dipolar Braiding Statistics in Two Dimensions

Yun-Tak Oh  Department of Physics, Korea Advanced Institute of Science and
Technology, Daejeon 34141, Republic of Korea  Division of Display and
Semiconductor Physics, Korea University, Sejong 30019, Korea  Jintae Kim
Department of Physics, Sungkyunkwan University, Suwon 16419, Korea  Jung Hoon
Han<EMAIL_ADDRESS>Department of Physics, Sungkyunkwan University, Suwon
16419, Korea

###### Abstract

A rank-2 toric code (R2TC) Hamiltonian in two dimensions can be constructed as
a Higgsed descendant of rank-2 U(1) lattice gauge theory. As noted by the
authors recently, [Y.-T. Oh, J. Kim, E.-G. Moon, and J. H. Han, Phys. Rev. B
105, 045128] the quasiparticles in that model shows unusual braiding
statistics that depends on the initial locations of the particles which
participate in the braiding. We show that this new kind of statistical phase
captures the total dipole moment of quasiparticles encompassed in the
braiding, in contrast to the conventional anyonic braiding seeing the total
charge. An Aharonov-Bohm interpretation of such dipolar braiding statistics is
made in terms of emergent, rank-1 vector potentials that are built out of the
underlying rank-2 gauge fields. Pertinent field theories of the quasiparticle
dynamics in the R2TC are developed, and the accompanying conservation laws
derived. A dipolar BF theory of the rank-2 gauge fields is constructed and
shown to correctly capture the dipolar braiding statistics, in contrast to the
conventional BF theory capturing the monopolar braiding statistics of anyons
in the rank-1 toric code.

##  I Introduction

The Chern-Simons (CS) field theory Deser  _et al._ (1982) made its entry into
condensed matter physics as an effective field theory of the fractional
quantum Hall system Zhang  _et al._ (1989); Zhang (1992) and encodes the
mutual statistics of anyonic quasiparticles Wilczek (1982a, b). A variant of
the CS theory known as the BF theory captures the mutual statistics among
different species of anyons, such as those of electric and magnetic
quasiparticles in toric code (TC) or of vortex and quasiparticle in
superfluids Hansson  _et al._ (2004). In both CS and BF theories, the phase
factor picked up in the process of braiding one anyon around other anyons is
proportional to the number of anyon charges enclosed and can be understood
within the flux attachment picture in which one anyon sees the effective
magnetic flux tied to other anyons and picks up the corresponding Aharonov-
Bohm (AB) phase.

Recently, sub-dimensional topological order embodied in various exactly
solvable spin models has emerged at the frontier of quantum matter research
Haah (2011); Vijay  _et al._ (2016); Pretko (2017); Slagle and Kim (2018);
Pretko and Radzihovsky (2018); Nandkishore and Hermele (2019); Pretko  _et
al._ (2020); Kim and Han (2021). These sub-dimensional topological orders are
characterized by excitations such as fractons or lineons with constrained
mobilities. Various exotic field theories attempting to capture aspects of
fracton physics are under vigorous pursuit Slagle and Kim (2017); Pretko
(2018); Shirley  _et al._ (2018); You  _et al._ (2020); Dubinkin  _et al._
(2020); Seiberg and Shao (2020, 2021a, 2021b); You  _et al._ (2021). It was
also realized that applying the Higgsing procedure to a lattice gauge theory
(LGT) yields an exactly solvable spin model while lowering U(1) gauge symmetry
to $\mathbb{Z}_{N}$ Ma  _et al._ (2018); Bulmash and Barkeshli (2018); Oh  _et
al._ (2022); Kim  _et al._ (2022). A well-known X-cube model of fractons Vijay
_et al._ (2016) is obtained by Higgsing the rank-2 LGT composed of gauge
fields with only off-diagonal indices Slagle and Kim (2017). On the other
hand, Higgsing the symmetric rank-2 LGT with both diagonal and off-diagonal
gauge fields on the two-dimensional square lattice yields the rank-2 toric
code (R2TC) Ma  _et al._ (2018); Bulmash and Barkeshli (2018); Oh  _et al._
(2022). In contrast to the parent gauge theory hosting immobile fractons or
lineons, the excitations in R2TC are mobile and can hop by $N$ lattice
spacings in the previously forbidden directions due to the condensation of
charge-$N$ excitations in the $\mathbb{Z}_{N}$ gauge theory. An interesting
consequence of it is the dependence of the ground state degeneracy on the
linear system size $(L_{x},L_{y})$ mod $N$, studied in detail in recent works
Oh  _et al._ (2022); Gorantla  _et al._ (2022); Pace and Wen (2022); Delfino
_et al._ (2022).

In this paper, we focus on the other exotic property of R2TC, i.e. mutual
braiding statistics among the quasiparticles that seems to break the
conventional wisdom of Abelian anyon braiding as well as the flux attachment
picture. In R2TC, the braiding phase was found to depend on the initial
positions of the quasiparticles participating in the braiding Oh  _et al._
(2022). Recall that the conventional anyon braiding counts the total charge of
anyons inside the loop, and the statistical phase is encoded in the BF field
theory. We call this framework the monopolar braiding statistics. Here we
develop a new field-theoretic formalism consistent with what we call the
dipolar braiding statistics and develop the associated dipolar BF theory of
braiding. The new field theory is constructed from the rank-2 gauge fields, as
the conventional BF theory is made out of rank-1 gauge fields. The flux
attachment picture is similarly modified to that of dipole attachment.

We begin by giving a quick review of the R2TC model in Sec. II. In Sec. III,
we define this anomalous position dependent statistics as the dipolar braiding
statistics. Construction of effective field theory is done in two steps.
First, in Sec. IV, we construct the effective Lagrangian for matter fields
minimally coupled to the rank-2 gauge fields. Then, in Sec. LABEL:sec:dipolar-
bf-theory, we construct the full dipolar BF theory and derive dipolar braiding
statistics from it. Summary and discussion follows in Sec. VI.

##  II Review of R2TC

The R2TC model is obtained by Higgsing the rank-2 U(1) LGT in two dimensions
employing a pair of canonically conjugate fields $A^{ab}_{i}$ and $E^{ab}_{i}$
obeying $[A,E]=i$ for the same site and indices. There are three independent
components $(A^{xx}_{i},A^{xy}_{i},A^{yy}_{i})$ per site $i$ (similarly for
$E$). Several mutually commuting operators called generators can be
constructed in the rank-2 LGT as:

$\displaystyle G^{x}_{i}=$ $\displaystyle
E_{i+\hat{x}}^{xx}-E_{i}^{xx}+E_{i+\hat{y}}^{xy}-E_{i}^{xy}\equiv\rho_{i}^{x}$
$\displaystyle\sim$
$\displaystyle\partial_{x}E^{xx}_{i}+\partial_{y}E^{xy}_{i},$ $\displaystyle
G^{y}_{i}=$ $\displaystyle
E_{i+\hat{x}}^{xy}-E_{i}^{xy}+E_{i+\hat{y}}^{yy}-E_{i}^{yy}\equiv\rho_{i}^{y}$
$\displaystyle\sim$
$\displaystyle\partial_{x}E^{xy}_{i}+\partial_{y}E^{yy}_{i},$ $\displaystyle
B_{i}=$ $\displaystyle
A^{xx}_{i+\hat{y}}+A^{xx}_{i-\hat{y}}-2A^{xx}_{i}+A^{yy}_{i+\hat{x}}+A^{yy}_{i-\hat{x}}-2A^{yy}_{i}$
$\displaystyle~{}~{}-A^{xy}_{i}+A^{xy}_{i-\hat{x}}+A^{xy}_{i-\hat{y}}-A^{xy}_{i-\hat{x}-\hat{y}}\equiv\frac{2\pi}{N}\rho_{i}^{m}$
$\displaystyle\sim$
$\displaystyle\partial_{y}^{2}A^{xx}_{i}+\partial_{x}^{2}A^{yy}_{i}-\partial_{x}\partial_{y}A^{xy}_{i}.$
(1)

The theory contains a vector of electric charges $(\rho_{i}^{x},\rho_{i}^{y})$
and a scalar magnetic charge $\rho_{i}^{m}$, tied to the gauge fields through
the Gauss law constraints $G_{i}^{a}=\rho_{i}^{a}$ ($a=x,y)$ and
$B_{i}\equiv(2\pi/N)\rho_{i}^{m}$ for some integer $N\geq 2$. They generate
the gauge transformations $A^{ab}_{i}\rightarrow
U_{A}^{\dagger}A_{i}^{ab}U_{A}$ and $E^{ab}_{i}\rightarrow
U_{E}^{\dagger}E_{i}^{ab}U_{E}$ which result in

$\displaystyle A^{aa}_{i}$ $\displaystyle\rightarrow
A^{aa}_{i}+f^{a}_{i}-f^{a}_{i-\hat{a}}~{}~{}(a=x,y)$ $\displaystyle
A^{xy}_{i}$ $\displaystyle\rightarrow
A^{xy}_{i}+f^{y}_{i}-f^{y}_{i-\hat{x}}+f^{x}_{i}-f^{x}_{i-\hat{y}}$
$\displaystyle E^{xx}_{i}$ $\displaystyle\rightarrow
E^{xx}_{i}+g_{i+\hat{y}}+g_{i-\hat{y}}-2g_{i}$ $\displaystyle E^{yy}_{i}$
$\displaystyle\rightarrow E^{yy}_{i}+g_{i+\hat{x}}+g_{i-\hat{x}}-2g_{i}$
$\displaystyle E^{xy}_{i}$ $\displaystyle\rightarrow
E^{xy}_{i}+g_{i+\hat{x}}+g_{i+\hat{y}}-g_{i}-g_{i+\hat{x}+\hat{y}}.$ (2)

The two unitary operators are

$\displaystyle
U_{A}=e^{i\sum_{i}(f_{i}^{x}G_{i}^{x}+f_{i}^{y}G_{i}^{y})},~{}~{}U_{E}=e^{i\sum_{i}g_{i}B_{i}}.$
(3)

Generalized Pauli operators are obtained from Higging the gauge fields as Ma
_et al._ (2018); Bulmash and Barkeshli (2018):

$\displaystyle X_{0,i}=e^{iA_{i}^{xy}},$ $\displaystyle Z_{0,i}=e^{2\pi
iE_{i}^{xy}/N},$ $\displaystyle X_{1,i}=e^{iA_{i}^{xx}},$ $\displaystyle
Z_{1,i}=e^{2\pi iE_{i}^{xx}/N},$ $\displaystyle X_{2,i}=e^{iA_{i}^{yy}},$
$\displaystyle Z_{2,i}=e^{2\pi iE_{i}^{yy}/N}.$ (4)

They operate on the $N$-dimensional local Hilbert space defined by

$\displaystyle
X|g\rangle=|g+1\rangle,~{}Z|g\rangle=\omega^{g}|g\rangle,~{}ZX=\omega XZ,$ (5)

where $|g\rangle=|0\rangle,~{}|1\rangle,~{}\cdots,~{}|N-1\rangle$ and
$\omega=e^{2\pi i/N}$. Applying the Higgsing formulas in Eq. (4) to the three
mutually commuting generators in Eq. (1) results in three types of commuting
spin operators:

$\displaystyle a_{i}^{x}$ $\displaystyle\equiv e^{(2\pi
i/N)G^{x}_{i}}=Z_{1,i}^{-1}Z_{1,i+\hat{x}}Z_{0,i}Z_{0,i-\hat{y}}^{-1},$
$\displaystyle a_{i}^{y}$ $\displaystyle\equiv e^{(2\pi
i/N)G^{y}_{i}}=Z_{2,i}^{-1}Z_{2,i+\hat{x}}Z_{0,i}Z_{0,i-\hat{x}}^{-1},$
$\displaystyle b_{i}$ $\displaystyle\equiv
e^{iB_{i}}=X_{0,i}^{-1}X_{0,i-\hat{x}}X_{0,i-\hat{y}}X_{0,i-\hat{x}-\hat{y}}^{-1}$
$\displaystyle\times
X_{2,i-\hat{x}}X_{2,i}^{-2}X_{2,i+\hat{x}}X_{1,i-\hat{y}}X_{1,i}^{-2}X_{1,i+\hat{y}}.$
(6)

The R2TC Hamiltonian is constructed as

$\displaystyle\mathbb{A}_{i}^{x}$
$\displaystyle=\frac{1}{N}\sum_{j=0}^{N-1}(a_{i}^{x})^{j},$
$\displaystyle\mathbb{A}_{i}^{y}$
$\displaystyle=\frac{1}{N}\sum_{j=0}^{N-1}(a_{i}^{y})^{j},$
$\displaystyle\mathbb{B}_{i}$
$\displaystyle=\frac{1}{N}\sum_{j=0}^{N-1}(b_{i})^{j},$ $\displaystyle H$
$\displaystyle=-\sum_{i}(\mathbb{A}^{x}_{i}+\mathbb{A}^{y}_{i}+\mathbb{B}_{i}).$
(7)

This model has the ground states characterized by $\mathbb{A}_{i}^{x}|{\rm
GS}\rangle=\mathbb{A}_{i}^{y}|{\rm GS}\rangle=\mathbb{B}_{i}|{\rm
GS}\rangle=|{\rm GS}\rangle$. An excited state is obtained when one of the
stabilizers $a_{i}^{x}$, $a_{i}^{y}$, or $b_{i}$ takes on the eigenvalue
$\omega^{j}$ with nonzero integer $j$. We use $e_{x}$, $e_{y}$, and $m$ to
denote the two electric and one magnetic charges associated with $a_{i}^{x}$,
$a_{i}^{y}$, and $b_{i}$ excitations, respectively. The $e_{x}$ ($e_{y}$)
monopole can move freely in the $x$-($y$-) direction, but can hop only by $N$
lattice spacings in the $y$\- ($x$-) direction; the magnetic monopole can only
hop by $N$ spacings in both directions Oh  _et al._ (2022).

##  III Dipolar Braiding Statistics

It was shown Oh  _et al._ (2022) that non-trivial braiding statistics arises
between $(e_{x},m)$ and $(e_{y},m)$ monopoles111We often refer to
$e_{x},e_{y},m$ quasiparticle excitations as ‘monopoles’ in order to
distinguish them from dipole excitations that are composites of the monopoles.
with phase factors given by

$\displaystyle
e^{i\phi_{e_{x},m}}=\omega^{(y_{x0}-y_{m0})},~{}~{}e^{i\phi_{e_{y},m}}=\omega^{(x_{m0}-x_{y0})}.$
(8)

The $e_{x}$, $e_{y}$, and $m$ monopoles are located at

$\displaystyle{\bm{r}}_{x0}=(x_{x0},y_{x0}),~{}{\bm{r}}_{y0}=(x_{y0},y_{y0}),~{}{\bm{r}}_{m}=(x_{m0},y_{m0})$
(9)

at the start of the braiding, respectively.

The following integral representations of the phase factors can be given for
the braiding of $e_{x}$ and $e_{y}$ particle around $m$:

$\displaystyle\phi_{e_{x},m}$ $\displaystyle=\frac{2\pi}{N}\int
d^{2}{\bm{r}}~{}(y_{x0}-y)\rho^{m}({\bm{r}}),$ $\displaystyle\phi_{e_{y},m}$
$\displaystyle=\frac{2\pi}{N}\int
d^{2}{\bm{r}}~{}(x-x_{y0})\rho^{m}({\bm{r}}).$ (10)

The magnetic monopole density inside the contour is $\rho^{m}({\bm{r}})$. The
phase factors in Eq. (8) is recovered for a single $m$ monopole
$\rho^{m}({\bm{r}})=\delta^{2}({\bm{r}}-{\bm{r}}_{m0})$. The initial position
of the $e_{x}$ or $e_{y}$ particle that does the braiding enters explicitly in
the above integrand. The familiar anyonic statistics is obtained instead from
integrating the monopole density; in contrast Eq. (10) can be thought to
integrate the dipole density of the enclosed particles.

The phases $\phi_{m,e_{x}}$ and $\phi_{m,e_{y}}$ for braiding the $m$ monopole
around $e_{x}$ and $e_{y}$ monopole, respectively, have similar integral
expressions:

$\displaystyle\phi_{m,e_{x}}=$ $\displaystyle\frac{2\pi}{N}\int
d^{2}{\bm{r}}(y-y_{m0})\rho^{x}({\bm{r}}),$ $\displaystyle\phi_{m,e_{y}}=$
$\displaystyle\frac{2\pi}{N}\int d^{2}{\bm{r}}(x_{m0}-x)\rho^{y}({\bm{r}}).$
(11)

The densities of $e_{x}$ and $e_{y}$ are $\rho^{x}({\bm{r}})$ and
$\rho^{y}({\bm{r}})$, respectively.

A lattice-model interpretation of the integral formulas in Eqs. (10) and (11)
can be given. Discretizing the phase integrals in Eqs. (10) and (11) and
invoking the Higgsing formula in Eq. (6),

$\displaystyle e^{i\phi_{e_{x},m}}$ $\displaystyle\rightarrow
W_{e_{x},m}=\exp\left(i\frac{2\pi}{N}\sum_{i\in{\cal
A}}(y_{x0}-y_{i})\rho_{i}^{m}\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=\prod_{i\in\mathcal{A}}\left(b_{i}\right)^{y_{x0}-y_{i}}$
$\displaystyle e^{i\phi_{e_{y},m}}$ $\displaystyle\rightarrow
W_{e_{y},m}=\exp\left(i\frac{2\pi}{N}\sum_{i\in{\cal
A}}(x_{i}-x_{y0})\rho_{i}^{m}\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=\prod_{i\in\mathcal{A}}\left(b_{i}\right)^{x_{i}-x_{y0}}$
$\displaystyle e^{i\phi_{m,e_{x}}}$ $\displaystyle\rightarrow
W_{m,e_{x}}=\prod_{i\in{\cal A}}\left(a_{i}^{x}\right)^{y_{i}-y_{m0}},$
$\displaystyle e^{i\phi_{m,e_{y}}}$ $\displaystyle\rightarrow
W_{m,e_{y}}=\prod_{i\in{\cal A}}\left(a_{i}^{y}\right)^{x_{m0}-x_{i}}.$ (12)

Each operator $W$ is given as a product of various stabilizers over the area
$\prod_{i\in\mathcal{A}}$ inside the braiding path. When acting on the ground
state this yields 1, but not so if some monopoles reside in the area.
Inserting explicit expressions for $b_{i}$ from Eq. (6) into $W_{e_{x},m}$ and
$W_{e_{y},m}$ shows that all the operators in the interior of ${\cal A}$
cancel out, leaving only the product of $X$ operators along the boundary hence
qualifying them as Wegner-Wilson operators of sorts.

In the case of a rectangular boundary ${\cal
A}=[x_{0},x_{1}]\times[y_{0},y_{1}]$ one finds

$\displaystyle W_{e_{x},m}$ $\displaystyle=T_{x}(x_{0}\rightarrow
x_{1}|y_{0})T_{x}(y_{0}\rightarrow y_{1}|x_{1})\left[T_{x}(x_{0}\rightarrow
x_{1}|y_{1})\right]^{-1}\left[T_{x}(y_{0}\rightarrow y_{1}|x_{0})\right]^{-1}$
$\displaystyle W_{e_{y},m}$ $\displaystyle=T_{y}(x_{0}\rightarrow
x_{1}|y_{0})T_{y}(y_{0}\rightarrow y_{1}|x_{1})\left[T_{y}(x_{0}\rightarrow
x_{1}|y_{1})\right]^{-1}\left[T_{y}(y_{0}\rightarrow y_{1}|x_{0})\right]^{-1}$
(13)

where the $e_{x},e_{y}$ translation operators are defined as

$\displaystyle
T_{x}(x_{1}\\!\rightarrow\\!x_{2}|y)\\!=\\!\prod_{x=x_{1}}^{x_{2}-1}\left(X_{1,i}\right)^{-1}\left(X_{1,i}X_{1,i\\!-\\!\hat{y}}^{-1}\right)^{y\\!-\\!y_{x0}}\\!=\\!\exp\left(i\sum_{x=x_{1}}^{x_{2}-1}\left[-A^{xx}(x,y)\\!+\\!(y\\!-\\!y_{x0})(A^{xx}(x,y)\\!-\\!A^{xx}(x,y\\!-\\!1))\right]\right),$
$\displaystyle
T_{x}(y_{1}\\!\rightarrow\\!y_{2}|x)\\!=\\!\prod_{y=y_{1}}^{y_{2}-1}\\!\left(X_{0,i}X_{0,i\\!-\\!\hat{y}}^{-1}X_{2,i}^{-1}X_{2,i\\!-\\!\hat{x}}\right)^{y\\!-\\!y_{x0}}\\!=\\!\exp\\!\left(\\!i\sum_{y\\!=\\!y_{1}}^{y_{2}-1}\\!\left(y\\!-\\!y_{x0}\right)\\!\left[A^{xy}(x,\\!y)\\!-\\!A^{xy}(x,\\!y\\!-\\!1)\\!-\\!A^{yy}(x,\\!y)\\!+\\!A^{yy}(x\\!-\\!1,\\!y)\right]\\!\right),$
$\displaystyle
T_{y}(x_{1}\\!\rightarrow\\!x_{2}|y)\\!=\\!\prod_{x=x_{1}}^{x_{2}-1}\\!\left(X_{0,i}X_{0,i\\!-\\!\hat{x}}^{-1}X_{1,i}^{-1}X_{1,i\\!-\\!\hat{y}}\right)^{x\\!-\\!x_{y0}}\\!=\\!\exp\\!\left(\\!i\sum_{x\\!=\\!x_{1}}^{x_{2}-1}\\!\left(x\\!-\\!x_{y0}\right)\\!\left[A^{xy}(x,\\!y)\\!-\\!A^{xy}(x\\!-\\!1,\\!y)\\!-\\!A^{xx}(x,\\!y)\\!+\\!A^{xx}(x,\\!y\\!-\\!1)\right]\\!\right),$
$\displaystyle
T_{y}(y_{1}\\!\rightarrow\\!y_{2}|x)\\!=\\!\prod_{y=y_{1}}^{y_{2}-1}\left(X_{2,i}\right)^{-1}\left(X_{2,i}X_{2,i\\!-\\!\hat{x}}^{-1}\right)^{x-x_{y0}}\\!=\\!\exp\left(i\sum_{y=y_{1}}^{y_{2}-1}\left[-A^{yy}(x,y)\\!+\\!(x\\!-\\!x_{y0})(A^{yy}(x,y)\\!-\\!A^{yy}(x\\!-\\!1,y))\right]\right).$
(14)

The $T_{x}(x_{1}\\!\rightarrow\\!x_{2}|y)$ translates the $e_{x}$
quasiparticle from $(x_{1},y)$ to $(x_{2},y)$. On the other hand,
$T_{x}(y_{1}\\!\rightarrow\\!y_{2}|x)$ translates $e_{x}$ from $(x,y_{1})$ to
$(x,y_{2})$ only if $y_{1}-y_{x0}~{}{\rm mod}~{}N=0$ and $y_{2}-y_{x0}~{}{\rm
mod}~{}N=0$ Oh  _et al._ (2022). Recall that the $e_{x}$ monopole is able to
hop only by $N$ lattice sites in the $y$-direction. Nevertheless, one can
define the operator $T_{x}(y_{1}\\!\rightarrow\\!y_{2}|x)$ for any pair of
coordinates $y_{1},y_{2}$, under a new interpretation: for general
$y_{1},y_{2}$, the translation of $e_{x}$ embodies the simultaneous creation
and motion of some auxiliary dipoles that are necessary to maintain the
overall dipole moment conservation mod $N$. Similar statements apply to
$T_{y}(x_{1}\rightarrow x_{2}|y)$. Details can be found in Appendix A.

On taking continuum limits of Eq. (14), the braiding phases become line
integrals

$\displaystyle\phi_{e_{x},m}$ $\displaystyle=\oint
d{\bm{r}}\cdot{\bm{a}}_{x},$ $\displaystyle\phi_{e_{y},m}$
$\displaystyle=\oint d{\bm{r}}\cdot{\bm{a}}_{y}$ (15)

with a pair of emergent vector potentials defined by

$\displaystyle{\bm{a}}_{x}=$
$\displaystyle\bigl{(}\\!-\\!A^{xx}\\!+\\!\left(y\\!-\\!y_{x0}\right)\\!\partial_{y}A^{xx},\\!(y\\!-\\!y_{x0})\\!\left(\partial_{y}A^{xy}\\!-\\!\partial_{x}A^{yy}\right)\bigr{)},$
$\displaystyle{\bm{a}}_{y}=$
$\displaystyle\bigl{(}\\!(x\\!-\\!x_{y0})\\!\left(\partial_{x}A^{xy}\\!-\\!\partial_{y}A^{xx}\right),\\!-\\!A^{yy}\\!+\\!\left(x\\!-\\!x_{y0}\right)\\!\partial_{x}A^{yy}\bigr{)}.$
(16)

The meaning of ${\bm{a}}_{x}$ and ${\bm{a}}_{y}$ is clarified on taking their
curl,

$\displaystyle{\bm{\nabla}}\times{\bm{a}}_{x}=$
$\displaystyle\frac{2\pi}{N}(y_{x0}-y)\rho^{m}$
$\displaystyle{\bm{\nabla}}\times{\bm{a}}_{y}=$
$\displaystyle\frac{2\pi}{N}(x-x_{y0})\rho^{m},$ (17)

where we used the constraint

$\displaystyle\partial_{x}^{2}A^{yy}+\partial_{y}^{2}A^{xx}-\partial_{x}\partial_{y}A^{xy}=(2\pi/N)\rho^{m}.$

In an ordinary flux attachment scenario, the curl of the emergent vector
potential equals the charge density. Here it gives either the $x$ or the $y$
component of the dipole density. The initial positions of the $e_{x}$ and
$e_{y}$ monopoles appear explicitly in the formulas as a result.

Next, we turn to the braiding operators $W_{m,e_{x}}$ and $W_{m,e_{y}}$. It
turns out $W_{e_{x},m}$ and $W_{e_{y},m}$ individually do not have the nice
cancellation of bulk terms but their product $W_{m,e}\equiv
W_{m,e_{x}}W_{m,e_{y}}$ does, and become a product of boundary operators

$\displaystyle W_{m,e}$ $\displaystyle=T_{m}(x_{0}\rightarrow
x_{1}|y_{0})T_{m}(y_{0}\rightarrow y_{1}|x_{1})\left[T_{m}(x_{0}\rightarrow
x_{1}|y_{1})\right]^{-1}\left[T_{m}(y_{0}\rightarrow
y_{1}|x_{0})\right]^{-1},$ (18)

with

$\displaystyle T_{m}(x_{1}\rightarrow
x_{2}|y)=\prod_{x=x_{1}\\!+\\!1}^{x_{2}}\left(Z_{2,i}\right)^{x\\!-\\!x_{m0}}\left(Z_{0,i-\hat{x}}\right)^{y_{m0}\\!-\\!y}\\!=\\!\exp\left[i\frac{2\pi}{N}\sum_{x=x_{1}\\!+\\!1}^{x_{2}}\bigl{(}(x\\!-\\!x_{m0})E^{yy}(x,y)\\!-\\!(y\\!-\\!y_{m0})E^{xy}(x\\!-\\!1,y)\bigr{)}\right],$
$\displaystyle T_{m}(y_{1}\rightarrow
y_{2}|x)=\prod_{y=y_{1}\\!+\\!1}^{y_{2}}\left(Z_{1,i}\right)^{y\\!-\\!y_{m0}}\\!\left(Z_{0,i-\hat{y}}\right)^{x_{m0}\\!-\\!x}\\!=\\!\exp\\!\left[i\frac{2\pi}{N}\\!\sum_{y=y_{1}\\!+\\!1}^{y_{2}}\\!\bigl{(}(y\\!-\\!y_{m0})E^{xx}(x,y)\\!-\\!(x\\!-\\!x_{m0})E^{xy}(x,y\\!-\\!1)\bigr{)}\\!\right].$
(19)

The overall braiding phase $\phi_{m,e}=\phi_{m,e_{x}}+\phi_{m,e_{y}}$ for $m$
around both $e_{x}$ and $e_{y}$ can be written $\phi_{m,e}=\oint
d{\bm{r}}\cdot{\bm{a}}_{m}$, where

$\displaystyle{\bm{a}}_{m}\\!=\\!\frac{2\pi}{N}$
$\displaystyle\bigl{(}(x\\!-\\!x_{m0})E^{yy}\\!-\\!(y\\!-\\!y_{m0})E^{xy},$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(y\\!-\\!y_{m0})E^{xx}\\!-\\!(x\\!-\\!x_{m0})E^{xy}\bigr{)},$
$\displaystyle{\bm{\nabla}}\times{\bm{a}}_{m}$
$\displaystyle=\frac{2\pi}{N}\left[(y-y_{m0})\rho^{x}-(x-x_{m0})\rho^{y}\right].$
(20)

The constraints

$\displaystyle\partial_{x}E^{xx}\\!+\\!\partial_{y}E^{xy}=\rho^{x},~{}~{}\partial_{x}E^{xy}\\!+\\!\partial_{y}E^{yy}=\rho^{y}$

were used in the last line.

To sum up, a dipolar AB phase formulation of the dipolar braiding statistics
is possible in terms of some emergent vector potential ${\bm{a}}_{x}$,
${\bm{a}}_{y}$ and ${\bm{a}}_{m}$, which are functions of the underlying
rank-2 fields $A^{ab},E^{ab}$ as well as the initial locations of the
quasiparticles. For completeness we mention that under the gauge
transformation

$\displaystyle A^{xx}\rightarrow$ $\displaystyle A^{xx}+\partial_{x}f_{x},$
$\displaystyle E^{xx}\rightarrow$ $\displaystyle E^{xx}+\partial_{y}^{2}g,$
$\displaystyle A^{xy}\rightarrow$ $\displaystyle
A^{xy}+\partial_{x}f_{y}\\!+\\!\partial_{y}f_{x},$ $\displaystyle
E^{xy}\rightarrow$ $\displaystyle E^{xy}-\partial_{x}\partial_{y}g,$
$\displaystyle A^{yy}\rightarrow$ $\displaystyle A^{yy}+\partial_{y}f_{y},$
$\displaystyle E^{yy}\rightarrow$ $\displaystyle E^{yy}+\partial_{x}^{2}g,$

of the rank-2 fields, the emergent vector potentials transform as

$\displaystyle{\bm{a}}_{x}\rightarrow{\bm{a}}_{x}+{\bm{\nabla}}F_{x},{\bm{a}}_{y}\rightarrow{\bm{a}}_{y}+{\bm{\nabla}}F_{y},{\bm{a}}_{m}\rightarrow{\bm{a}}_{m}+{\bm{\nabla}}G,$

where

$\displaystyle(F_{x},F_{y})$
$\displaystyle=\bigl{(}-f_{x}+(y-y_{x0})\partial_{y}f_{x},f^{y}-(x-x_{y0})\partial_{x}f_{y}\bigr{)}$
$\displaystyle G$
$\displaystyle=(x-x_{m0})\partial_{x}g+(y-y_{m0})\partial_{y}g-g.$

##  IV Matter Fields with Minimal Coupling

There are three matter fields in the R2TC, represented by field operators
$(\psi_{x},\psi_{y},\psi_{m})$, which are related to the monopole density by
$\rho^{a}=\psi_{a}^{\dagger}\psi_{a}$ $(a=x,y,m)$. Appropriate field theory
for these can be constructed by exploiting the gauge transformation properties
in Eqs. (2) and (3) and identifying proper covariant derivatives of the matter
fields.

Taking the continuum limit of the unitary operators in Eq. (3) and invoking
the constraints $G_{i}^{a}=\rho_{i}^{a}$ and $B_{i}=(2\pi/N)\rho_{i}^{m}$ from
Eq. (1), one can show that the matter fields $\psi_{a}$ $(a=x,y,m)$ and the
gauge fields $(A^{ab},E^{ab})$ transform as

$\displaystyle(\psi_{x},\psi_{y},\psi_{m})$
$\displaystyle\rightarrow(e^{if_{x}}\psi_{x},e^{if_{y}}\psi_{y},e^{i\frac{2\pi}{N}g}\psi_{m})$
$\displaystyle A^{ab}$ $\displaystyle\rightarrow
A^{ab}\\!+\\!(\partial_{x}f_{x},\partial_{x}f_{y}\\!+\\!\partial_{y}f_{x},\partial_{y}f_{y})$
$\displaystyle E^{ab}$ $\displaystyle\rightarrow
E^{ab}\\!+\\!(\partial_{y}^{2}g,-\partial_{x}\partial_{y}g,\partial_{x}^{2}g).$
(21)

Identifying the covariant derivatives and constructing a gauge-invariant
Lagrangian ia a straightforward exercise Kim  _et al._ (2022). For $e_{x}$
$(\psi_{x})$ and $e_{y}$ ($\psi_{y}$) particles they are

$\displaystyle D_{x}\psi_{x}\equiv(\partial_{x}-iA^{xx})\psi_{x},$
$\displaystyle D_{y}\psi_{y}\equiv(\partial_{y}-iA^{yy})\psi_{y},$
$\displaystyle
D_{xy}\psi_{xy}\equiv\psi_{x}\partial_{x}\psi_{y}+\psi_{y}\partial_{y}\psi_{x}-iA^{xy}\psi_{x}\psi_{y},$
(22)

and

$\displaystyle L_{e}=\sum_{a=x,y}\psi_{a}^{\dagger}\partial_{t}\psi_{a}$
$\displaystyle-\frac{1}{2}\left(\alpha_{e}\left|D_{x}\psi_{x}\right|^{2}+\alpha_{e}\left|D_{y}\psi_{y}\right|^{2}+\beta_{e}\left|D_{xy}\psi_{xy}\right|^{2}\right).$
(23)

Similarly for the $m$ particles,

$\displaystyle
D_{xx}\psi_{m}\equiv\psi_{m}\partial_{x}^{2}\psi_{m}-(\partial_{x}\psi_{m})^{2}-i\frac{2\pi}{N}E^{yy}\psi_{m}^{2},$
$\displaystyle
D_{yy}\psi_{m}\equiv\psi_{m}\partial_{y}^{2}\psi_{m}-(\partial_{y}\psi_{m})^{2}-i\frac{2\pi}{N}E^{xx}\psi_{m}^{2},$
$\displaystyle
D_{xy}\psi_{m}\equiv\psi_{m}\partial_{x}\partial_{y}\psi_{m}-(\partial_{x}\psi_{m})(\partial_{y}\psi_{m})+i\frac{2\pi}{N}E^{xy}\psi_{m}^{2}$
(24)

and

$\displaystyle L_{m}=i\psi_{m}^{\dagger}\partial_{t}\psi_{m}$
$\displaystyle-\frac{1}{4}\left(\alpha_{m}\left|D_{xx}\psi_{m}\right|^{2}+\alpha_{m}\left|D_{yy}\psi_{m}\right|^{2}+\beta_{m}\left|D_{xy}\psi_{m}\right|^{2}\right).$
(25)

define the appropriate matter field theory. Various constants ($\alpha$’s
$\beta$’s) appear, which do not change the generic features of the field
theory.

Continuity equations for the quasiparticles particles follow as

$\displaystyle\partial_{t}{\rho}^{x}+\partial_{x}J_{e}^{xx}+\partial_{y}J_{e}^{xy}$
$\displaystyle=0,$
$\displaystyle\partial_{t}{\rho}^{y}+\partial_{x}J_{e}^{xy}+\partial_{y}J_{e}^{yy}$
$\displaystyle=0,$
$\displaystyle\partial_{t}{\rho}^{m}+\partial_{x}^{2}J_{m}^{xx}+\partial_{x}\partial_{y}J_{m}^{xy}+\partial_{y}^{2}J_{m}^{yy}$
$\displaystyle=0,$ (26)

where

$\displaystyle J_{e}^{aa}\equiv$
$\displaystyle-\frac{i}{2}\alpha_{e}\Big{[}\psi_{a}^{\dagger}\partial_{a}\psi_{a}-\left(\partial_{a}\psi_{a}^{\dagger}\right)\psi_{a}\Big{]},~{}~{}(a=x,y)$
$\displaystyle J_{e}^{xy}\equiv$
$\displaystyle-\frac{i}{2}\beta_{e}\psi_{x}^{\dagger}\psi_{y}^{\dagger}\left(\psi_{x}\partial_{x}\psi_{y}+\psi_{y}\partial_{y}\psi_{x}\right)+h.c.=J_{e}^{yx},$
$\displaystyle J_{m}^{aa}=$
$\displaystyle\\!-\\!\frac{i}{4}\alpha_{m}\\!\left(\psi_{m}^{\dagger}\partial_{a}^{2}\psi_{m}^{\dagger}\\!-\\!\partial_{a}\psi_{m}^{\dagger}\partial_{a}\psi_{m}^{\dagger}\right)\psi_{m}^{2}\\!+\\!h.c.~{}(a\\!=\\!x,y)$
$\displaystyle J_{m}^{xy}=$
$\displaystyle\\!-\\!\frac{i}{4}\beta_{m}\\!\left(\psi_{m}^{\dagger}\partial_{x}\partial_{y}\psi_{m}^{\dagger}\\!-\\!\partial_{x}\psi_{m}^{\dagger}\partial_{y}\psi_{m}^{\dagger}\right)\psi_{m}^{2}+h.c.\\!=\\!J_{m}^{yx}.$
(27)

Both the covariant derivatives and the conservation laws depart significantly
from those of particles in the rank-1 U(1) gauge fields.

Several conserved quantities can be identified from the continuity equations.
Writing the three monpole charges $Q^{a}=\int d^{2}{\bm{r}}\rho^{a}$
$(a=x,y,m)$, and three dipole charges

$\displaystyle{\bm{\mu}}^{m}$ $\displaystyle=\int
d^{2}{\bm{r}}{\bm{r}}\rho^{m},$ $\displaystyle{\mu}^{e}$ $\displaystyle=\int
d^{2}{\bm{r}}(x\rho^{y}-y\rho^{x}),$ (28)

one can show their time derivatives vanish identicallly under the appropriate
boundary conditions at infinity. The existence of six conserved quantities in
the theory is closely tied to the existence of (up to) six independent
holonomies and (up to) $N^{6}$ ground state degeneracy in the R2TC222There is
an intricate relation between the six conserved quantities mentioned here and
the holonomies that generate the ground state degeneracy. They will be
discussed in an upcoming article Oh et al. (2022).

Finally, the covariant derivatives and the accompanying Lagrangians we
constructed can be cast in the lattice mode. Writing
($\psi_{x},\psi_{y},\psi_{m}\rightarrow X_{i},Y_{i},M_{i}$), we obtain

$\displaystyle|D_{x}\psi_{x}|^{2}$ $\displaystyle\rightarrow
X^{\dagger}_{i}X_{i-\hat{x}}e^{iA^{xx}_{i}}+h.c.$
$\displaystyle|D_{y}\psi_{y}|^{2}$ $\displaystyle\rightarrow
Y^{\dagger}_{i}Y_{i-\hat{y}}e^{iA^{yy}_{i}}+h.c.$
$\displaystyle|D_{xy}\psi_{xy}|^{2}$ $\displaystyle\rightarrow
Y^{\dagger}_{i}Y_{i-\hat{x}}X^{\dagger}_{i}X_{i-\hat{y}}e^{iA^{xy}_{i}}+h.c.$
$\displaystyle|D_{xx}\psi_{m}|^{2}$ $\displaystyle\rightarrow
M_{i}^{2}M^{\dagger}_{i+\hat{x}}M^{\dagger}_{i-\hat{x}}e^{iE^{xx}_{i}}+h.c.$
$\displaystyle|D_{yy}\psi_{m}|^{2}$ $\displaystyle\rightarrow
M_{i}^{2}M^{\dagger}_{i+\hat{y}}M^{\dagger}_{i-\hat{y}}e^{iE^{yy}_{i}}+h.c.$
$\displaystyle|D_{xy}\psi_{m}|^{2}$ $\displaystyle\rightarrow
M_{i}M_{i+\hat{x}+\hat{y}}M^{\dagger}_{i+\hat{x}}M^{\dagger}_{i+\hat{y}}e^{iE^{xy}_{i}}+h.c.$
(29)

Similar expressions have been suggested in the tight-binding model of higher-
order topological insulator You  _et al._ (2021).

##  V Dipolar BF theory

Comparing the temporal derivatives of the constraints in Eq. (1) with the
continuity equation in Eq. (26) allows the identification of several
identities:

$\displaystyle
J_{e}^{ab}=-\partial_{t}{E}^{ab},\frac{2\pi}{N}J_{m}^{xx}=-\partial_{t}{A}^{yy},$
$\displaystyle\frac{2\pi}{N}J_{m}^{xy}=\partial_{t}{A}^{xy},\frac{2\pi}{N}J_{m}^{yy}=-\partial_{t}{A}^{xx}.$
(30)

The Lagrangian that encodes all these constraints as well as the commutation
$[A^{ab}({\bm{r}}),E^{ab}({\bm{r}}^{\prime})]=i\delta^{2}({\bm{r}}-{\bm{r}}^{\prime})$
is

$\displaystyle{\cal L}_{\rm dBF}$
$\displaystyle=A^{0x}\left(\partial_{x}E^{xx}\\!+\\!\partial_{y}E^{xy}-\rho^{x}\right)+A^{0y}\left(\partial_{x}E^{xy}\\!+\\!\partial_{y}E^{yy}-\rho^{y}\right)+E^{0}\left(\partial_{x}^{2}A^{yy}\\!+\\!\partial_{y}^{2}A^{xx}\\!-\\!\partial_{x}\partial_{y}A^{xy}-\frac{2\pi}{N}\rho^{m}\right)$
$\displaystyle+E^{xx}\left(\partial_{t}{A}^{xx}\\!+\\!\frac{2\pi}{N}J_{m}^{yy}\right)+E^{xy}\left(\partial_{t}{A}^{xy}\\!-\\!\frac{2\pi}{N}J_{m}^{xy}\right)+E^{yy}\left(\partial_{t}{A}^{yy}\\!+\\!\frac{2\pi}{N}J_{m}^{xx}\right)-A^{xx}J_{e}^{xx}\\!-\\!A^{xy}J_{e}^{xy}\\!-\\!A^{yy}J_{e}^{yy}.$
(31)

This action is the dipolar BF theory (dBF) in contrast to the BF theory of the
rank-1 TC. The Lagrange multipliers transform as

$\displaystyle A^{0a}$
$\displaystyle\\!\rightarrow\\!A^{0a}\\!+\\!\partial_{t}f^{a},$ $\displaystyle
E^{0}\\!\rightarrow\\!E^{0}\\!+\\!\partial_{t}g.$ (32)

The conservation laws of Eq. (26) are recovered from the gauge invariance of
the action.

The dipolar AB phase factors follow straightforwardly from the dBF action, but
before doing so one must first address the important conceptual question: how
can one adiabatically move a particle when its motion is constrained? For
instance, the $m$ particle cannot move at all in the limit
$N\rightarrow\infty$ due to the dipole conservation, and $e_{x}$ ($e_{y}$) can
only move smoothly in the $x$ ($y$) direction. Nevertheless, by making use of
the translation operators defined in Eqs. (14) and (19), it seems as though
the monopoles can be freely moved around. Upon careful analysis of the lattice
translation operators, we realize that such motion becomes possible at the
expense of creating an auxiliary dipole that compensates for the changes in
the dipole moment during the monopole motion, as explicitly demonstrated in
Appendix B. The composite particle of a monopole and an auxiliary dipole
carries a net dipole moment of zero, and are free from constraint.

We show how to derive the dipolar braiding phase explicitly in the the
continuum theory. As an $m$ particle initially at
${\bm{r}}_{m0}=(x_{m0},y_{m0})$ is moved to
${\bm{r}}_{m}(t)=(x_{m}(t),y_{m}(t))$, the monopole and auxiliary dipole
composite has the density

$\displaystyle\rho^{m}({\bm{r}})=[1+({\bm{r}}_{m}(t)-{\bm{r}}_{m0})\cdot{\bm{\nabla}}]\delta^{2}({\bm{r}}-{\bm{r}}_{m}(t)).$
(33)

The second part reflects the dipole contribution. Details of the derivation
can be found in Appendix B. The total dipole moment ${\bm{\mu}}_{m}=\int
d^{2}{\bm{r}}{\bm{r}}\rho^{m}({\bm{r}})$ is indeed conserved. When the
position ${\bm{r}}_{m}$ varies over time, the time derivative of the total
density $\rho^{m}$ becomes

$\displaystyle\partial_{t}\rho^{m}({\bm{r}})=$
$\displaystyle-\dot{x}_{m}\left[(x_{m}-x_{m0})\partial_{x}^{2}\delta^{2}({\bm{r}}-{\bm{r}}_{m})+(y_{m}-y_{m0})\partial_{x}\partial_{y}\delta^{2}({\bm{r}}-{\bm{r}}_{m})\right]$
$\displaystyle-\dot{y}_{m}\left[(x_{m}-x_{m0})\partial_{x}\partial_{y}\delta^{2}({\bm{r}}-{\bm{r}}_{m})+(y_{m}-y_{m0})\partial_{y}^{2}\delta^{2}({\bm{r}}-{\bm{r}}_{m})\right].$
(34)

Invoking the continuity equation for $\rho^{m}$ in Eq. (26), one can deduce
the accompanying current density

$\displaystyle J_{m}^{xx}$
$\displaystyle=\dot{x}_{m}(x_{m}-x_{m0})\delta^{2}({\bm{r}}-{\bm{r}}_{m}),$
$\displaystyle J_{m}^{xy}$
$\displaystyle=\left(\dot{x}_{m}(x_{m}-x_{m0})+\dot{y}_{m}(y_{m}-y_{m0})\right)\delta^{2}({\bm{r}}-{\bm{r}}_{m}),$
$\displaystyle J_{m}^{yy}$
$\displaystyle=\dot{y}_{m}(y_{m}-y_{m0})\delta^{2}({\bm{r}}-{\bm{r}}_{m}).$
(35)

Finally, inserting the obtained current density into the dBF action gives

$\displaystyle\frac{2\pi}{N}\int dt\int
d^{2}{\bm{r}}\left[E^{xx}J_{m}^{yy}-E^{xy}J_{m}^{xy}+E^{yy}J_{m}^{xx}\right]$
$\displaystyle=\oint d{\bm{r}}_{m}\cdot{\bm{a}}_{m}({\bm{r}_{m}}),$ (36)

precisely equal to the statistical phase $\phi_{m,e}$ from braiding $m$ around
$e_{x},e_{y}$ quasiparticles.

Likewise when the $e_{x}$ monopole is translated from
${\bm{r}}_{x0}=(x_{x0},y_{x0})$ to ${\bm{r}}_{x}(t)=(x_{x}(t),y_{x}(t))$, or
the $e_{y}$ monopole from ${\bm{r}}_{y0}=(x_{y0},y_{y0})$ to
${\bm{r}}_{y}(t)=(x_{y}(t),y_{y}(t))$, the net monopole and auxiliary dipole
composite has the density

$\displaystyle\rho^{x}({\bm{r}})$
$\displaystyle=\delta^{2}({\bm{r}}\\!-\\!{\bm{r}}_{x}(t))+(y_{x}(t)\\!-\\!y_{x0})\partial_{y}\delta^{2}({\bm{r}}-{\bm{r}}_{x}(t)),$
$\displaystyle\rho^{y}({\bm{r}})$
$\displaystyle=\delta^{2}({\bm{r}}\\!-\\!{\bm{r}}_{y}(t))+(x_{y}(t)\\!-\\!x_{y0})\partial_{x}\delta^{2}({\bm{r}}\\!-\\!{\bm{r}}_{y}(t)).$
(37)

The net electric dipole moment $\mu_{e}=\int
d^{2}{\bm{r}}[x\rho^{y}({\bm{r}})-y\rho^{x}({\bm{r}})]$ is conserved during
the motion.

The current density accompanying the adiabatic motion of $e_{x}$ is

$\displaystyle J^{xx}_{e}=$
$\displaystyle\dot{x}_{x}\left[1+(y_{x}-y_{x0})\partial_{y}\right]\delta^{2}({\bm{r}}-{\bm{r}}_{x}),$
$\displaystyle J^{xy}_{e}=$
$\displaystyle\dot{y}_{x}(y_{x}-y_{x0})\partial_{y}\delta^{2}({\bm{r}}-{\bm{r}}_{x}),$
$\displaystyle J^{yy}_{e}=$
$\displaystyle-\dot{y}_{x}(y_{x}-y_{x0})\partial_{x}\delta^{2}({\bm{r}}-{\bm{r}}_{x}),$
(38)

while that of the $e_{y}$ is

$\displaystyle J^{xx}_{e}=$
$\displaystyle-\dot{x}_{y}(x_{y}-x_{y0})\partial_{y}\delta^{2}({\bm{r}}-{\bm{r}}_{y}),$
$\displaystyle J^{xy}_{e}=$
$\displaystyle\dot{x}_{y}(x_{y}-x_{y0})\partial_{x}\delta^{2}({\bm{r}}-{\bm{r}}_{y}),$
$\displaystyle J^{yy}_{e}=$
$\displaystyle\dot{y}_{y}\left[1+(x_{y}-x_{y0})\partial_{x}\right]\delta^{2}({\bm{r}}-{\bm{r}}_{y}).$
(39)

Inserting the current densities of the $e_{x}$ or $e_{y}$ monopoles in the dBF
action results in

$\displaystyle-\int dt\int
dxdy\left[A^{xx}J_{e}^{xx}+A^{xy}J_{e}^{xy}+A^{yy}J_{e}^{yy}\right]$ (40)

equal to either $\oint
d{\bm{r}}_{x}\cdot{\bm{a}}_{y}({\bm{r}}_{x})=\phi_{e_{x},m}$ or $\oint
d{\bm{r}}_{y}\cdot{\bm{a}}_{y}({\bm{r}}_{y})=\phi_{e_{y},m}$. The dBF action
is thus consistent with the dipolar braiding phases obtained from the analysis
of the lattice model. Importantly, both the monopole and the auxiliary dipole
components in the density must be kept in order to derive the correct braiding
statistics.

##  VI Discussion

Field-theoretic formulation for the dipolar braiding statistics first found in
the R2TC model is developed. The statistical phases count the net dipole
moment rather than the charge of the quasiparticles enclosed in the braiding.
The adiabatic motion of a quasiparticle is possible when accompanied by the
simultaneous motion of auxiliary dipole, which helps conserve the total dipole
moment. The dipolar braiding statistics is derived by solving the equation of
motion for quasiparticles.

Various field theories of dipolar nature have been proposed in the past in
three dimensions You  _et al._ (2020, 2021). Our paper proposes a dipolar
field theory in two dimensions along with the accompanying dipolar statistics.
While this work was under review, several papers analyzed the ground state
degeneracy and the dipolar braiding statistics from a complementary
perspective of anyon lattice Pace and Wen (2022); Delfino  _et al._ (2022),
and proposed a Chern-Simons-type field (Ref. Pace and Wen (2022)) or a
combination of Chern-Simons and dipolar Chern-Simons theory descriptions (Ref.
Delfino  _et al._ (2022)). It remains to explore how the dipolar BF theory
proposed in this work relates to these other theories. An important part of
the specification of the topological field theory is that of the global
structure under the so-called large gauge transformation. We have not
completely worked out the global structure of the dipolar BF theory or the
related ground state degeneracy count in this paper. These problems will be
addressed in a forthcoming article Oh et al. (2022).

###### Acknowledgements.

Y.-T. O. was supported by National Research Foundation of Korea under Grant
NRF-2014R1A6A1030732, NRF-2021R1A2C4001847, NRF-2020R1A4A3079707, and
NRF-2022R1I1A1A01065149. JHH was supported by NRF-2019R1A6A1A10073079. He also
acknowledges financial support from EPIQS Moore theory centers at MIT and
Harvard. J. K. was supported by the quantum computing technology development
program of the National Research Foundation of Korea(NRF) funded by the Korean
government (Ministry of Science and ICT(MSIT)). (No.2021M3H3A103657312)
Additional support came from National Research Foundation of Korea under Grant
NRF-2021M3E4A1038308 and NRF-2021M3H3A1038085. JHH acknowledges informative
discussion with Ho Tat Lam, Salvatore Pace, and Yizhi You.

##  Appendix A Emergent gauge fields from lattice consideration

Figure 1:  Schematic diagram for the operators (a) $W_{e_{x},m}$ (c) and
$W_{e_{y},m}$ in Eq. (41). The boundary of WW operators are expressed by the
navy lines.

In this section, we elucidate the procedure by which the emergent gauge fields
are obtained from the lattice consideration. Initially the discretized phase
integrals were written as areal products $W_{e_{x},m}$ and $W_{e_{y},m}$:

$\displaystyle W_{e_{x},m}=$
$\displaystyle\prod_{i\in\cal{A}}\left(b_{i}\right)^{y_{x0}-y_{i}},$
$\displaystyle W_{e_{y},m}=$
$\displaystyle\prod_{i\in\cal{A}}\left(b_{i}\right)^{x_{i}-x_{y0}}.$ (41)

By substituting the definition of $b_{i}$ in terms of spin operators, we
obtain equivalent, Wegner-Wilson (WW) loop expressions as a product of various
$X$ operators along the boundary. We consider the case of a rectangular
boundary ${\cal A}=[x_{0},x_{1}]\times[y_{0},y_{1}]$ as depicted in Fig. 1. We
can further decompose them into products of four line operators as

$\displaystyle W_{e_{x},m}$ $\displaystyle=T_{x}(x_{0}\rightarrow
x_{1}|y_{0})T_{x}(y_{0}\rightarrow y_{1}|x_{1})\left[T_{x}(x_{0}\rightarrow
x_{1}|y_{1})\right]^{-1}\left[T_{x}(y_{0}\rightarrow
y_{1}|x_{0})\right]^{-1},$ $\displaystyle W_{e_{y},m}$
$\displaystyle=T_{y}(x_{0}\rightarrow x_{1}|y_{0})T_{y}(y_{0}\rightarrow
y_{1}|x_{1})\left[T_{y}(x_{0}\rightarrow
x_{1}|y_{1})\right]^{-1}\left[T_{y}(y_{0}\rightarrow
y_{1}|x_{0})\right]^{-1}.$ (42)

For instance, $T_{x}(x_{0}\rightarrow x_{1}|y_{0})$ involves the product of
operators at $x=x_{0}$ through $x=x_{1}$ for a fixed $y=y_{0}$. Specifically,

$\displaystyle T_{x}(x_{1}\rightarrow
x_{2}|y)\\!=\\!\prod_{x=x_{1}}^{x_{2}-1}\left(X_{1,i}\right)^{-1}\left(X_{1,i}X_{1,i\\!-\\!\hat{y}}^{-1}\right)^{y-y_{x0}},$
$\displaystyle T_{x}(y_{1}\rightarrow
y_{2}|x)\\!=\\!\prod_{y=y_{1}}^{y_{2}-1}\left(X_{0,i}X_{0,i\\!-\\!\hat{y}}^{-1}X_{2,i}^{-1}X_{2,i\\!-\\!\hat{x}}\right)^{y\\!-\\!y_{x0}},$
(43)

and,

$\displaystyle T_{y}(x_{1}\rightarrow
x_{2}|y)\\!=\\!\prod_{x=x_{1}}^{x_{2}-1}\left(X_{0,i}X_{0,i\\!-\\!\hat{x}}^{-1}X_{1,i}^{-1}X_{1,i\\!-\\!\hat{y}}\right)^{x\\!-\\!x_{y0}},$
$\displaystyle T_{y}(y_{1}\rightarrow
y_{2}|x)\\!=\\!\prod_{y=y_{1}}^{y_{2}-1}\left(X_{2,i}\right)^{-1}\left(X_{2,i}X_{2,i\\!-\\!\hat{x}}^{-1}\right)^{x-x_{y0}}.$
(44)

The shape of the boundary can be relaxed. For any type of boundary made by
putting together several rectangles of arbitrary sizes, one still finds
complete cancellation of operators except those at the boundary, which can be
decomposed into the products of $T$ operators given in Eqs. (43) and (44).

Figure 2: Graphical illustration for (a) $T_{x}(x_{0}\rightarrow x_{1}|y_{0})$
and (b) $T_{x}(y_{0}\rightarrow y_{1}|x_{1})$. The navy line represents the
domain of $T_{x}$. The cyan square represents $e_{x}$ monopole with charge
$\pm 1$ together with other charges comprising the auxiliary dipoles. The
orange square represents the other pair that make up the auxiliary dipole. The
auxiliary dipoles are connected by the gray line.

The first operator in $W_{e_{x},m}$ is $T_{x}(x_{0}\rightarrow x_{1}|y_{0})$,
which is illustrated in Fig. 2(a). It reduces to
$\prod_{x=x_{0}}^{x_{1}-1}\left(X_{1,i}\right)^{-1}$ when $y_{0}-y_{x0}~{}{\rm
mod}~{}N=0$ and creates $e_{x}$ monopoles with charge $-1$ at $(x_{0},y_{0})$
and charge $+1$ at $(x_{1},y_{0})$. In other words, the operator translates
the $e_{x}$ monopole of charge +1 from $(x_{0},y_{0})$ to $(x_{1},y_{0})$. For
$y_{0}-y_{x0}~{}{\rm mod}~{}N\neq 0$, we need to take account of the action by
the auxiliary factor
$\left(X_{1,i}X_{1,i\\!-\\!\hat{y}}^{-1}\right)^{y_{0}-y_{x0}}$ as well, which
is the creation of a pair of $e_{x}$ dipoles. One of the dipoles created
consists of two $e_{x}$ monopoles with charges $(y_{0}-y_{x0})$ at
$(x_{0},y_{0})$ and $-(y_{0}-y_{x0})$ at $(x_{0},y_{0}-1)$. The other dipole
consists of two $e_{x}$ monopoles with charges $-(y_{0}-y_{x0})$ at
$(x_{1},y_{0})$ and $(y_{0}-y_{x0})$ at $(x_{1},y_{0}-1)$.

We next look into the action of $T_{x}(y_{0}\rightarrow y_{1}|x_{1})$
operator, which is illustrated in Fig. 2(b). First of all one can write
$T_{x}(y_{0}\rightarrow y_{1}|x_{1})=T_{x}(y_{0}\rightarrow
y_{x0}|x_{1})T_{x}(y_{x0}\rightarrow y_{1}|x_{1})$. When both
$y_{0}-y_{x0}~{}{\rm mod}~{}N=0$ and $y_{1}-y_{x0}~{}{\rm mod}~{}N=0$, one can
think of it as a product of two translation operators of $e_{x}$, namely
$T_{x}(y_{0}\rightarrow y_{x0}|x_{1})$ and $T_{x}(y_{x0}\rightarrow
y_{1}|x_{1})$, moving $e_{x}$ from $(x_{1},y_{0})$ to $(x_{1},y_{x0})$ and
then subsequently from $(x_{1},y_{x0})$ to $(x_{1},y_{1})$. When
$y_{0}-y_{x0}~{}{\rm mod}~{}N\neq 0$ or $y_{1}-y_{x0}~{}{\rm mod}~{}N\neq 0$,
either the action of $T_{x}(y_{0}\rightarrow y_{x0}|x_{1})$ or
$T_{x}(y_{x0}\rightarrow y_{1}|x_{1})$ is accompanied by the creation of an
auxiliary dipole near $(x_{1},y_{0})$ or $(x_{1},y_{1})$, respectively. We can
understand the action of $T_{x}$ and its inverse $T_{x}^{-1}$ as performing
the translation of $e_{x}$ monopole while either creating or annihilating an
auxiliary dipole. At the end of the $W_{e_{x},m}$ operation, all the auxiliary
dipoles disappear and the $e_{x}$ monopole has completed a loop.

By following the same procedure, one can interpret $T_{y}$ as the translation
operator moving the $e_{y}$ monopole from one point to another and with or
without the accompanying dipole depending on the modality of the coordinates
with respect to $N$. The auxiliary dipoles are aligned in the $x$-direction
and disappear at the completion of the loop.

Figure 3:  Graphical illustrations for the operators (a) $W_{m,e_{x}}$ and (b)
$W_{m,e_{y}}$ in Eq. (45). The boundary of WW operators are expressed by the
navy lines.

The operators $W_{m,e_{x}}$ and $W_{m,e_{y}}$ are defined as

$\displaystyle W_{m,e_{x}}=$ $\displaystyle\prod_{i\in{\cal
A}}\left(a_{i}^{x}\right)^{y_{i}-y_{m0}},$ $\displaystyle W_{m,e_{y}}=$
$\displaystyle\prod_{i\in{\cal A}}\left(a_{i}^{y}\right)^{x_{m0}-x_{i}}.$ (45)

In the case of the rectangular boundary ${\cal
A}=[x_{0},x_{1}]\times[y_{0},y_{1}]$, $W_{m,e_{x}}$ and $W_{m,e_{y}}$ are
expressed in Fig. 3. The presence of $Z_{0}$ operators in the interior of the
boundary makes it impossible to decompose $W_{m,e_{x}}$ and $W_{m,e_{y}}$ into
the product of line operators as was the case in $W_{e_{x},m}$ and
$W_{e_{y},m}$. On the other hand, on inspecting Fig. 3, one realizes that the
$Z_{0}$ operators in the interior cancel out by multiplying $W_{m,e_{x}}$ and
$W_{m,e_{y}}$. Hence, we can decompose $W_{m,e}\equiv W_{m,e_{x}}W_{m,e_{y}}$
as the product of line operators:

$\displaystyle W_{m,e}$ $\displaystyle=T_{m}(x_{0}\rightarrow
x_{1}|y_{0})T_{m}(y_{0}\rightarrow y_{1}|x_{1})\left[T_{m}(x_{0}\rightarrow
x_{1}|y_{1})\right]^{-1}\left[T_{m}(y_{0}\rightarrow
y_{1}|x_{0})\right]^{-1}.$ (46)

Here, the $T_{m}$ operator is defined by

$\displaystyle T_{m}(x_{1}\rightarrow
x_{2}|y)=\prod_{x=x_{1}\\!+\\!1}^{x_{2}}\left(Z_{2,i}\right)^{x-x_{m0}}\left(Z_{0,i-\hat{x}}\right)^{y_{m0}-y},$
$\displaystyle T_{m}(y_{1}\rightarrow
y_{2}|x)=\prod_{y=y_{1}\\!+\\!1}^{y_{2}}\left(Z_{1,i}\right)^{y-y_{m0}}\left(Z_{0,i-\hat{y}}\right)^{x_{m0}-x}.$
(47)

The complete cancellation of $Z_{0}$ takes place for other shapes of the
boundary as well, and $W_{m,e}$ can always be decomposed into the product of
$T_{m}$’s.

Figure 4: Graphical illustration for (a) $T_{m}(x_{0}\rightarrow x_{1}|y_{0})$
and (b) $T_{m}(y_{0}\rightarrow y_{1}|x_{1})$. The navy line represents the
domain of $T_{x}$. The cyan square represents $m$ monopole with charge $\pm 1$
together with other charges composing the auxiliary dipole. The orange square
represents the other pair that make up the auxiliary dipole. The auxiliary
dipoles are connected by the gray line.

The first operator in $W_{m,e}$ is $T_{m}(x_{0}\rightarrow x_{1}|y_{0})$,
which is illustrated in Fig. 4(a). Its role is to translate the $m$ monopoles
from $(x_{0},y_{0})$ to $(x_{1},y_{0})$ and creating the auxiliary dipoles
depending on the modality of the coordinates with respect to $N$. When
$y_{0}-y_{m0}~{}{\rm mod}~{}N\neq 0$, two auxiliary $m$ dipoles aligned in the
$y$-direction are created near $(x_{0},y_{0})$ and $(x_{1},y_{0})$,
respectively. The dipole near the point $(x_{0},y_{0})$ has dipole moment
$y_{0}-y_{m0}$, and the other dipole near the point $(x_{1},y_{0})$ has dipole
moment $-(y_{0}-y_{m0})$. When $x_{0}-x_{m0}~{}{\rm mod}~{}N\neq 0$,
$T_{m}(x_{0}\rightarrow x_{1}|y_{0})$ creates an $x$-directed auxiliary dipole
near $(x_{0},y_{0})$ with the dipole moment $x_{0}-x_{m0}$. When
$x_{1}-x_{m0}~{}{\rm mod}~{}N\neq 0$, $T_{m}(x_{0}\rightarrow x_{1}|y_{0})$
creates an $x$-directed auxiliary dipole near $(x_{1},y_{0})$ with the dipole
moment is $-(x_{1}-x_{m0})$. In contrast to $T_{x}$ or $T_{y}$ creating the
auxiliary dipole only with the $x$\- or $y$-component dipole moment, the
auxiliary dipole created by $T_{m}$ has both $x$\- and $y$-components and it
consists of three monopoles. For example, the auxiliary dipole created by
$T_{m}(x_{0}\rightarrow x_{1}|y_{0})$ near $(x_{1},y_{0})$ consists of three
monopoles with charge $-(y_{0}-y_{m0})$ at $(x_{1},y_{0}+1)$, charge
$-(x_{1}-x_{m0})$ at $(x_{1}+1,y_{0})$, and charge
$(y_{0}-y_{m0})+(x_{1}-x_{m0})$ at $(x_{1},y_{0})$.

The graphical illustration of $T_{m}(y_{0}\rightarrow y_{1}|x_{1})$ is given
in Fig. 4(b). After the application of $T_{m}(x_{0}\rightarrow x_{1}|y_{0})$
to the vacuum, applying $T_{m}(y_{0}\rightarrow y_{1}|x_{1})$ in succession
translates the $m$ monopole with charge $+1$ and the auxiliary dipole from
$(x_{1},y_{0})$ to $(x_{1},y_{1})$. Subsequently, applying
$\left[T_{m}(x_{0}\rightarrow x_{1}|y_{1})\right]^{-1}$ and
$\left[T_{m}(y_{0}\rightarrow y_{1}|x_{0})\right]^{-1}$ moves the $m$
monopoles and the auxiliary dipole from $(x_{1},y_{1})$ to $(x_{0},y_{1})$ to
$(x_{0},y_{0})$, which completes the entire loop. In summary, during the
operation of $W_{m,e}$, the consecutive $T_{m}$’s operation braids the $m$
monopole and the accompanied auxiliary dipole around the boundary of the
$W_{m,e}$ operator.

##  Appendix B Auxiliary dipoles and continuity equations

In this section, we show how to write the charge densities of $e_{x}$,
$e_{y}$, and $m$ quasiparticles in the presence of their auxiliary dipoles in
the continuum theory, and solve the appropriate continuity equations to find
the corresponding current densities that match the charge densities. In Sec. A
we have shown that in the spin system, the translation of $e_{x}$, $e_{y}$,
and $m$ monopoles are accompanied by the auxiliary dipoles. The auxiliary
dipole moment for $e_{x}$ ($e_{y}$) varies in a way that conserves the
$y$-component ($x$-component) of the total dipole moment as the $e_{x}$
($e_{y}$) monopole is translated. The auxiliary dipole for $m$ varies in a way
that conserves both the $x$\- and $y$-component for the $m$ monopole and the
accompanying dipole.

First, consider the $m$ monopole of unit charge moving around by the
applications of $T_{m}$ operators. When the $m$ particle is translated from
the initial point ${\bm{r}}_{m0}=(x_{m0},y_{m0})$ to some point
${\bm{r}}_{m}=(x_{m},y_{m})$, it is accompanied by the auxiliary dipole. The
auxiliary dipole, according to the lattice calculation, consists of three $m$
monopoles with charge $-(x_{m}-x_{m0})$ at $(x_{m}+1,y_{m})$, charge
$(x_{m}-x_{m0})+(y_{m}-y_{m0})$ at $(x_{m},y_{m})$, and charge
$-(y_{m}-y_{m0})$ at $(x_{m},y_{m}+1)$. The combined dipole moment of the
auxiliary dipole is

$\displaystyle{\bm{\mu}}^{m}=-(x_{m}-x_{m0},y_{m}-y_{m0}),$ (48)

which precisely compensates for the dipole moment incurred by the motion of
the $m$ particle from ${\bm{r}}_{m0}$ to ${\bm{r}}_{m}$.

To take the continuum limit we switch the lattice constant from 1 to $\Delta$
and take the limit $\Delta\rightarrow 0$. The positions of three monopoles
composing the auxiliary dipole become $(x_{m}+\Delta,y_{m})$, $(x_{m},y_{m})$,
and $(x_{m},y_{m}+\Delta)$. The monopole charges are modified accordingly to
$-\frac{1}{\Delta}(x_{m}-x_{m0})$,
$\frac{1}{\Delta}(x_{m}-x_{m0})+\frac{1}{\Delta}(y_{m}-y_{m0})$, and
$-\frac{1}{\Delta}(y_{m}-y_{m0})$, respectively. We can write the charge
density for the auxiliary dipole as

$\displaystyle\rho^{{\bm{\mu}}^{m}}({\bm{r}})=$
$\displaystyle\frac{\left(x_{m}\\!-\\!x_{m0}\right)}{\Delta}\Big{[}-\\!\delta\left(x\\!-\\!x_{m}\\!-\\!\Delta\right)\delta\left(y\\!-\\!y_{m}\right)\\!+\\!\delta\left(x\\!-\\!x_{m}\right)\delta\left(y\\!-\\!y_{m}\right)\Big{]}$
$\displaystyle+\\!\frac{\left(y_{m}\\!-\\!y_{m0}\right)}{\Delta}\Big{[}-\\!\delta\left(x\\!-\\!x_{m}\right)\delta\left(y\\!-\\!y_{m}\\!-\\!\Delta\right)\\!+\\!\delta\left(x\\!-\\!x_{m}\right)\delta\left(y\\!-\\!y_{m}\right)\Big{]},$
(49)

which becomes, in the $\Delta\rightarrow 0$ limit,

$\displaystyle\rho^{{\bm{\mu}}^{m}}({\bm{r}})=$
$\displaystyle(x_{m}-x_{m0})\partial_{x}\delta(x-x_{m})\delta(y-y_{m})$
$\displaystyle+(y_{m}-y_{m0})\delta(x-x_{m})\partial_{y}\delta(y-y_{m})$
$\displaystyle=$
$\displaystyle({\bm{r}}_{m}-{\bm{r}}_{m0})\cdot{\bm{\nabla}}\delta^{2}({\bm{r}}-{\bm{r}}_{m}),$
(50)

where $\delta^{2}({\bm{r}}-{\bm{r}}_{m})=\delta(x-x_{m})\delta(y-y_{m})$.
Together with the original $m$ monopole density
$\rho^{m,0}({\bm{r}})=\delta^{2}({\bm{r}}-{\bm{r}}_{m})$, the net density
becomes

$\displaystyle\rho^{m}({\bm{r}})=$
$\displaystyle\rho^{m,0}({\bm{r}})+\rho^{{\bm{\mu}}^{m}}({\bm{r}})$
$\displaystyle=$
$\displaystyle\delta^{2}({\bm{r}}-{\bm{r}}_{m})+({\bm{r}}_{m}-{\bm{r}}_{m0})\cdot{\bm{\nabla}}\delta^{2}({\bm{r}}-{\bm{r}}_{m}).$
(51)

Similar consideration applies to the two electric monopoles. When the $e_{x}$
monopole is translated from ${\bm{r}}_{x0}=(x_{x0},y_{x0})$ to
${\bm{r}}_{x}=(x_{x},y_{x})$, it is accompanied by an auxiliary dipole
consisting of $e_{x}$ monopole with charge $(x_{x}-x_{x0})$ at
$(x_{x},y_{x}-1)$ and charge $-(x_{x}-x_{x0})$ at $(x_{x},y_{x})$. By
following the same procedure for taking the continuum limit as before, we
obtain the charge density of $e_{x}$ particle as

$\displaystyle\rho^{x}=$
$\displaystyle\rho^{x,0}+\rho^{\mu^{x}}=\delta^{2}({\bm{r}}-{\bm{r}}_{x})+(y_{x}-y_{x0})\partial_{y}\delta^{2}({\bm{r}}-{\bm{r}}_{x}),$
$\displaystyle\rho^{y}=$ $\displaystyle 0,$ (52)

where $\rho^{x,0}({\bm{r}})=\delta^{2}({\bm{r}}-{\bm{r}}_{x})$ is the $e_{x}$
monopole density, and $\rho^{\mu^{x}}$ is the density of the auxiliary dipole.
Note that we need to consider the density of the $e_{y}$ monopole
simultaneously, albeit $\rho^{y}=0$, since it is coupled to $\rho^{x}$ by the
continuity equations in Eq. (26). For the motion of $e_{y}$ we find

$\displaystyle\rho^{y}=$
$\displaystyle\rho^{y,0}+\rho^{\mu^{y}}=\delta^{2}({\bm{r}}-{\bm{r}}_{x})+(x_{y}-x_{y0})\partial_{x}\delta^{2}({\bm{r}}-{\bm{r}}_{x}),$
$\displaystyle\rho^{x}=$ $\displaystyle 0.$ (53)

## References

  * Deser  _et al._ (1982) S. Deser, R. Jackiw, and S. Templeton, Annals of Physics 140, 372 (1982).
  * Zhang  _et al._ (1989) S. C. Zhang, T. H. Hansson, and S. Kivelson, Phys. Rev. Lett. 62, 82 (1989).
  * Zhang (1992) S. C. Zhang, International Journal of Modern Physics B 6, 25 (1992).
  * Wilczek (1982a) F. Wilczek, Phys. Rev. Lett. 48, 1144 (1982a).
  * Wilczek (1982b) F. Wilczek, Phys. Rev. Lett. 49, 957 (1982b).
  * Hansson  _et al._ (2004) T. Hansson, V. Oganesyan, and S. L. Sondhi, Annals of Physics 313, 497 (2004).
  * Haah (2011) J. Haah, Phys. Rev. A 83, 042330 (2011).
  * Vijay  _et al._ (2016) S. Vijay, J. Haah, and L. Fu, Phys. Rev. B 94, 235157 (2016).
  * Pretko (2017) M. Pretko, Phys. Rev. B 95, 115139 (2017).
  * Slagle and Kim (2018) K. Slagle and Y. B. Kim, Phys. Rev. B 97, 165106 (2018).
  * Pretko and Radzihovsky (2018) M. Pretko and L. Radzihovsky, Phys. Rev. Lett. 120, 195301 (2018).
  * Nandkishore and Hermele (2019) R. M. Nandkishore and M. Hermele, Annual Review of Condensed Matter Physics 10, 295 (2019).
  * Pretko  _et al._ (2020) M. Pretko, X. Chen, and Y. You, International Journal of Modern Physics A 35, 2030003 (2020).
  * Kim and Han (2021) J. Kim and J. H. Han, Phys. Rev. B 104, 115128 (2021).
  * Slagle and Kim (2017) K. Slagle and Y. B. Kim, Phys. Rev. B 96, 195139 (2017).
  * Pretko (2018) M. Pretko, Phys. Rev. B 98, 115134 (2018).
  * Shirley  _et al._ (2018) W. Shirley, K. Slagle, Z. Wang, and X. Chen, Phys. Rev. X 8, 031051 (2018).
  * You  _et al._ (2020) Y. You, T. Devakul, S. L. Sondhi, and F. J. Burnell, Phys. Rev. Research 2, 023249 (2020).
  * Dubinkin  _et al._ (2020) O. Dubinkin, A. Rasmussen, and T. L. Hughes, Annals of Physics 422, 168297 (2020).
  * Seiberg and Shao (2020) N. Seiberg and S.-H. Shao, SciPost Phys. 9, 46 (2020).
  * Seiberg and Shao (2021a) N. Seiberg and S.-H. Shao, SciPost Phys. 10, 27 (2021a).
  * Seiberg and Shao (2021b) N. Seiberg and S.-H. Shao, SciPost Phys. 10, 3 (2021b).
  * You  _et al._ (2021) Y. You, F. J. Burnell, and T. L. Hughes, Phys. Rev. B 103, 245128 (2021).
  * Ma  _et al._ (2018) H. Ma, M. Hermele, and X. Chen, Phys. Rev. B 98, 035111 (2018).
  * Bulmash and Barkeshli (2018) D. Bulmash and M. Barkeshli, Phys. Rev. B 97, 235112 (2018).
  * Oh  _et al._ (2022) Y.-T. Oh, J. Kim, E.-G. Moon, and J. H. Han, Phys. Rev. B 105, 045128 (2022).
  * Kim  _et al._ (2022) J. Kim, Y.-T. Oh, and J. H. Han, arXiv preprint arXiv:2203.14917 (2022).
  * Gorantla  _et al._ (2022) P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, Phys. Rev. B 106, 045112 (2022).
  * Pace and Wen (2022) S. D. Pace and X.-G. Wen, Phys. Rev. B 106, 045145 (2022).
  * Delfino  _et al._ (2022) G. Delfino, W. B. Fontana, P. R. S. Gomes, and C. Chamon, arXiv preprint arXiv:2207.00409 (2022).
  * Oh et al. (2022) Y. T. Oh et al., to appear (2022).