# A Probabilistic Approach to The Perfect Sum Problem

Kristof Pusztai
Columbia University Department of Statistics
<EMAIL_ADDRESS>

###### Abstract

The subset sum problem is known to be an NP-hard problem in the field of
computer science with the fastest known approach having a run-time complexity
of $O(2^{0.3113n})$. A modified version of this problem is known as the
perfect sum problem and extends the subset sum idea further. This extension
results in additional complexity, making it difficult to compute for a large
input. In this paper, I propose a probabilistic approach which approximates
the solution to the perfect sum problem by approximating the distribution of
potential sums. Since this problem is an extension of the subset sum, our
approximation also grants some probabilistic insight into the solution for the
subset sum problem. We harness distributional approximations to model the
number of subsets which sum to a certain size. These distributional
approximations are formulated in two ways: using bounds to justify normal
approximation, and approximating the empirical distribution via density
estimation. These approximations can be computed in $O(n)$ complexity, and can
increase in accuracy with the size of the input data making it useful for
large-scale combinatorial problems. Code is available at
https://github.com/KristofPusztai/PerfectSum.

_Keywords_ Subset Sum Problem $\cdot$ Perfect Sum Problem $\cdot$ Probability
Theory $\cdot$ Central Limit Theorem $\cdot$ Combinatorics

##  1 Introduction

Combinatorial problems have a notorious reputation for being computationally
hard to solve efficiently and, as a result, understanding the underlying
properties proves useful in many real life scenarios. The difficulty of
solving such problems has made them particularly relevant in the field of
cryptography with ciphers such as the well known ENIGMA machine used by
Germans during World War II being a prime example(Prasad and Kumari, 2020).
Such encryptions make use of of the immense solution space associated with
such combinatorial problems, making cracking them a complex task(Sarvate and
Seberry, 1986).

Additionally, we can model lots of natural events via a combinatorial
approach. Examples where this is especially relevant can be found in
combinatorial biology, specifically, the study of gene re-combination, where
it is important to understand the dynamics of different possible combinations
of gene alterations(Waterman, 1995)(Reményi et al., 2004). Phenotypic
convergence refers to multiple different combinations of genomic events
leading to a similar biological outcome, or trait(Rosenblum et al., 2014). How
many different ways (i.e. through what combinations of individual genomic
events) a particular phenotype may emerge is of interest in disease genetics.
Ultimately, any task which involves finding subgroups of choices out of a set
of total possible choices will be well modelled by problems such as the one we
are examining.

Many past works have focused on finding the most efficient solutions to such
problems or investigating specific algorithm behavior. Despite many
sophisticated approaches to this problem, the fastest deterministic solution
runs in $O(2^{0.3113n})$ complexity(Howgrave-Graham and Joux, 2010). There
exists a faster algorithm that implements some form of randomization and so
inherently relies on probabilistic aspects, but this work does not focus
heavily on the topic and instead is interested in the underlying dynamics of
their devised algorithm(Cao and Liu, 2018). In this work, we will focus
instead on exploring the purely probabilistic aspects which will allow us to
approximate a solution with relative computational ease rather then trying to
find the exact solution.

##  2 Past Works

Recent works in this area are mainly focused on fast deterministic and
randomized approaches. In fact, such scenarios are commonly used to introduce
dynamic programming, specifically, the Bellman method(Bellman, 1961). Other
solution methods include exhaustive search and divide and-conquer approaches
with different modifications yielding sophisticated and efficient algorithms.
One such modification to general exhaustive search by (Cao and Liu, 2018) uses
a novel data arrangement to find a solution with a computation complexity of
$O(2^{n})$ and is the first published algorithm which returns all subsets. As
a result, this method finds the exact solution to the perfect sum. In this
paper, a probabilistic algorithm is also introduced which implements random
permutations and a truncation. However, this algorithm only works well for
large $n$ and is designed to answer the decision problem.

Some more statistically oriented works also exist but are heavily focused on
analyzing specific algorithms. The most notable of these is the work of
(D’Atri and Puech, 1982) which investigates the limiting behavior of variants
of the GOLOSONE algorithm. This work obtains exact distributions for the
output parameters of the algorithms. While the result is certainly
interesting, it does not reveal any direct properties of the problem itself.
Our work focuses directly on the problem itself and does not deal with any
specific solution algorithm.

##  3 Problem Formulation

First we will define the subset sum problem and then the perfect sum problem
to show that the latter is an extension of the former.

###  3.1 Subset Sum

The subset sum problem decision problem is defined formally as: given a set
$S=\\{x_{1},x_{2},...,x_{n}\\}$ and some value $T$, the subset sum problem is
to decide if there is a subset of $S$ that sums up to $T$(Koiliaris and Xu,
2018). That is, there exists some subset $S_{k}\subseteq S$, such that $\sum
S_{k}=\sum_{x\in S_{k}}x=T$

###  3.2 Perfect Sum

The perfect sum can be stated as follows: given a set of values and a target
sum, we want to find the number of all subsets which sum to our target value.
Formally, this can be defined in the following set of conditions.

Given a set of size n, $S=\\{x_{1},x_{2},...,x_{n}\\}$, and some value $T$ we
want to find all such subsets $S_{k}=\\{x_{i},...,x_{j}\\}$, $1\leq i<j\leq
n$, $S_{k}\subseteq S$ so that $\sum S_{k}=\sum_{x\in S_{k}}x=T$ for all
$k=1,...,n$. In this paper, we will focus on finding the _total_ number of all
such subsets instead of the actual subsets themselves.

Clearly this is an extension of subset sum since the solution to the perfect
sum gives us sufficient information to solve the decision problem while also
providing additional information.

##  4 Probabilistic Analysis

Using our definitions from section 3.2, we can treat any subset $S_{k}$ as a
sample taken without replacement from the set S. In order to inspect the
dynamics of the subset further, we can calculate the probability that any
value of the set S, $x\in S$, will be in the subset $S_{k}$. These
calculations are not novel, and their proofs can be found in the appendix for
rigor. The novelty will come in the application of these calculations in
approximating a solution.

####  Lemma 2:

The event that value $x\in S$ is chosen to be in $S_{k}$ has probability
$P(x\in S_{k})=k/n$ where $k$ denotes the cardinality of the subset $S_{k}$,
$k=|S_{k}|$, and $n$ denotes the cardinality of the set $S$, $n=|S|$.

This subtle property is simple yet useful as we can use this to show that the
expected average of any subset will be the same as the average of the original
set.

####  Theorem 1:

The expected value of the sum of any subset of size $k$, $\sum{S_{k}}$, is
equal to $k$ times the mean of the total set, $\bar{S}$.

Formally:

$E[\sum{S_{k}}]=k\bar{S}$

We are now faced with a more difficult task which is to calculate the variance
of the sums of some subset of size k. This is made complex by the fact that we
are sampling without replacement which induces non-zero covariance among each
member of the sample.

From the formula for variance we have:

$Var[\sum{S_{k}}]=Var[\sum_{x\in S_{k}}x]=\sum_{x\in
S_{k}}Var[x]+\sum_{x_{1},x_{2}\in S_{k},x_{1}\neq x_{2}}Cov[x_{1},x_{2}]$

Hence, we know that the variance of the sum of random variables $x\in S_{k}$
will include covariance terms due to the lack of replacement. Intuitively, we
know that this covariance should be negative since if we have two chosen
values from $S_{k}$, call them $x_{1}$ and $x_{2}$ then we know that if
$x_{1}$ has a large value we can no longer pick this large value when choosing
$x_{2}$. Additionally, we know that the total size of the set we are choosing
from should also be present since as the set size approaches infinity, this
limit approaches sampling with replacement. With these in mind, we can take a
look at the formula for covariance and notice a few things.

$Cov[x_{1},x_{2}]=E[x_{1}*x_{2}]-E[x_{1}]*E[x_{2}]\textrm{ \ \ \ for
$x_{1},x_{2}\in S_{k}$ and $x_{1}\neq x_{2}$}$

Firstly, we note that $E[x_{1}]=E[x_{2}]$ for any $x_{1},x_{2}\in S_{k}$ since
all values $x\in S$ have the same probability of being in $S_{k}$ from our
calculations in Theorem 1, and all values in $x\in S_{k}$ have the same
probability, $\frac{1}{k}$, of being chosen. In fact, we can show that
$E[x]=\bar{S}$ where $\bar{S}$ denotes the mean of the original size $n$ set.
A short proof of this is left in the appendix section 9.1. Now we are left to
find the $E[x_{1}*x_{2}]$ term which happens to be something quite convenient
as we see in below in Lemma 2.

####  Lemma 2:

$E[x_{1}*x_{2}]=-\frac{\sigma^{2}}{n-1}+\bar{S}^{2}$

We now have all of our terms in covariance accounted for which allows us to
conclude Lemma 3 by simply plugging in the values from our previous
derivations.

####  Lemma 3:

$\textrm{for any $x_{1},x_{2}\in S_{k}$ we have that
}Cov[x_{1},x_{2}]=\frac{-\sigma^{2}}{n-1}$

We now have all the components to find the variance of the sum of a subset and
by simply plugging in our calculated values from the above sections we arrive
at Theorem 2.

####  Theorem 2:

$Var[\sum{S_{k}}]=k\sigma^{2}(1-\frac{(k-1)}{n-1})$

This conclusion provides us with significant information about the
distribution of possible sum values and is vital to our application for
approximating a solution. Note that as $k$ approaches $n$, the variance goes
to 0.

##  5 Application to Finding Solutions

In the following sections, we provide applications for the above probabilistic
properties in finding a solution. We will explore various distributional
approximation methods and apply them to calculate a concrete value which we
compare to the true answer for solvable randomly generated problems.
Psuedocode for our general algorithm structure can be found in the appendix
and can be seen to have $O(n)$ complexity.

###  5.1 Exact Solution

With the previous calculations in mind, we can try and apply these conclusions
to finding a solution to the perfect sum problem. From a probabilistic
perspective, if we know the exact discrete distribution of the larger set $S$,
then we can calculate the exact distribution of each of the subsets $S_{k}$.
As a result, we can find the percentage of subsets of each size $k$ which sum
up to a given value. In fact, with this distributional perspective we gain
immense flexibility as we can also find the percentage of subsets which sum to
greater, or less then a specified value. Once we have the percentage, we can
multiply this percentage by the total number of different possible subsets of
size $k$, which is just a combination formula computation, and we are left
with the exact number of subsets.

The usefulness of this approach is limited, though, by the fact that we must
calculate exact discrete probability density’s. For example, let’s examine $k$
= 3. Let’s denote the sum of 3 random variables $x_{1}=j,x_{2}=l,x_{3}=p$ and
$X=\sum S_{3}=x_{1}+x_{2}+x_{3}$ then we have:

$P(X=x)=\sum_{j\in S}\sum_{p\in S\backslash
j}P(x_{1}=j)*P(x_{2}=x-j-p|x_{1})*P(x_{3}=x-j-l|x_{1},x_{2})*P(x_{1},x_{2},x_{3}\in
S_{k})$

which leaves us once again, with a combinatorial problem. Ultimately, our
previous results are not useful in the pursuit of an exact solution, but what
they grant us is the ability to calculate accurate approximations easily when
$n$ is large and certain conditions are met. As $n$ gets large, the exact
solution becomes infeasible to calculate but distributional approximation
accuracy increases and is significantly faster to calculate with $O(n)$
complexity. Note that in the above formula, as the size of $k$ increases, the
number of combinations of different possible values taken on by the sum
increases.

###  5.2 Naive Approximation via Normal Distribution

The most naive approximation approach takes form in a modified version of the
central limit theorem (CLT)(Fischer, 2010), allowing us to approximate the
distribution of the sum of any subset as a normal distribution. This
approximation hinges on the assumption that our values in set $S$ are i.i.d
random samples taken from some underlying latent distribution. If these
assumptions do not hold, we cannot invoke the Berry-Esseen bounds shown below
and the accuracy of our approximation will not necessarily increase as set
size and threshold value increase. However, if these assumptions do hold, then
the accuracy of this approximation will increase the larger $T$and $S$ are. We
can see this empirically in Figure 1. For smaller $T$, solutions will contain
smaller subset groupings and can cause inaccuracy. These cases will depend
heavily on the actual distribution of values in S. However, the larger the
size of set $S$, $n$, is and the larger our target, $T$, is, the more subsets
of larger size $k$, which results in a more accurate approximation via this
method as these larger subsets will contain most of the desired sums. More
formally, we can theoretically bound the deviation between the distribution we
are approximating and the normal distribution. Many of the different proofs of
the classical CLT rely on such bounding, examining the absolute difference
between the normal distribution and distribution of the sum of i.i.d random
variables(Fischer, 2010). However, for our case, we are looking at a bound for
the sum of i.i.d samples from a finite population which we assume to be a set
of independent random variables.

With this in mind, we can harness the bound presented in ”Berry–Esseen Bound
for a Sample Sum from a Finite Set of Independent Random Variables”(Zhao et
al., 2004):

$\sup_{x}|P(\frac{S_{k}}{\sqrt{kb}}\leq x)-\Phi(x)|\leq
C\min(\delta_{1},\delta_{2}+1/\sqrt{kq})$

where $C$ is an absolute constant, $\frac{1}{0}:=\infty$,
$p=\frac{k}{n}$,
$q=1-p$,
$b=1-pn^{-1}\sum_{i=1}^{n}(E(x_{i}))^{2}$,
$\delta_{1}=n^{-1}\sum_{i=1}^{n}\frac{E(|x_{i}|^{3})}{\sqrt{(}k)b^{3/2}}$,
$\delta_{2}=n^{-1}\sum_{i=1}^{n}\frac{|E(x_{i})|^{3}}{\sqrt{nb}}+n^{-1}\sum_{i=1}^{n}\frac{E(|x_{i}-pE(x_{i})|^{3})}{\sqrt{n}b^{3/2}}$

This bound gives us a good idea of the theoretical worst case convergence of
our sums to the normal distribution and provides a guarantee that these sums
do, in fact, approach the normal distribution as $n$ increases. However,
calculating this theoretical deviation does not give us much information on
how well the normal approximation is doing with a specific data-set. For a
more concrete and specific result, we can use the discrete Jensen-Shannon
divergence metric(Lin, 1991) over the desired data which quantifies how well
the normal approximation is doing. The discrete Jensen-Shannon divergence
metric over probability space $\Omega$ is defined as:

$JSD(P||Q)=\frac{1}{2}\sum_{x\in\Omega}P(x)log(\frac{P(x)}{M(x)})+\frac{1}{2}\sum_{x\in\Omega}Q(x)log(\frac{Q(x)}{M(x)})$

where $M(x)=\frac{1}{2}(P(x)+Q(x))$

Note that there are many different forms of CLT’s all with their own
assumptions and convergence rates. In fact, there are also general methods for
bounding the difference between distributions. Specifically, Stein’s
method(Stein, 1972) provides a way to calculate bounds of the difference
between normal distribution and any other distribution we might assume our
data comes from. As a result, this naive approximation can still useful and
provide accurate enough approximations with minimal computation. In section
5.4 we will introduce a naive approximation which does not hinge on these
assumptions and instead tries to estimate the density directly via sampling.

(a) JS-Divergence of Normal Approximation: 0.04613
JS-Divergence of Irwin-Hall: 0.03013

(b) JS-Divergence of Normal Approximation: 0.03118
JS-Divergence of Irwin-Hall: 0.02225

(c) JS-Divergence of Normal Approximation: 0.03164
JS-Divergence of Irwin-Hall: 0.02237

(d) JS-Divergence of Normal Approximation: 0.01849
JS-Divergence of Irwin-Hall: 0.01832

Figure 1: In these graphs we see that the normal distribution does well at
approximating the distribution of our subset sums, which come from a uniform
distribution. Notice, however, that the normal distribution is not such a
great approximation for low values of $k$, namely graph (a),$k=1$ and (b)$k=2$
show a larger deviation between the normal approximation and the true values.
Graph (c), $k=3$, and (d), $k=4$, are much better approximated. As the subset
size gets larger, the normal pdf will become a more accurate approximation as
shown by the Berry-Esseen Bounds(Zhao et al., 2004).

###  5.3 Improved Approximation

For scenarios where our set $S$ can be approximated well via some other
continuous distribution besides the normal, we can make our approximation even
more accurate. Formally, if the original set is an i.i.d sample of some known
continuous distribution, then we can use this to find a better distribution of
the sums of each different subset size. This will result in a more accurate
approximation, especially for subsets of smaller size $k$, since this is where
the CLT approximation is the least accurate.

However, we must be careful with this as the normal distribution can still be
a better approximation if the size of $S$ is small. We want the variance of
the sample mean to not be too large to ensure that this alternate distribution
is truly better as even the slightest deviation in sample values can cause
this alternate distribution to fail at accurately describing larger subset
sums. This is due to the fact that smaller samples will have higher variance
in means. We can see this directly through a simple examination of the formula
for calculating sample means:

$\bar{X}=\frac{1}{n}\sum_{i=1}^{n}x_{i}\rightarrow
Var(\bar{X})=\frac{Var(x_{i})}{n}$

However, since the normal approximation takes the mean and standard deviation
of the samples themselves into account, it is able to better correct for these
deviations. We can clearly see this in Figure 2 which demonstrates the effects
of set size on accuracy of using the underlying distribution from which it was
sampled.

(a) JS-Divergence of Normal Approximation: 0.09837
JS-Divergence of Chi-Square: 0.02314

(b) JS-Divergence of Normal Approximation: 0.01370
JS-Divergence of Chi-Square: 0.15208

(c) JS-Divergence of Normal Approximation: 0.09587
JS-Divergence of Chi-Square: 0.01658

(d) JS-Divergence of Normal Approximation: 0.03234
JS-Divergence of Chi-Square: 0.01357

Figure 2: Notice in the simulated set of size 200 that there is a slight
deviation between the theoretical distribution and sample distribution in
picture (a). This causes the approximation to be off for larger subset sums as
seen in (b). In (c) and (d) the set size is much larger, 20,000 samples, which
leads to a better approximation and justifies the use of the alternate
distribution over the normal.

(a) JS-Divergence of Normal Approximation: 0.12712
JS-Divergence of Irwin-Hall: 0.12986

(b) JS-Divergence of Normal Approximation: 0.08121
JS-Divergence of Irwin-Hall: 0.10295

(c) JS-Divergence of Normal Approximation: 0.05145
JS-Divergence of Irwin-Hall: 0.09945

(d) JS-Divergence of Normal Approximation: 0.02421
JS-Divergence of Irwin-Hall: 0.10691

Figure 3: The Normal approximation works better for our discrete example since
it directly takes into account mean and standard deviation of the discrete
set. This discrete set was drawn from a discrete uniform distribution and one
would think trying to approximate via a continuous uniform distribution would
yield relatively good results. However, the Bates distribution (sum of uniform
distributions) does not do a good job at approximation and cannot account for
certain properties resulting from the discreteness of the data.

To quantify when using a distribution other then the normal would perform
better, we can again take an empirical approach and use the discrete Jensen-
Shannon divergence metric(Lin, 1991) to measure how well the two distributions
approximate our subsets. This will then allow us to make exact comparisons
between the two distributions approximation abilities for our specific data.
Additionally, we can take a theoretical approach and note that we want the
mean of our sample set $S$ to be close to the actual mean of the distribution.
Hence, we want some some size $n$ sample which results in a low variance of
$\bar{X}$. The value of $n$ will completely depend on the desired threshold,
and the variance of the underlying distribution from which samples are drawn.

###  5.4 Empirical Approximation Via Density Estimation

For situations where any form of CLT may not be applicable and normal
distribution does not provide a good approximation of the densities, we can
estimate densities directly through sampling and using this sampled empirical
distribution. There are several methods both parametric and non-parametric for
density approximation. Parametric approaches include fitting a Pearson
distribution(Lahcene, 2013), Johnson-$S_{u}$ distribution(Slifker and Shapiro,
1980) or a variety of distributions which can be fitted via either bayesian or
frequentist methods. Non-parametric approaches include the well known method
of kernel density estimation, which is what we will employ in our
implementation due to its flexible nature. Specifically, we will use a tophat
kernel with bandwidth, h, specified by the 10% quantile of the difference
between data points. The estimated probability for some threshold value $t$ is
then calculated as follows:

$P(X=t)=\frac{1}{nh}\sum_{i}^{n}K(\frac{t-x_{i}}{h})$

where $K(x)$ signifies the kernel function with the appropriate properties of
$\int_{-\infty}^{\infty}K(x)\,dx$ = 1 and symmetry, $K(-x)=K(x)$.

We use this to estimate the probability density based on sampling without
replacement from the set $S$ and taking the sums. As a result, we are directly
estimating the distribution of the sum. The choice of Kernel Density
Estimation is due to the flexibility that this estimation method provides, as
it makes no assumptions about the data unlike some other parametric
techniques. However, this flexibility comes at a cost as the number of
samples, choice of bandwidth and choice of kernel will heavily affect the
shape and accuracy of the estimation. Finding the optimal parameters will
likely require extensive exploratory data analysis. Additionally, calculating
such kernel density estimations depends on the number of data points and more
data increases run-time.

##  6 Simulations

There are three cases to test our approximation accuracy on generated sets.
The first case is where the naive approximation does well. We explore both
continuous and discrete cases. We find that in discrete case an improved
approximation method may not exist and that our naive method performs the
best. Specifically, we simulated sets from a discrete scaled uniform
distribution, $U(0,20)$, to create our initial set. In the next section (6.1)
we show that the naive normal approximation here provides better results then
by trying to model this discrete distribution by its continuous counterpart.

However, for cases in which $S$ is drawn from a continuous underlying
distribution, we can use this continuous distribution which is a better
approximation than the normal. In Section 6.2 we will explore the case where
$S$ is drawn from chi-square distribution which is known to be heavily skewed
and so will render our naive normal approximation less accurate. From our
simulations, we can see where our approximations work well, and where they
don’t.

Lastly, in Section 6.3, we will explore cases where independence does not hold
and observe how these effect our normal approximation and kernel density
approximations.

Figure 4: This graph plots the absolute error percentage between ground truth
and approximation values for 20 set simulations. Notice a decreasing trend in
absolute approximation error as set size increases. This is due to the fact
that subsets with larger size begin to dominate the solution and allow for
better approximation. As an increasing number of solutions are found within
larger subsets, our approximation becomes more accurate. Orange deviation
lines signify 1 SD.

(a)

(b)

(c)

(d)

Figure 5: The Normal approximation works better for smaller $n$ due to small
sample variation which the chi-squared approximation is not able to adjust for
despite the samples coming from an underlying chi-square distribution.
However, as $n$ increases, the set follows a chi-squared distribution more
closely and the roles reverse, with the chi-squared approximation becoming
optimal.

###  6.1 Cases for Naive Approximation

For cases where the values of set $S$ are constrained to discrete values, it
may not be feasible to find a more accurate continuous approximation then the
normal distribution. A demonstration of this can be found in appendix section
8.4. Additionally, if no information about the underlying distribution is
known, then the Naive approximation is certainly one approach to modelling the
distribution of sums and, theoretically, will be extremely effective for large
sets and a large target threshold.

In our example, we use a scaled discrete uniform distribution, $U(0,20)$, and
investigate the approximation abilities of the normal distribution. One would
think initially that using a continuous uniform distribution to approximate
this would provide good results, but we must be careful with this. For
example, let’s consider a uniform (0,1) discrete distribution. Note that 0 and
1 both have probability of 0.5 of being picked. Then if we sum two independent
samples from this distribution, the probability of the sum being 2 will be
0.25. If we tried to approximate this via a continuous uniform distribution,
we would estimate the probability of the sum being 2 to be 0, hence not a very
useful approximation.

Due to computational constraints, the ground truth solution for the perfect
sum problem can only be calculated up to set $S$ sizes of 26. However, even
with such a small set size, we see some trends emerge. Specifically, in Figure
4 we see that as the set size increases, our absolute approximation error
decreases.

###  6.2 Cases for Improved approximation

For cases where the values of set $S$ can take on continuous values, we may
find a continuous distribution which approximates these values better then the
normal distribution, but this only works well for a significantly large n.
This is especially the case if the distribution is highly skewed. Another
aspect to consider when approaching the continuous case is that finding the
number of exact sums results in an exact solution of 0. In fact, We get more
insight from looking at how many sets sum up to greater/less then a certain
threshold value. We performed various simulations which explore the validity
of our approximation approach and found promising results. This approximation
does well when the threshold is higher up relative to the set size, since
larger subsets will contribute more to the solution.

Additionally, the normal approximation did better on smaller input set sizes
when compared to the chi-squared approximation. This is likely due to the fact
that the normal approximation directly takes the mean and standard deviation
of the set as parameters making it highly flexible when faced with sample
deviations which occur in small samples sizes. However, as the set size
increased, the chi-squared approximation performed better, with lower
approximation error. This is depicted in Figure 6. Note, however, that despite
despite the normal approximation becoming less accurate for larger n, it is
still a relatively good approximation, with an maximum absolute error of 0.08,
while the chi squared approximation was off by 0.4 in the worst case. The
normal approximation is versatile as it’s parameters allow for sample
deviation adjustments.

###  6.3 Cases for Empirical Approach

For cases where no form of the CLT can be invoked and we have no a priori
knowledge of the distribution, we are faced with two options: continue using
the normal, or estimate the distribution empirically. There are several
benefits to using the normal approximation including that it will not depend
on the accuracy/validity of taking a sample to estimate the empirical
distribution. As a result, no randomness will be involved and this approach
will always give back the same answer for a given set of inputs. Additionally,
the existence of bounds on the normal approximation error for dependent
samples has been shown in Stein(Stein, 1972). However, if no bounds can be
shown and the normal will not be a good approximation then an empirical
approach might best suit the situation. Specifically in situations where the
data is very spread out as seen in Figure 7 and our empirical method
significantly outperforms the naive normal approximation.

(a) JS-Divergence of Normal Approximation: 0.79811
JS-Divergence of KDE: 0.27923

(b) JS-Divergence of Normal Approximation: 0.82409
JS-Divergence of KDE: 0.10890

(c) JS-Divergence of Normal Approximation: 0.79985
JS-Divergence of KDE: 0.13563

(d) JS-Divergence of Normal Approximation: 0.78635
JS-Divergence of KDE: 0.23210

Figure 6: Note the KDE’s ability to account for large gaps between the
possible sums such as cases like a.) and d.). We see that the KDE provides a
better estimation of such distributions then the normal as shown by the JS-
Divergences.

##  7 Applications To Other Problems

This solution approach can be applied to any problem which can be boiled down
to a perfect sum-type question. These questions include the classic ”knapsack”
and ”change-making” problems(Goodrich and Tamassia, 2015)(Silvano, 1990).
However, we must bear in mind that, as mentioned in the problem statements
section, our approach is focused on determining the _total_ number of possible
solutions, and not actually finding these solutions themselves. This leads to
the major drawback of this method in that, while we may have an accurate
approximation of this total value, we are not any closer to knowing what these
sets actually are. In fact, there are no truly efficient algorithms to find
these exact solution sets and so in big-data settings, finding such sets
becomes infeasible, not only computation-wise, but also memory wise as the
number of solutions can grow to hundreds of millions very rapidly(i.e. 100
choose 5 = 75,287,520 different possible combinations). With this in mind, one
can apply our approach to help narrow down the computational scale, focusing
on desired subset-sizes which are likely to contain solutions. For example, to
find the minimum subset size which sums to a desired number, our approach can
be used to see whether a subset of size $k$ is likely to sum up to this
desired value. From here, we can either choose the minimum set that’s
predicted to have at least one subset of this size sum to the desired value,
or we can choose the minimum set that has a probability of such a sum
occurring based on the estimated distribution. If we want to find this actual
solution set, then we can pair this information up with an exact solution
finding algorithm that actually examines the sets to find a solution and focus
on on just subsets of the sizes we’ve identified from the previous step. In
this way, we’ve reduced the computation time significantly as we will only
need to focus on a limited range of different subset sizes instead of all
possible sizes.

Examining real world applications of this approach, it is immediate that this
is particularly relevant to the field of genomics research due to the
combinatorial nature of gene inheritance/recombination. There is a whole sub-
field of biology which focuses on combinatorial methods for gene interaction(J
and S, 2001). Research such as ”Combinatorial control of gene expression”
demonstrates the relevance of such ideas in understanding genomics(Reményi et
al., 2004). Specifically, the ideas presented in this paper are useful if one
wants to find the number of different possible genes which could result in a
specific event occurring. We need only to assign genes a numerical value in
their contribution to this event, and then we can apply our approach to find
the total number/percentage of potentially threshold surpassing combinations.
This could be especially useful for cancer research, where we can think about
each gene combination as contributing a certain amount, and once a stability
level is crossed, a cell turns cancerous. These applications I will leave for
further thought and exploration to the reader.

##  8 Conclusion

In this paper, we introduce a novel approximation approach to the perfect sum
problem which is focused on finding the total number of sub-sets which sum
to/surpass a desired threshold. This approach harnesses probability theory to
estimate the distribution of the sums of subsets. Specifically, we examine the
use of the normal distribution as well as non-parametric density estimation to
estimate the distribution of the sums of different sized subsets when no a
priori distribution information is known. Our algorithm, presented in appendix
section Perfect Sum Psuedocode runs in $O(n)$ complexity and can increase in
accuracy as set size increases, depending on the threshold level and problem
set up. As a result, it directly addresses the need for good approximation in
big data cases where finding the exact solution is infeasible.

##  9 Appendix

###  9.1 Proof Of Expected Value For Members of Subset Size k

for any $x\in S_{K}$ we have:

$E[x]=\sum_{x\in S}x*P(x)$

Note that $P(x)$ denotes the probability of picking the value x from all
possible values in subset $S_{k}$ which is the joint event of this value
ending up in $S_{k}$ first and then being picked out of all other k-1 values

$=\sum_{x\in S}x*1/k*k/n$ $=\sum_{x\in S}x*1/n$ $=1/n\sum_{x\in S}x=\bar{S}$

∎

###  9.2 Proof of Lemma 1

Notice that in a subset of size $k$, for some value $x\in S$ there are
$n-1\choose k-1$ different possible subsets $S_{k}$ which satisfy $x\in
S_{k}$. There are a total of $n\choose k$ subsets of size k. Thus, the
probability of choosing a subset which contains $x$ out of all possible
subsets of size $k$ is:

$P(x\in S_{k})=\frac{{n-1\choose k-1}}{{n\choose
k}}=\frac{\frac{(n-1)!}{(k-1)!(n-k)!})}{\frac{n!}{k!(n-k)!}}=\frac{k!(n-1)!}{(k-1)!n!}=k/n$

∎

###  9.3 Proof of Theorem 1

$\sum{S_{k}}=\sum_{x\in S_{k}}x=\sum_{x\in S}x*\mathbbm{1}\\{x\in
S_{k}\\}\textrm{ \ \ \ where }\mathbbm{1}\\{x\in
S_{k}\\}=\left\\{\begin{array}[]{l}1\textrm{ if }x\in S_{k}\\\ 0\textrm{
otherwise}\end{array}\right.$

Then we have that:

$E[\sum{S_{k}}]=E[\sum_{x\in S}x*\mathbbm{1}\\{x\in S_{k}\\}]$ $=\sum_{x\in
S}x*E[\mathbbm{1}\\{x\in S_{k}\\}]$ $=\sum_{x\in S}x*P(x_{i}\in S_{k})$
$=\sum_{x\in S}x*k/n$ $=\frac{k}{n}\sum_{x\in S}x=k\bar{S}$

∎

###  9.4 Proof of Lemma 2

$E[x_{1}*x_{2}]=\sum_{x_{1},x_{2}\in S,x_{1}\neq
x_{2}}x_{1}*x_{2}*P(x_{1},x_{2}\in S_{k})*P(x_{1},x_{2})$

Following the notation from Lemma 1, $P(x_{1},x_{2}\in S_{k})$ denotes the
probability of the event that values $x_{1}$ and $x_{2}$ get chosen to be in
the subset $S_{k}$. $P(x_{1},x_{2})$ denotes the probability of the event that
values $x_{1}$ and $x_{2}$ get chosen from the values in the subset $S_{k}$.
We need these two events because in order for the ultimate event of $x_{1}$
and $x_{2}$ getting chosen from $S_{k}$, we need to ensure that they actually
are in $S_{k}$ to begin with, otherwise they cannot be chosen so we need to
capture the joint event. Expanding on these probabilities further:

$P(x_{1},x_{2}\in S_{k})=\frac{k}{n}*\frac{k-1}{n-1}$
$P(x_{1},x_{2})=\frac{1}{k(k-1)}$

Hence,

$\sum_{x_{1},x_{2}\in S,x_{1}\neq x_{2}}x_{1}*x_{2}*P(x_{1},x_{2}\in
S_{k})*P(x_{1},x_{2})=\sum_{x_{1},x_{2}\in S,x_{1}\neq
x_{2}}x_{1}*x_{2}*\frac{1}{n(n-1)}$ $=\frac{1}{n(n-1)}(\sum_{x_{1}\in
S}x_{1}*\sum_{x_{2}\in S,x_{2}\neq x_{1}}x_{2})$
$=\frac{1}{n(n-1)}(\sum_{x_{1}\in S}x_{1}*(\sum_{x_{2}\in S}x_{2}-x_{1}))$
$=\frac{1}{n(n-1)}(\sum_{x_{1}\in S}x_{1}*(\sum_{x_{2}\in S}x_{2}-x_{1}))$
$=\frac{1}{n(n-1)}(\sum_{x_{1}\in S}x_{1}*\sum_{x_{2}\in
S}x_{2}-\sum_{x_{1}\in S}x_{1}^{2})\textrm{ \ \ \ note that \ $\sum_{x_{1}\in
S}x_{1}*\sum_{x_{2}\in S}x_{2}=nE[x_{1}]*nE[x_{2}]$}$
$=\frac{1}{(n-1)}(n^{2}E[x_{1}]E[x_{2}]-E[x_{1}^{2}])$
$=\frac{1}{(n-1)}(n\bar{S}^{2}-\sigma^{2}-\bar{S}^{2})$
$=\frac{n}{(n-1)}\bar{S}^{2}-\frac{1}{(n-1)}\sigma^{2}-\frac{1}{(n-1)}\bar{S}^{2}$
$=\bar{S}^{2}-\frac{\sigma^{2}}{(n-1)}$

∎

###  9.5 Proof of Lemma 3

Plugging everything into our formula for covariance we have:

$Cov[x_{1},x_{2}]=E[x_{1}*x_{2}]-E[x_{1}]*E[x_{2}]\textrm{ \ \ \ for
$x_{1},x_{2}\in S_{k}$ and $x_{1}\neq x_{2}$}$
$=-\frac{\sigma^{2}}{n-1}+\bar{S}^{2}-\bar{S}^{2}$ $=-\frac{\sigma^{2}}{n-1}$

∎

###  9.6 Proof of Theorem 2

$Var[\sum_{x\in S_{k}}x]=\sum_{x_{1},x_{2}\in S_{k}}Cov[x_{1},x_{2}]$
$=\sum_{x\in S_{k}}Var[x]+\sum_{x_{1},x_{2}\in S_{k},x_{1}\neq
x_{2}}Cov[x_{1},x_{2}]$
$=k\sigma^{2}+k(k-1)\frac{-\sigma^{2}}{n-1}=k\sigma^{2}(1-\frac{(k-1)}{n-1})$

∎

###  9.7 Perfect Sum Psuedocode

Algorithm 1 Perfect Sum Probabilistic Approximation

1:Compute mean, $\mu$, and variance, $\sigma^{2}$, of distribution(or
approximate distribution) of $S$

2:

3:Create variable, $R=0$, to keep track of total number of subsets which
satisfy sum conditions regarding $T$

4:for $k=1,2,\ldots\textrm{length}(S)$ do

5: Compute values (i.e. $\hat{\sigma^{2}}$, $\hat{\mu}$ via Theorem 1, Theorem
2 ) necessary for distribution approximation choice

6:

7: Compute distribution approximation of sum of subsets size $k$

8:

9: Use distribution to find desired probability, $P=P(\sum
S_{k}\\{=,\geq,\leq\\}T)$, regarding given sum threshold, $T$

10:

11: Convert probability to number of subsets, $s=integer(P*{n\choose k})$

12:

13: Add to variable $R$ to keep track of total over all different sized
subsets, $R=R+s$

14:end for

15:return R

## References

  * Prasad and Kumari (2020) Kalika Prasad and Munesh Kumari.  A review on mathematical strength and analysis of enigma.  Technical report, Department of Mathematics, Central University of Jharkhand, India, 2020.  URL https://arxiv.org/pdf/2004.09982.pdf.
  * Sarvate and Seberry (1986) Dinesh G. Sarvate and Jennifer Seberry.  Encryption methods based on combinatorial designs.  Technical report, University of Wollongong, 1986.  URL https://ro.uow.edu.au/infopapers/1019.
  * Waterman (1995) Michael S. Waterman.  _Handbook of Combinatorics_ , chapter 39.  Elsevier Science B.V., 1995.
  * Reményi et al. (2004) Attila Reményi, Hans R Schöler, and Matthias Wilmanns.  Combinatorial control of gene expression.  _Nature Structural and Molecular Biology_ , 812–815, 2004.
  * Rosenblum et al. (2014) Erica Bree Rosenblum, Christine E. Parent, and Erin E. Brandt.  The molecular basis of phenotypic convergence.  _Annual Review of Ecology, Evolution, and Systematics_ , 45:203-226, 2014.
  * Howgrave-Graham and Joux (2010) Nick Howgrave-Graham and Antoine Joux.  New generic algorithms for hard knapsacks.  _Advances in Cryptology – EUROCRYPT 2010, vol 6110_ , pages 235–256, 2010.
  * Cao and Liu (2018) Zhengjun Cao and Lihua Liu.  New algorithms for subset sum problem.  _CoRR_ , abs/1807.02611, 2018.
  * Bellman (1961) Richard Bellman.  On the approximation of curves by line segments using dynamic programming.  _Communications of the ACM_ , 4 No. 6, 1961.
  * D’Atri and Puech (1982) Gianfranco D’Atri and Claude Puech.  Probabilistic analysis of the subset-sum problem.  _Discrete Applied Mathematics_ , Volume 4, Issue 4, 1982.
  * Koiliaris and Xu (2018) Konstantinos Koiliaris and Chao Xu.  Subset sum made simple, 2018.
  * Fischer (2010) Hans Fischer.  _A History of the Central Limit Theorem_.  Springer New York Dordrecht Heidelberg London, 2010.
  * Zhao et al. (2004) L. C. Zhao, C. Q. Wu, and Qiying Wang.  Berry–esseen bound for a sample sum from a finite set of independent random variables.  _Journal of Theoretical Probability_ , 2004.
  * Lin (1991) Jianhua Lin.  Divergence measures based on the shannon entropy.  _IEEE TRANSACTIONS ON INFORMATION THEORY_ , 37, NO.1, 1991.
  * Stein (1972) Charles Stein.  A bound for the error in the normal approximation to the distribution of a sum of dependent random variables.  _Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability_ , 583-602, 1972.
  * Lahcene (2013) Bachioua Lahcene.  On pearson families of distributions and its applications.  _African Journal of Mathematics and Computer Science Research_ , 6(5), 2013.
  * Slifker and Shapiro (1980) James F. Slifker and Samuel S. Shapiro.  The johnson system: Selection and parameter estimation.  _Technometrics_ , Vol. 22, No. 2, 1980.
  * Goodrich and Tamassia (2015) Michael T Goodrich and Roberto Tamassia.  _Algorithm Design and Applications_.  Wiley, 2015.
  * Silvano (1990) Martello Silvano.  _Knapsack problems : algorithms and computer implementations_.  Wiley, 1990.
  * J and S (2001) Pelletier J and Sidhu S.  Mapping protein-protein interactions with combinatorial biology methods.  _Curr. Opin. Biotechnol._ , 340–7, 2001.
  * Ferrari (2019) Alfredo Ferrari.  A note on sum and difference of correlated chi-squared variables.  _arXiv_ , 1906.09982, 2019.


# Entanglement entropy in the Ising model with topological defects

Ananda Roy<EMAIL_ADDRESS>Department of Physics and
Astronomy, Rutgers University, Piscataway, NJ 08854-8019 USA  Hubert Saleur
Institut de Physique Théorique, Paris Saclay University, CEA, CNRS, F-91191
Gif-sur-Yvette  Department of Physics and Astronomy, University of Southern
California, Los Angeles, CA 90089-0484, USA

###### Abstract

Entanglement entropy (EE) contains signatures of many universal properties of
conformal field theories (CFTs), especially in the presence of boundaries or
defects. In particular, topological defects are interesting since they reflect
internal symmetries of the CFT, and have been extensively analyzed with field-
theoretic techniques with striking predictions. So far, however, no lattice
computation of EE has been available. Here, we present an ab-initio analysis
of EE for the Ising model in the presence of a topological defect. While the
behavior of the EE depends, as expected, on the geometric arrangement of the
subsystem with respect to the defect, we find that zero-energy modes give rise
to crucial finite-size corrections. Importantly, contrary to the field-theory
predictions, the universal subleading term in the EE when the defect lies at
the edge of the subsystem arises entirely due to these zero-energy modes and
is not directly related to the modular S-matrix of the Ising CFT.

Entanglement plays a central role in the development of long-range
correlations in quantum critical phenomena. Thus, quantification of the
entanglement in a quantum-critical system provides a way to characterize the
universal properties of the critical point. The von-Neumann entropy is a
natural candidate to perform this task. For zero-temperature ground-states of
1+1D quantum-critical systems described by conformal field theories (CFTs),
the von-Neumann entropy [equivalently, entanglement entropy (EE)] for a
subsystem exhibits universal logarithmic scaling with the subsystem size [1,
2]. The coefficient of this scaling determines a fundamental property of the
bulk CFT: the central charge. For finite systems with boundaries at a
conformal critical point, the EE receives universal, subleading, boundary-
dependent contributions, the so-called ‘boundary entropy’ [3] – a central
concept in a variety of physical problems both in condensed matter physics and
string theory. This boundary contribution to the EE provides a valuable
diagnostic for identifying the different boundary fixed points of a given CFT
[2, 4, 5, 6, 7].

While entanglement has been analyzed extensively in CFTs with and without
boundaries, its behavior is much less understood in the presence of defects.
The question is particularly intriguing since entanglement measures may
provide an alternate way to classify defects in CFTs. Of particular interest
are topological (perfectly-transmissive) defects [8, 9, 10, 11, 12]. These
defects commute with the generators of conformal transformations and thus, can
be deformed without affecting the values of the correlation functions as long
as they are not taken across field insertions (hence the moniker topological).
They reflect the internal symmetries of the CFT and relate the order-disorder
dualities of the CFT to the high-low temperature dualities of the
corresponding off-critical model [10, 13, 14]. They also play an important
role in the study of anyonic chains and in the correspondence between CFTs and
three-dimensional topological field theories [15].

It is natural to analyze EE in the presence of topological defects. Two
distinct geometries have been considered in the literature: i) where the
defect is entirely within the subsystem, and ii) where the defect is located
precisely at the interface between the subsystem and the rest. While both
cases exhibit identical leading-order logarithmic scaling with subsystem size,
they differ when subleading, i.e., $O(1)$, corrections are taken into account.
In the first case, after usual folding maneuvers, the subleading term can be
equated to a boundary entropy with double the bulk degrees of freedom [16, 17,
18]. In this way, the subleading correction to symmetric EE can be computed
analytically for all rational CFTs [19]. In the second case, the subleading
term for the interface EE is much more difficult to obtain. Field-theory
computations based on the replica trick for the free, real boson [20], the
free, real fermion [21] (see also [22, 23, 24]) and generalization to all
rational CFTs [19] using the corresponding twisted torus partition functions
[8] provide results whose validity have never been tested with ab-initio,
lattice computations. Such a computation is particularly important since the
mapping to the twisted torus partition function does not faithfully capture
the geometric arrangement of the subsystem with respect to the defect. This is
important for subleading terms in the interface EE which are, in fact, the
signatures of the topological defect.

Our goal, in this paper, is to investigate the question in the simplest
possible case, the Ising model, using ab-initio calculations. These
calculations are more complex than is usually the case for EE due to the
presence of zero-energy modes. The latter are a salient feature of the
topological nature of the defect. While the effects of zero-modes have been
extensively quantified in gapped systems due to their relevance in
topologically ordered systems, the same is much less understood in the
critical systems. Here, we show that these zero-modes give rise to nontrivial
contributions to the EE when the subsystem size is comparable to the size of
the whole system, similar to the case of a periodic ring of free fermions [25,
26]. Thus, these zero-modes profoundly affect corrections to scaling. After
the role of these zero-modes is properly taken into account, we find results
that are as expected for case (i), but are not compatible with the formulas
proposed in the string/field theory literature [21, 19] for case (ii).

Figure 1:  Schematic of a single defect in the critical Ising chain with
periodic boundary conditions in the spin (left) and fermionic (right) picture.
Away from the defect, the Hamiltonian contains nearest-neighbor ferromagnetic
$\sigma_{i}^{x}\sigma_{i+1}^{x}$ interaction (black lines) and onsite
transverse field $\sigma_{i}^{z}$ (blue lines). The defect (purple box)
concerns the spins located at sites $L$ and 1. The energy (duality) defect
part of the Hamiltonian is: $b_{\epsilon}\sigma_{L}^{x}\sigma_{1}^{x}/2$
$+\sigma_{L}^{z}/2$ ($b_{\sigma}\sigma_{L}^{y}\sigma_{1}^{x}/2$). After
Jordan-Wigner (JW) transformation, the ferromagnetic coupling corresponds to
an interaction of the form $\gamma_{2i}\gamma_{2i+1}$, while the transverse
field corresponds to $\gamma_{2i-1}\gamma_{2i}$. For the duality defect, there
is a zero-mode localized at $\gamma_{2L}$. Another zero-mode is delocalized
throughout the system (see below).

The Ising model provides a particularly appropriate testbed for the EE in the
presence of topological defects due to the following. First, the lattice
Hamiltonians are well-understood [27, 16]. Second, the zero-mode structure is
the simplest and yet, sufficient to explain the main concepts. Third, the
defect Hamiltonian can be mapped, after Jordan-Wigner (JW) transformation, to
a bilinear fermionic Hamiltonian. The latter can be diagonalized semi-
analytically and leads to very accurate predictions of EE [28, 29, 30]. To
emphasize the nontrivial nature of the topological defect, we compare our
results with the non-topological defects in the model.

There are two nontrivial classes of defects in the Ising CFT: energy
($\epsilon$) and duality ($\sigma$). We consider a periodic Ising chain with
the defect concerning the spins at sites $L$ and 1 [31, 32, 27]. The
Hamiltonians are: $H_{\epsilon,\sigma}=H_{0}$$-h_{\epsilon,\sigma}$ (see Fig.
1). Here, the bulk Hamiltonian term,
$H_{0},=-\sum_{i=1}^{L-1}\sigma_{i}^{x}\sigma_{i+1}^{x}/2$
$-\sum_{i=1}^{L-1}\sigma_{i}^{z}/2$. The two defect terms are:
$h_{\epsilon}=b_{\epsilon}\sigma_{L}^{x}\sigma_{1}^{x}/2$ $+\sigma_{L}^{z}/2$
and $h_{\sigma}=b_{\sigma}\sigma_{L}^{y}\sigma_{1}^{x}/2$ 111Note that
$h_{\epsilon,\sigma}$ are marginal perturbations which, in general, lead to
continuously varying critical exponents [45].. The energy defect consists of
one bond with altered strength $b_{\epsilon}$ connecting spins at sites $L$
and 1. In particular, $b_{\epsilon}=0,+1$ and $-1$ correspond to Ising models
with open, periodic and antiperiodic boundary conditions (bcs) respectively.
Unlike the energy defect, the duality defect consists of a
$\sigma_{L}^{y}\sigma_{1}^{x}$ interaction 222This duality defect Hamiltonian
is related by a local unitary rotation on the $L^{\rm th}$ spin to the one
considered in Refs. [16], which has $\sigma_{L}^{z}\sigma_{1}^{x}$
interaction. We do not use this alternate form since it no longer leads to a
bilinear Hamiltonian under JW transformation and cannot be solved by free-
fermion techniques. . Equally important, there is no transverse field at the
$L^{\rm th}$ site. The duality defect for $b_{\sigma}=1$ (equivalently
$b_{\sigma}=-1$, which is related by a local unitary rotation) is the
topological defect for the Ising CFT.

Next, we perform a JW transformation:
$\gamma_{2k-1}=\sigma_{k}^{x}\prod_{j=1}^{k-1}\sigma_{j}^{z}$,
$\gamma_{2k}=\sigma_{k}^{y}\prod_{j=1}^{k-1}\sigma_{j}^{z}$, where
$\gamma_{j}$-s are real, Majorana fermion operators obeying
$\\{\gamma_{j},\gamma_{k}\\}=2\delta_{j,k}$. In the fermionic language, the
defect Hamiltonians are
$H_{\epsilon,\sigma}^{f}=H_{0}^{f}-h_{\epsilon,\sigma}^{f}$, where

$\displaystyle H_{0}^{f}$
$\displaystyle=\frac{i}{2}\sum_{j=1}^{L-1}\gamma_{2j}\gamma_{2j+1}+\frac{i}{2}\sum_{j=1}^{L-1}\gamma_{2j-1}\gamma_{2j},$
(1) $\displaystyle h_{\epsilon}^{f}$
$\displaystyle=\frac{ib_{\epsilon}}{2}\gamma_{2L}\gamma_{1}-\frac{i}{2}\gamma_{2L-1}\gamma_{2L},\
h_{\sigma}^{f}=-\frac{ib_{\sigma}}{2}\gamma_{2L-1}\gamma_{1}.$ (2)

Here we have restricted ourselves to the symmetry sector
$Q=\prod_{j=1}^{L}\sigma_{j}^{z}=1$ 333Equivalently, we analyze the mixed
sector Hamiltonian [27]:
$H_{\epsilon,\sigma}^{m}=P_{+}H_{\epsilon,\sigma}^{f}(Q=+1)+P_{-}H_{\epsilon,\sigma}^{f}(Q=-1)$.
For the energy defect, from the definition of $H_{\epsilon}^{f}$, we recover
the well-known fact that the periodic (antiperiodic) coupling at the boundary
for the spin model corresponds to antiperiodic (periodic) coupling for the
fermionic model. In particular, the periodic fermionic model
($b_{\epsilon}=-1$) contains two nonlocal Majorana zero-modes, which together
are responsible for the two-fold degenerate ground-state of the fermionic
model. For the duality defect, the operator $\gamma_{2L}$ does not occur in
$H_{\sigma}^{f}$. It commutes with the Hamiltonian:
$[\gamma_{2L},H_{\sigma}^{f}]=0$, and anticommutes with the conserved
$\mathbb{Z}_{2}$ charge: $\\{\gamma_{2L},Q\\}=0$. Thus, it is a zero-mode of
the model which is perfectly localized in space. It has a partner zero-mode
which is completely delocalized:
$\Lambda(b_{\sigma})=\sum_{k=1}^{L}\gamma_{2k-1}+b_{\sigma}\sum_{k=1}^{L-1}\gamma_{2k}$.
Note that the zero-modes exist for all values of $b_{\sigma}$ and are not
special features of the topological point. The fermionic Hamiltonian also
reaffirms a CFT result [27]: $H_{\sigma}^{f}$ describes a chain of $2L-1$
Majorana fermions or equivalently, $L-1/2$ spins. This is important for
quantifying finite-size effects.

Now, we compute symmetric and interface EEs for the ground states of the
fermionic Hamiltonians $H_{\epsilon,\sigma}^{f}$. Since the latter are
bilinear in the fermionic operators, we compute the relevant EEs from the
ground-state correlation matrix [28, 29, 30]. The latter is calculated from
the ground state by filling up the negative energy states. The method is
unambiguous in the absence of zero-energy modes. However, in the presence of
the latter (e.g., $b_{\epsilon}=-1$ or any $b_{\sigma}$), it raises the
question: are the zero-energy states empty or occupied in the ground state?
Yet another possibility is to consider an incoherent superposition of filled
and empty states. This leads to the total system being in a mixed state, but
is appropriate when taking the zero-temperature limit of a thermal ensemble
[25]. The question is crucial to the computation since zero-energy modes
nontrivially affect the EE. For $b_{\epsilon}=-1$, the zero-modes give rise to
nontrivial corrections to the EE of a subsystem of size $r$ within a total
system of size $L$ [25, 26]. The correction $\Delta S(r/L)=$ $\frac{\pi
r}{L}\int_{0}^{\infty}\tanh(\pi rh/L)[\coth(\pi h)-1]$. For $r\ll L$, the EE
is oblivous to the existence of the two nonlocal zero-modes spread throughout
the system: $\Delta S\simeq\pi^{2}r^{2}/12L^{2}\rightarrow 0$. The situation
changes as the subsystem occupies appreciable fraction of the total system
($r\sim L$) culminating in $\Delta S(r=L)=\ln 2$, the latter being the entropy
of the two-fold degenerate ground state of a periodic chain of fermions.
Below, we present analogous results for the topological defect considering
separately the cases when the total system is in a pure and mixed state for
the relevant fermionic models 444See Supplementary Material for details.

Figure 2:  Results for the symmetric EE ($S_{s}$) for a periodic chain of size
$L=3000$ with a single defect and $10<r<100$ (left panel) and $100<r<L/2$
(right panels). The blue crosses are obtained for $b_{\epsilon}=1$ (i.e., no
defect). The maroon (green) hexagons are for the energy defect
$b_{\epsilon}=-1$ (duality defect $b_{\sigma}=1$), when the total system is in
a mixed state. The differences with $S_{s}(b_{\epsilon}=1)$ are $\Delta
S(r/L)$ and $\Delta S(r/L)/2+(\ln 2)/2$ for the two cases. The corresponding
predictions are denoted by maroon and green pluses. Compared to
$b_{\epsilon}=-1$, the EE for $b_{\sigma}=1$ has an additional offset $(\ln
2)/2$ even for $r/L\rightarrow 0$ due to the zero-mode $\gamma_{2L}$ localized
at the center of the interval. As $r/L$ increases, the EE for $b_{\sigma}=1$,
due to the single nonlocal zero mode, $\Lambda(b_{\sigma}=1)$, receives a
contribution half as large as that for $b_{\epsilon}=-1$, which has two such
modes. The maroon squares (green diamonds) present the results when the total
system is in a pure state (the zero-energy state being filled or empty) as
opposed to a mixed state. The results for $b_{\epsilon}=\pm 1$ are
indistinguishable. However, relative to the energy defect, the $b_{\sigma}=1$
shows an offset $\Delta S(1-r/L)/2$, which again reduces to $(\ln 2)/2$ for
$r/L\ll 1$.  Figure 3:  (Left panel) Results for the interface EE ($S_{\cal
I}$) for a periodic chain of $L=500$ with a single defect. One end of the
subsystem is at the defect and the other end sweeps through the system. The
$b_{\epsilon}=-1$ case, both when the system is in an incoherent superposition
(maroon hexagons) and a pure state (maroon squares), can be directly
understood from the symmetric EE (Fig. 2). The topological defect results are
identical for mixed (green hexagons) and pure (green diamonds) total-system
state. Then, the interface EE is best compared to that of a spin-chain of size
$L-1/2$ without defects [see Fig. 1 and discussion below Eq. (1)]. The
difference between the two EEs is $\Delta S[r/(L-1/2)]+\delta_{r,L}(\ln 2)/2$
(orange pluses). For comparison, the corresponding predictions obtained by
comparing to size-$L$ chain without defect (green pluses) are also plotted.
For $r$ not too close to $L$, the difference is negligible. But, zooming in
around $r\sim L$ (inset) makes the difference manifest. (Right panel) We show
the scaling of the interface EE for $r=L/2$ as a function of $\ln L$. Every
curve exhibits the same leading-order scaling yielding a central charge
$c\simeq 0.5$. However, the offset with respect to no-defect for
$b_{\epsilon}=-1$ is $\Delta S(1/2)=-1/2+\ln 2$ when the total system is mixed
and 0 when pure. On the other hand, the corresponding offsets for both pure
and mixed cases for $b_{\sigma}=1$ are $\Delta S(1/2)/2=-1/4+(\ln 2)/2$.

First, we compute the symmetric EE (Fig. 2). The results for $b_{\epsilon}=1$
(no defect) are shown with blue crosses. The symmetric EE exhibits the
expected logarithmic scaling: $S_{s}(b_{\epsilon}=1)$
$=\frac{c}{3}\ln\big{[}\frac{L}{\pi}\sin\big{(}\frac{\pi r}{L}\big{)}\big{]}$
$+S_{0}$ for all values of $r/L$. We set the lattice spacing to 1 throughout
this work. Fitting this expression yields the expected central charge $c\simeq
0.5$ and $S_{0}\simeq 0.478$. The maroon (green) hexagons corresponds to the
symmetric EE for $b_{\epsilon}=-1$ ($b_{\sigma}=1$) when the total system is
in an incoherent superposition of the zero-energy states being filled and
empty. Compared to $S_{s}(b_{\epsilon}=1)$, the symmetric EEs ($S_{s}^{m},m$
denoting the total system being mixed) for the two cases get an additional
contribution of $\Delta S(r/L)$ [26] and $\Delta S(r/L)/2+(\ln 2)/2$. Thus,
when $r/L\ll 1$, for both $b_{\epsilon}=-1$ and $b_{\sigma}=1$, $S_{s}^{m}$
exhibits the expected logarithmic dependence with $c\simeq 0.5$ (left panel).
However, the offset, $S_{0}$, for $S_{s}^{m}(b_{\sigma}=1)$ is $(\ln 2)/2$
higher. This higher offset unambiguously distinguishes the topological defect
from the energy defect and is because of the localized unpaired Majorana zero
mode at the center of the subsystem (the result will be the same as long as
the defect lies within and not at the edge of the subsystem). This is
consistent with the identification of this defect problem to a boundary CFT
problem at the ‘continuous Neumann boundary fixed-point’ after folding, with
the corresponding $g$-function $=\sqrt{2}$ [16] (see also Ref. [37]). This
should be compared with the ‘continuous Dirichlet boundary fixed point’
($b_{\epsilon}=-1$ case), which has $g$-function $=1$ and thus, no additional
boundary entropy contribution. For $b_{\sigma}=1$, increasing $r/L$ leads to a
further offset of $\Delta S(r/L)/2$ due to the contribution from the second
nonlocal zero-mode. Here, the factor of 1/2 accounts for the difference in the
number of nonlocal zero-modes in the $b_{\epsilon}=-1$ and $b_{\sigma}=1$
models. The maroon squares (green diamonds) correspond to the symmetric EEs
obtained by keeping the zero-energy state empty (denoted by $S_{s}^{p},p$
denoting the total system being pure; the results are identical for the filled
case). Now, the results for $b_{\epsilon}=\pm 1$ are indistinguishable.
However, compared to the case without defects, $S_{s}^{p}(b_{\sigma}=1)$
exhibits a $\Delta S(1-r/L)/2$ offset. For $r\ll L$, this again leads to the
offset, $S_{0}$, being $(\ln 2)/2$ higher than the other cases. For $r\sim L$,
the $S_{s}^{p}$ diminishes compared to the case when the total system is mixed
due to the purity of the total system.

Next, we compute the interface EE (Fig. 3), where one end of the subsystem (of
size $r$) is located at the defect and the other end sweeps through the system
(of size $L$) 555Note that the symmetric EE computation required much larger
total system sizes in order to get a good result for offset $(\ln 2)/2$.. For
$b_{\epsilon}=-1$, the results are identical to the symmetric case. This is
expected since the resulting model is just a fermionic model with periodic bc
and there is no difference between symmetric and interface EEs. The results
for $b_{\sigma}=1$ are the same for the total system in a mixed (green
hexagons) and pure (green diamonds) state. Then, $S_{\cal
I}^{m,p}(b_{\sigma}=1)$ is given by the EE of an Ising chain of length $L-1/2$
without any defects [$S_{{\cal I},L-1/2}(b_{\epsilon}=1)$] together with an
offset. The first contribution to this offset is $\Delta S[r/(L-1/2)]$. It
arises as the subsystem size grows and becomes aware of the nonlocal zero-mode
$\Lambda(b_{\sigma})$ in the chain of length $L-1/2$. The second contribution,
$\delta_{r,L}(\ln 2)/2$, arises only when $r=L$. This is due the localized
zero-mode, $\gamma_{2L}$, which contributes only when the subsystem covers the
entire system. The resulting prediction is shown with orange pluses. For
comparison, we have also shown the curves obtained by computing the
corresponding offsets for the system-size $L$ (green pluses). For $r$ not
close to $L$, both predictions work well. However, for $L-r\sim 1$, only the
computation with system-size $L-1/2$ leads to the correct predictions (see
inset). The field-theory computations are usually done for $1\ll r,L$ with $r$
not too close to $L$. For definiteness, we consider the scaling of $S_{\cal
I}(r=L/2)$ with $\ln L$ [19, 21]. Fitting to $S_{\cal I}(r=L/2)=(c/3)\ln
L+\tilde{S}_{0}$ 666Note that for a periodic (open) system, the coefficient of
$\ln L$ is $c/3$ $(c/6)$. yields the expected central charge $c\simeq 0.5$ for
all the curves. Recall that for $0\leq b_{\sigma}<1$, the coefficient of $\ln
L$ is $c_{\rm eff}/3$, where the ‘effective central charge’ $c_{\rm
eff}\in[0,1/2)$ [22, 21]. The difference of the offsets, $\tilde{S}_{0}$,
between $b_{\epsilon}=-1$ and $b_{\epsilon}=1$ cases is $2\Delta_{I}=\Delta
S(1/2)$ $=-1/2+\ln 2$ or 0 depending on the total system being in a mixed or
pure state. For the topological case, both pure and mixed states, the
corresponding offset difference is $\Delta_{I}=\Delta S(1/2)/2=-1/4+(\ln
2)/2$. Importantly, this offset occurs entirely due to a ‘finite-size effect’
correction arising due to the existence of nonlocal zero-modes and bears no
relationship to the specific modular S-matrix elements predicted in Refs. [19,
21].

To summarize, we computed the symmetric and interface EEs for the Ising CFT
with a topological defect taking into account the subtle effects of the zero-
modes on the EEs. We showed that while both the EEs exhibit identical leading-
order logarithmic scaling, the subleading $O(1)$ corrections are of completely
different origin. The subleading term [$=(\ln 2)/2$] for the symmetric EE is
related to the $g$-function of the corresponding defect at the boundary in the
folded picture. However, the corresponding term [$\Delta S(r/L)/2$] in the
interface EE arises only when the subsystem occupies a finite fraction of the
total system and is entirely due to the local and nonlocal zero-modes of the
topological defect Hamiltonian. In the limit $r\ll L$, there is no additional
offset compared with the case without defect. The interface EE result is in
sharp contrast with the existing predictions in terms of the modular S-matrix
of the CFT [21, 19], which predict an offset equal to $-\ln 2$ instead. We
also computed the EEs in the case of open/free boundary conditions with a
defect at the center of the chain 777See Supplementary Material for details.
The results for symmetric EE were identical to that obtained for the periodic
chain. The offset interface EE is the same as for the periodic chain when the
total system was pure. The total system being mixed contributes another $(\ln
2)/2$ to the offset.

Several features of the topological defect persist away from the topological
point, i.e., $b_{\sigma}\neq 1$, and even away from the conformal critical
point. We plan to address some of these questions elsewhere. Defect
Hamiltonians for other rational CFTs contain more complicated set of zero
modes [41]. The question of subleading corrections in EEs for these models
remains open. Finally, advancements in measurement of Renyi entropies in
engineered quantum systems [42, 43] and quantum simulation [44] can lead to
potential verification of our analytical predictions.

We thank Natan Andrei, Pasquale Calabrese, Christopher Herzog and Ingo Peschel
for discussions. We also thank David Rogerson and Frank Pollmann for
discussions and collaborations on a related project. AR acknowledges support
from a grant from the Simons Foundation (825876, TDN). HS was supported by ERC
Advanced Grant NuQFT.

## References

  * Holzhey  _et al._ [1994] C. Holzhey, F. Larsen, and F. Wilczek, Geometric and renormalized entropy in conformal field theory, Nuclear Physics B 424, 443 (1994).
  * Calabrese and Cardy [2004] P. Calabrese and J. Cardy, Entanglement entropy and quantum field theory, Journal of Statistical Mechanics: Theory and Experiment 2004, P06002 (2004).
  * Affleck and Ludwig [1991] I. Affleck and A. W. W. Ludwig, Universal noninteger “ground-state degeneracy” in critical quantum systems, Phys. Rev. Lett. 67, 161 (1991).
  * Calabrese and Cardy [2009] P. Calabrese and J. Cardy, Entanglement entropy and conformal field theory, J. Phys. A42, 504005 (2009), arXiv:0905.4013 [cond-mat.stat-mech] .
  * Affleck  _et al._ [2009] I. Affleck, N. Laflorencie, and E. S. Sørensen, Entanglement entropy in quantum impurity systems and systems with boundaries, Journal of Physics A: Mathematical and Theoretical 42, 504009 (2009).
  * Roy  _et al._ [2020] A. Roy, F. Pollmann, and H. Saleur, Entanglement Hamiltonian of the 1+1-dimensional free, compactified boson conformal field theory, J. Stat. Mech. 2008, 083104 (2020), arXiv:2004.14370 [cond-mat.stat-mech] .
  * Roy  _et al._ [2021] A. Roy, D. Schuricht, J. Hauschild, F. Pollmann, and H. Saleur, The quantum sine-Gordon model with quantum circuits, Nucl. Phys. B 968, 115445 (2021), arXiv:2007.06874 [quant-ph] .
  * Petkova and Zuber [2001] V. B. Petkova and J. B. Zuber, Generalized twisted partition functions, Phys. Lett. B 504, 157 (2001), arXiv:hep-th/0011021 .
  * Bachas  _et al._ [2002] C. Bachas, J. de Boer, R. Dijkgraaf, and H. Ooguri, Permeable conformal walls and holography, JHEP 06, 027, arXiv:hep-th/0111210 .
  * Fröhlich  _et al._ [2004] J. Fröhlich, J. Fuchs, I. Runkel, and C. Schweigert, Kramers-wannier duality from conformal defects, Phys. Rev. Lett. 93, 070601 (2004).
  * Frohlich  _et al._ [2007] J. Frohlich, J. Fuchs, I. Runkel, and C. Schweigert, Duality and defects in rational conformal field theory, Nucl. Phys. B 763, 354 (2007), arXiv:hep-th/0607247 .
  * Aasen  _et al._ [2016] D. Aasen, R. S. K. Mong, and P. Fendley, Topological Defects on the Lattice I: The Ising model, J. Phys. A 49, 354001 (2016), arXiv:1601.07185 [cond-mat.stat-mech] .
  * Kramers and Wannier [1941] H. A. Kramers and G. H. Wannier, Statistics of the two-dimensional ferromagnet. part i, Phys. Rev. 60, 252 (1941).
  * Savit [1980] R. Savit, Duality in field theory and statistical systems, Rev. Mod. Phys. 52, 453 (1980).
  * Buican and Gromov [2017] M. Buican and A. Gromov, Anyonic Chains, Topological Defects, and Conformal Field Theory, Commun. Math. Phys. 356, 1017 (2017), arXiv:1701.02800 [hep-th] .
  * Oshikawa and Affleck [1997] M. Oshikawa and I. Affleck, Boundary conformal field theory approach to the critical two-dimensional ising model with a defect line, Nuclear Physics B 495, 533 (1997).
  * Saleur [1998] H. Saleur, Lectures on nonperturbative field theory and quantum impurity problems, (1998), arXiv:cond-mat/9812110 .
  * Saleur [2000] H. Saleur, Lectures on nonperturbative field theory and quantum impurity problems: Part 2, (2000), arXiv:cond-mat/0007309 .
  * Gutperle and Miller [2016] M. Gutperle and J. D. Miller, A note on entanglement entropy for topological interfaces in RCFTs, JHEP 04, 176, arXiv:1512.07241 [hep-th] .
  * Sakai and Satoh [2008] K. Sakai and Y. Satoh, Entanglement through conformal interfaces, JHEP 12, 001, arXiv:0809.4548 [hep-th] .
  * Brehm and Brunner [2015] E. M. Brehm and I. Brunner, Entanglement entropy through conformal interfaces in the 2D Ising model, JHEP 09, 080, arXiv:1505.02647 [hep-th] .
  * Eisler and Peschel [2010] V. Eisler and I. Peschel, Entanglement in fermionic chains with interface defects, Annalen der Physik 522, 679 (2010).
  * Peschel and Eisler [2012] I. Peschel and V. Eisler, Exact results for the entanglement across defects in critical chains, Journal of Physics A: Mathematical and Theoretical 45, 155301 (2012).
  * Calabrese  _et al._ [2012] P. Calabrese, M. Mintchev, and E. Vicari, Entanglement Entropy of Quantum Wire Junctions, J. Phys. A 45, 105206 (2012), arXiv:1110.5713 [cond-mat.stat-mech] .
  * Herzog and Nishioka [2013] C. P. Herzog and T. Nishioka, Entanglement Entropy of a Massive Fermion on a Torus, JHEP 03, 077, arXiv:1301.0336 [hep-th] .
  * Klich  _et al._ [2017] I. Klich, D. Vaman, and G. Wong, Entanglement hamiltonians for chiral fermions with zero modes, Phys. Rev. Lett. 119, 120401 (2017).
  * Grimm [2002] U. Grimm, Spectrum of a duality twisted Ising quantum chain, J. Phys. A 35, L25 (2002), arXiv:hep-th/0111157 .
  * Vidal  _et al._ [2003] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Entanglement in quantum critical phenomena, Phys. Rev. Lett. 90, 227902 (2003), arXiv:quant-ph/0211074 .
  * Peschel [2003] I. Peschel, Calculation of reduced density matrices from correlation functions, Journal of Physics A: Mathematical and General 36, L205 (2003).
  * Latorre  _et al._ [2004] J. I. Latorre, E. Rico, and G. Vidal, Ground state entanglement in quantum spin chains, Quantum Inf. Comput. 4, 48 (2004).
  * Henkel  _et al._ [1989] M. Henkel, A. Patkós, and M. Schlottmann, The ising quantum chain with defects (i). the exact solution, Nuclear Physics 314, 609 (1989).
  * Baake  _et al._ [1989] M. Baake, P. Chaselon, and M. Schlottmann, The ising quantum chain with defects (ii). the so(2n)kac-moody spectra, Nuclear Physics B 314, 625 (1989).
  * Note [1] Note that $h_{\epsilon,\sigma}$ are marginal perturbations which, in general, lead to continuously varying critical exponents [45].
  * Note [2] This duality defect Hamiltonian is related by a local unitary rotation on the $L^{\rm th}$ spin to the one considered in Refs. [16], which has $\sigma_{L}^{z}\sigma_{1}^{x}$ interaction. We do not use this alternate form since it no longer leads to a bilinear Hamiltonian under JW transformation and cannot be solved by free-fermion techniques.
  * Note [3] Equivalently, we analyze the mixed sector Hamiltonian [27]: $H_{\epsilon,\sigma}^{m}=P_{+}H_{\epsilon,\sigma}^{f}(Q=+1)+P_{-}H_{\epsilon,\sigma}^{f}(Q=-1)$.
  * Note [4] See Supplementary Material for details.
  * Alba  _et al._ [2017] V. Alba, P. Calabrese, and E. Tonni, Entanglement spectrum degeneracy and the cardy formula in 1+1 dimensional conformal field theories, Journal of Physics A: Mathematical and Theoretical 51, 024001 (2017).
  * Note [5] Note that the symmetric EE computation required much larger total system sizes in order to get a good result for offset $(\mathop{ln}\nolimits 2)/2$.
  * Note [6] Note that for a periodic (open) system, the coefficient of $\mathop{ln}\nolimits L$ is $c/3$ $(c/6)$.
  * Note [7] See Supplementary Material for details.
  * Belletête  _et al._ [2020] J. Belletête, A. M. Gainutdinov, J. L. Jacobsen, H. Saleur, and T. S. Tavares, Topological defects in periodic rsos models and anyonic chains (2020), arXiv:2003.11293 [math-ph] .
  * Islam  _et al._ [2015] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and M. Greiner, Measuring entanglement entropy in a quantum many-body system, Nature 528, 77 (2015).
  * Brydges  _et al._ [2019] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Probing rényi entanglement entropy via randomized measurements, Science 364, 260–263 (2019).
  * Monroe  _et al._ [2021] C. Monroe, W. Campbell, L.-M. Duan, Z.-X. Gong, A. Gorshkov, P. Hess, R. Islam, K. Kim, N. Linke, G. Pagano, and et al., Programmable quantum simulations of spin systems with trapped ions, Reviews of Modern Physics 93, 10.1103/revmodphys.93.025001 (2021).
  * Cardy [1987] J. L. Cardy, Continuously varying exponents and the value of the central charge, Journal of Physics A: Mathematical and General 20, L891 (1987).