# Evaluation of Local Model-Agnostic Explanations Using Ground Truth

Amir Hossein Akhavan Rahnama
KTH Royal Institute of Technology
<EMAIL_ADDRESS>
&Jüdith Bütepage
Spotify
<EMAIL_ADDRESS>
&Pierre Geurts
Universite de Liège, Belgium
<EMAIL_ADDRESS>
&Henrik Böstrom
KTH Royal Institute of Technology
<EMAIL_ADDRESS>

###### Abstract

Explanation techniques are commonly evaluated using human-grounded methods,
limiting the possibilities for large-scale evaluations and rapid progress in
the development of new techniques. We propose a functionally-grounded
evaluation procedure for local model-agnostic explanation techniques. In our
approach, we generate ground truth for explanations when the black-box model
is Logistic Regression and Gaussian Naive Bayes and compare how similar each
explanation is to the extracted ground truth. In our empirical study,
explanations of Local Interpretable Model-agnostic Explanations (LIME),
SHapley Additive exPlanations (SHAP), and Local Permutation Importance (LPI)
are compared in terms of how similar they are to the extracted ground truth.
In the case of Logistic Regression, we find that the performance of the
explanation techniques is highly dependent on the normalization of the data.
In contrast, Local Permutation Importance outperforms the other techniques on
Naive Bayes, irrespective of normalization. We hope that this work lays the
foundation for further research into functionally-grounded evaluation methods
for explanation techniques.

##  1 Introduction

Local model-agnostic explanation techniques estimate the importance of
features for the prediction made by a black-box model given a single input
instance. Explanation techniques are usually evaluated using two different
types of methods. Human-grounded evaluation methods use the judgment of human
subjects to evaluate the accuracy of explanations. The subjectivity of human
judgment, time, and costs make human-ground evaluation often infeasible. These
challenges can be mitigated by functionally-grounded evaluation methods in
which different proxy measures are used to evaluate explanation techniques
systematically (see Doshi-Velez and Kim (2017) for more details). The major
challenge of functionally-grounded evaluations of explanation techniques is
the lack of the ground truth feature importance scores (see Hooker et al.
(2018) for a detailed discussion). In this work, we propose a procedure that
extracts the ground truth of feature importance scores for the Logistic
Regression and the Naive Bayes models. The extracted ground truth is used to
compare different explanation techniques in terms of how similar the
explanations are to the ground truth.

Local additive model-agnostic explanation techniques are a category of
explanations in which the black-box model’s predicted probability of a single
instance for a given class is decomposed into an additive sum of importance
scores for all features. Local Interpretable Model-agnostic Explanations
(LIME) (Ribeiro et al. (2016)) and SHapley Additive exPlanations (SHAP)
(Lundberg and Lee (2017)) are examples of widely used techniques to obtain
additive explanations. Local Permutation Importance (LPI) defines the
importance score to be the decrease in black-box score when a single feature
value is randomly shuffled (Casalicchio et al. (2018)).

The need to find ground truth for explanations is an important step toward
evaluating different local additive explanation techniques. One major
challenge in finding the ground truth for local additive explanations is that
the additive decomposition of a predicted value into the sum of feature
importance scores is imposed by the explanation techniques while the black-box
model is usually more complex. There are two categories of functionally-
grounded methods that evaluate the explanations without the need for ground
truth: ablation and robustness evaluation techniques. Ablation techniques
(Hooker et al. (2018)) measure the sensitivity of the model’s predicted value
after the removal of features with the largest importance scores. The
rationale behind this evaluation method is that the more substantial the
change in the model’s predicted value, the more accurate are the explanations.
Another set of evaluation methods analyze the robustness of explanations
following the introduction of noise or redundant features into the explained
instance (Alvarez-Melis and Jaakkola (2018); Camburu et al. (2019)). The main
rationale behind this category of evaluation techniques is that a robust
explanation technique should be resilient with respect to noise and assign
relatively low importance scores to redundant features. The main limitation of
ablation and robustness evaluation techniques is that changing feature values
can shift the distribution of data which can lead to highly uncertain
predictions by the black-box model. In addition, there is no objective way to
conclude that a robust explanation is accurate or to measure the influence of
correlated features in ablation techniques.

In our study, we propose an alternative approach to evaluate explanation
techniques that directly involve the extraction of ground truth for
explanations. Our evaluation procedure extracts the ground truth for a single
instance from the log odds ratio of predictions. To this end, we focus on two
models, namely Logistic Regression and Naive Bayes. The extracted ground truth
is then used to compare different explanation techniques using Spearman’s rank
correlation. We study two local additive explanation techniques, i.e. Local
Interpretable Model-agnostic Explanations (LIME) and SHapley Additive
exPlanations (SHAP), along with a widely used explanation technique, namely
Local Permutation Importance (LPI) (Casalicchio et al. (2018)).

In the next section, a summary of related work with a focus on the evaluation
of different explanation techniques is presented. In Section 3, our evaluation
procedure is presented. The result from the empirical study including the
comparison of different techniques is discussed in Section 4. In Section 5,
discussions and concluding remarks including future directions for possible
extensions of this study are presented.

##  2 Related Work

Functionally-grounded evaluation methods of explanation techniques fall into
two main categories: ablation methods and analysis of robustness. In Hooker et
al. (2018), RemOve And Retrain (ROAR) is proposed to evaluate different
explanations techniques that output saliency maps for neural network models.
In this framework, the black-box model is retrained on an the auxiliary
dataset that stripes pixels deemed important by explanations from the input.
The accuracy of the explanations are calculated as the difference between the
predicted class probabilities of the original and the auxiliary instance. The
results show that explanations of VarGrad (Adebayo et al. (2018)) and
SmoothGrad-Squared (Smilkov et al. (2017)) outperform other explanation
techniques. ROAR is computationally exhaustive as the underlying model needs
to be retrained on a large auxiliary dataset. Another limitation of this
evaluation methods is that the new auxiliary data and retrained black-box
model deviate from the original data and model and can inhibit new structures
and properties.

Explanations are evaluated with regards to their robustness following the
introduction of Lipschitz noise to the explained instances in Alvarez-Melis
and Jaakkola (2018). The authors empirically show that some explanation
techniques including Local Interpretable Model-agnostic Explanations (LIME)
are not robust to the Lipschitz noise and therefore have low robustness. In
Camburu et al. (2019), the robustness of Local Interpretable Model-agnostic
Explanations (LIME), SHapley Additive exPlanations (SHAP) and Learn To eXplain
(L2X) (Chen et al. (2018)) are studied following the introduction of
irrelevant features to the explained instance in a sentiment analysis task.
According to the results, L2X explanations outperform other explanation
techniques across different measures of robustness that are proposed in the
study. One limitation of this category of evaluation techniques is that
changes in the feature values of instances can introduce out-of-distribution
data. As a result, predictions by the black-box model can have high
uncertainty as discussed in Lakkaraju et al. (2020); Rahnama and Boström
(2019).

In contrast to the studies mentioned in this section, we propose a
functionally-grounded procedure that directly measures the accuracy of
explanations. Our proposed procedure extracts the ground truth systematically
and measures the similarity between an explanation and the ground truth.

##  3 Methodology

We begin by introducing the notation in Section 3.1. In Section 3.2, we
describe our proposed approach to extract the ground truth for Logistic
Regression and Naive Bayes from the intrinsic additive decomposition in their
log odds ratios. The proposed evaluation procedure of explanation techniques
is discussed in Section 3.3.

###  3.1 Notation

In this section, we formally define the notation used throughout this paper.
Let $X\in\mathbb{R}^{m\times n}$ be a data matrix of $m$ instances with $n$
corresponding features. Suppose $x\in X$ is an instance where
$x\in\mathbb{R}^{n}$. In local additive explanations111Throughout this paper
and for brevity, the notion of local additive explanation and local additive
model-agnostic explanation are used interchangeably., the conditional
probability $P$ of model $f$ decomposes the predicted value for a single
instance given a specific class into a linear additive combination. Let
$x\in\mathbb{R}^{n}$ and $p_{x}^{C}$ be the predicted probability score of the
instance $x$ by model $f$ with respect to class $C$, hence:
$p_{x}^{C}=\phi_{1}+...+\phi_{n}$ where explanations can be represented as
$\Phi=\\{(\textrm{feature}_{j},\phi_{i})|j=1,...,N\\}$ where $\phi_{j}$ is an
importance score for feature $j$ with respect to $p_{x}^{C}$. In practice and
due to cognitive limitations of humans, a subset of $\Phi$ is often selected
including the top-$K$ ranked elements by absolute value of importance scores
where $K\ll N$. In addition, the extracted ground truth is denoted a
$\Lambda=(\lambda_{1},...,\lambda_{N})\in\mathbb{R}^{n}$ where $\lambda_{j}$
represents the ground truth for feature $j$.

###  3.2 Ground Truth

As stated earlier, additive explanations decompose class prediction
probability scores into an additive sum of feature importance scores. However,
the black-box model usually has no intrinsic additive structure in their
prediction functions. In this section, we present an alternative solution to
this problem by using the log odds ratios of an instance with respect to a
specific class instead of class prediction probability. Since the additive
structures are present in the log odds predictions of the Logistic Regression
and the Naive Bayes model, contributions of each features to the predicted log
odds can be extracted and used as the ground truth.

####  3.2.1 Logistic Regression

Given a Logistic Regression model in form of
$y^{(i)}=w_{0}+w_{1}x_{i,1}+...+w_{n}x_{i,N}$ and instance
$x_{i}=(x_{i,1},...,x_{i,n})$,

$P(y^{(i)}=1|x_{i},W)=\frac{1}{1+\textrm{e}^{{-\sum_{j=0}^{i=N}w_{j}x_{i,j}}}}$
(1)

where $W=(w_{1},...,w_{n})$. Equation 2 show that the prediction probability
function of Logistic Regression has no additive structure. However, we can
derive a inherently additive decomposition when deriving the log odds ratio
for $X$ with respect to class 1:

$\textrm{log}\frac{P(y^{(i)}=1|x_{i},W)}{P(y^{(i)}=0)|x_{i},W)}=\sum_{j=0}^{i=N}w_{j}x_{i,j}$
(2)

In Equation 2, it is clear that the contribution of feature $j$ to the log
odds of an instance with respect to class 1 is $w_{j}x_{i,j}$ and we consider
this to be the ground truth importance score $\lambda_{j}$ of feature $j$.

####  3.2.2 Gaussian Naive Bayes

Given input $X$ and class $C$, the prediction function of Gaussian Naive Bayes
model222The notion of Gaussian Naive Bayes and Naive Bayes are used
interchangeably throughout the paper for brevity., for class $C$ and feature
$j$ where $j=1,...,N$ and $i=1,...,M$,

$P(x_{i,j}|Y=C)=\frac{1}{\sqrt{2\pi\sigma_{j,C}}}e^{-\frac{(x_{i,j}-\mu_{j}^{C})^{2}}{2{\sigma_{j}^{C}}^{2}}}\sim\mathcal{N}(x_{i,j}|\mu_{j}^{C},\sigma_{j}^{C})$
(3)

Similar to the case of Logistic Regression, no additive decomposition can be
seen in the probability function. The log odds ratios of an instance $x_{i}$
with respect to class 1 has an intrinsic natural additive decomposition,

$\textrm{log}\frac{P(y^{(i)}=1|x_{i})}{P(y^{(i)}=0|x_{i})}=\textrm{log}\frac{P(y^{(i)}=1)}{P(y^{(i)}=0)}+\sum_{j=1}^{N}\textrm{log}\frac{\mathcal{N}(x_{i,j}|\mu_{j}^{1},\sigma_{j}^{1})}{\mathcal{N}(x_{i,j}|\mu_{j}^{0},\sigma_{j}^{0})}$
(4)

In Equation 4, the contribution of feature $j$ to the log odds ratio of
instance $x_{i}$ with respect to class 1,
$\textrm{log}\frac{\mathcal{N}(x_{i,j}|\mu_{j,1},\sigma_{j,1})}{\mathcal{N}(x_{i,j}|\mu_{j,0},\sigma_{j,0})}$
is used as the ground truth $\lambda_{j}$ for that feature. Figure 1 shows an
example of the ground truth explanation obtained for a single instance from
Pima Indians dataset for Logistic Regression and Naive Bayes respectively.

Figure 1: Ground truth obtained for Logistic Regression model and Naive Bayes
on a single instance from Pimal Indians dataset.

###  3.3 Evaluation procedure

The outline of our proposed evaluation procedure that measures the similarity
between ground truth and explanations is as follows (see Algorithm 1). Given
an input $x$, and a black-box model $f$ (Logistic Regression and Naive
Bayes)333We use the term black-box for Logistic Regression and Naive Bayes to
be consistent with the terminology used in the explainability literature., we
extract the ground truth $\Lambda$ from the log odds ratio of $x$ with respect
to class $C$. In the next step, explanation technique $g$ is selected and the
feature importance $\Phi$ of the log odds ratios of $x$ with respect to class
$C$ is obtained444In Section 4.4, we have performed experiments in which we
show the effect of obtaining feature importance $\Phi$ from the prediction
probability scores of the black-box model on our evaluation procedure.. Given
a similarity function $k$, the algorithm calculates the similarity of
explanations and ground truth and returns the score $r_{x,f,g}$.

Algorithm 1 Evaluating Explanations

Input $x$: instance, $f$: black-box, $g$: explanation method
$k$: Spearman’s rank correlation, $\Lambda$: ground truth
Output $r_{x,f,g}$: correlation value

1:$\Phi\leftarrow g(f,x)$

2:$\Lambda\leftarrow t(f,x)$

3:$r_{x,f,g}\leftarrow k(\Phi,\Lambda)$

The choice of similarity function is essential in our evaluation procedure.
Several studies have used measures such as cosine similarity or Euclidean
distance (see Montavon et al. (2018), Yang and Kim (2019)) for measuring
similarity between explanations. In our approach, we use Spearman’s Rank
correlation (Zar (2005)) as the selected measure of similarity. The rationale
behind our choice is that we are interested in measuring the monotonic
relationship between the ground truth and explanations. In simpler terms,
small changes in the importance score of a single feature which do not result
in changes of the relative rank of that feature with respect to all other
features should not be taken into account by our choice of similarity
function. In addition, Spearman’s rank correlation is interpretable with
respect to the existence of a significance threshold for correlation values
and the direction of the correlation. 555Negative Spearman’s rank correlation
indicate the decreasing monotonic relationship of the two vectors.

We are interested in comparing different techniques based on their similarity
to the ground truth across all datasets. Let $d_{i}$ be the $i$-th dataset
($i=1,...,T$) where $T$ is the total number of datasets. Suppose
$R^{i}_{f,g}=\\{r_{x_{k},f,g}|k=1,...,m\\}$ is obtained by applying Algorithm
1 for $m$ test instances in the $i$-th dataset. For dataset $i$, the median
value of $R^{i}_{f,g}$ scores from each explanation technique is calculated,
and techniques are then ranked in a descending using this value for model $f$.
In simpler terms, the techniques that generate higher median rank correlation
values will be ranked lower in the $i$-th dataset when explaining model $f$.
The overall ranks for each technique are the average ranks over $T$ datasets
and therefore explanation techniques with a lower average ranks outperform the
techniques with a higher average ranks for model $f$.

##  4 Empirical Investigation

###  4.1 Experimental setup

The proposed evaluated procedure is applied to explanations of Logistic
Regression and Naive Bayes models. We compare the explanations extracted with
LIME, SHAP, and LPI to the ground truth explanations. More specifically, we
are using TabularLIME when obtaining LIME explanations for LIME and KernelShap
for SHAP explanations. We evaluate our procedure on 20 different tabular
datasets, all concerning binary classification tasks666All datasets are
publicly available and follow the ethics and data protection guidelines.
Information about the datasets is provided in the Appendix. Code to replicate
the experiments is attached as supplementary material.. The experiments were
run on a Ubuntu 18.10 Intel(R) Core(TM) i9-7980XE 2.60GHz CPU with 18 cores
and 20GB memory. In the experiments, numerical features are standardized and
categorical features are encoded as one-hot vectors. For datasets for which no
separate test set has been provided at the source, a hold-out set of 25% was
used. The random search was used for tuning the hyper-parameters of the
Logistic Regression model 777Hyper parameters were chosen after 100 trials
with the hyper-parameter space consisting of L1 and L2 regularization with the
regularization parameter selected from a uniform distribution of values
between 0 to 4.. The number of samples generated for LIME and SHAP is 5000. In
the case of LPI, it is the size of the training set as suggested in
Casalicchio et al. (2018). This means the sample size of LPI is half the size
of LIME and SHAP samples on average.

Table 5888The table is in the Appendix. includes the test accuracy of Logistic
Regression and Naive Bayes models along with the abbreviated names of each
dataset used in the empirical experiments of this study. As stated earlier,
since our ground truth are extracted from log odds ratios of instances with
respect to a specific class, all examined explanation techniques are adjusted
to decompose the log odds ratios of a given instance. LPI explanations are
adjusted to rank importance scores for the difference in the log odds
prediction following the random permutations of features.

In Section 4.2, we compare the selected explanations techniques against our
ground truth. As stated earlier, we compare the log odds ratio of a single
instance with respect to a specific class instead of class probability scores.
This is necessary as the extracted ground truth is based the log odds ratio
(see Equation 2 and 4). In 4.3, we show that the comparison is substantially
affected by the choice of pre-processing technique that is applied on the
numerical features in each dataset. In section 4.4, we compare all techniques
with the ground truth when explanations are obtained from predicted
probability scores instead of log odds ratio.

Table 1: Average rank of explanation techniques based on the median rank correlation across all datasets  | Logistic Regression | Naive Bayes
---|---|---
| LIME | SHAP | LPI | LIME | SHAP | LPI
Average Rank (All datasets) | 2.1 | 2.08 | 1.82 | 2 | 2.7 | 1.3
Figure 2: Box-plots of Spearman’s rank correlation of explanation techniques
when the underlying model is Logistic Regression (top) and Naive Bayes
(bottom). The shaded area visualizes the threshold where $\textrm{corr}>0.7$.

###  4.2 General performance

In this section, we perform experiments that evaluate explanation techniques
against the extracted ground truth using rank correlation. To the best of our
knowledge, this is the first study that examines the explanations of log odds
ratios of explanation techniques.

The box-plot including quartiles of the Spearman rank correlation values for
each technique is shown for using logistic regression in Figure 2. The shaded
area in the figure represent the significant correlation threshold where
$\textrm{corr}>0.7$. In some cases, the explanation techniques provide
explanations that correlate significantly with the ground truth, e.g. SHAP
explanations on the Banknote dataset for the Naive Bayes model and LPI
explanations on the Haberman dataset for the case of the Logistic Regression
model. To our surprise, LIME explanations show median values below or close to
zero across many datasets with large standard deviations, e.g. breast cancer,
Pima Indians, Banknote, Iris, and Haberman. One important pattern that can be
seen in Figure 2 is that the Spearman’s rank correlation values from
explanation techniques are more likely to be significant where the black-box
model is Logistic Regression compared to Naive Bayes across numerous datasets,
e.g. Adult, Heart Disease, Churn and HR.

We present the average performance of the different explanation techniques in
Table 1.999The rank of explanation techniques based on median Sparman’s rank
correlation for each dataset can be found in Appendix.. As stated in Section
3.3, the explanation techniques with lower average ranks over all datasets
outperform techniques with higher average ranks. It can be seen that LPI
explanations outperform other techniques for both black-box models. In
addition, the minimum difference between average ranks of LPI and other
techniques is 3.5 times larger in the case of the Naive Bayes model compared
to the Logistic Regression model.

Figure 3: Box-plots of Spearman’s rank correlation of explanation techniques
when the black-box model is Logistic Regression (top) and Naive Bayes (bottom)
under min-max transformation. The shaded area visualizes the threshold where
$\textrm{corr}>0.7$. Figure 4: Box-plots of Spearman’s rank correlation of
explanation techniques when the black-box model is Logistic Regression (top)
and Naive Bayes (bottom) under inter-quartile transformation. The shaded area
visualizes the threshold where $\textrm{corr}>0.7$.

###  4.3 Pre-processing Effect

As stated earlier, numerical features were normalized and the encoded
categorical values are processed as one-hot vectors in the experiments of
Section 4.2. In this section, we show that different pre-processing
transformations on numerical features have a direct effect on the ranks of
explanation techniques over multiple datasets. We have included the results of
applying two other pre-processing techniques including min-max scaling and
inter-quartile pre-processing in this section.

Table 2: The average rank of explanation techniques across all datasets for each pre-processing transformation | Logistic Regression | Naive Bayes
---|---|---
| LIME | SHAP | LPI | LIME | SHAP | LPI
Normalization | 2.1 | 2.075 | 1.825 | 2 | 2.7 | 1.3
Min-max Scaling | 1.525 | 1.925 | 2.55 | 1.925 | 2.65 | 1.425
Inter-quartile Normalization | 2.275 | 1.375 | 2.35 | 1.975 | 2.65 | 1.375

The box-plots of rank correlations for each explanation techniques are
visualized after min-max (Figure 3) and inter-quartile scaling (Figure 4).
Table 2 contains the summary of the rank of all explanation techniques for
each pre-processing transformation averaged over all datasets. SHAP and LIME
explanations marginally outperform the other techniques in the case of min-max
and inter-quartile scaling respectively. Examples of this trend for LIME
explanations can be seen in datasets such Attrition, Breast Cancer, Pima
Indians and Banknote data sets when explaining Logistic Regression as can be
seen in Figure 3. Similarly, SHAP explanations show substantial improvements
in explaining Logistic Regression for data sets such as Breast Cancer,
Spambase, Seismic and Thera as shown in Figure 4. LPI explanations outperform
all explanation techniques across all pre-processing techniques applied when
examining the Naive Bayes model.

###  4.4 Comparison based on prediction probability

In Section 3.3, we derived the feature importance values $\Phi$ from
explanation technique $g$ based on the log odds ratio of instance $x$ with
respect to class $C$. In this section, we investigate whether obtaining
feature importance $\Phi$ from explanation techniques based on the prediction
probability score instead of the log odds ratio value affects our evaluation
of explanation techniques. Table 3 shows the results for average rank of
explanation techniques over all datasets. As it can be seen from the table,
using prediction probability has insignificant effects on the average ranks of
explanations across all datasets and therefore similar conclusions can be made
for the most outperforming technique when explaining Naive Bayes and Logistic
Regression.

Table 3: The average rank of explanation techniques with normalization where explanation techniques decompose probability scores | Logistic Regression | Naive Bayes
---|---|---
| LIME | SHAP | LPI | LIME | SHAP | LPI
Standardization | 2.1 | 2.12 | 1.78 | 1.98 | 2.75 | 1.27
Min-max Scaling | 1.48 | 1.92 | 2.6 | 2 | 2.7 | 1.3
Inter-quartile Normalization | 2.22 | 1.48 | 2.3 | 2.02 | 2.7 | 1.27

##  5 Concluding remarks

In this study, we propose an evaluation procedure for comparing local additive
explanations to ground truth local explanations. We showed two examples of
extracting ground truth from white-box models that have intrinsic additive
representations, namely Logistic Regression and Naive Bayes. Our comparison
procedure is based on Spearman’s Rank Correlation between the explanation
vectors obtained from an explanation technique against our ground truth. Using
our method, we were able to conclude what relative performance can be expected
of the explanation techniques for the considered model types. The empirical
experiments compared explanations obtained from different techniques such as
LIME, SHAP and LPI on 20 tabular datasets. Based on our comparison, each
explanation technique can outperform other techniques depending on the choice
of pre-processing transformation applied to the datasets. In the case of Naive
Bayes, the explanations of LPI outperform the other techniques across all pre-
processing techniques.

The limitations of our proposed approach are two-fold. Firstly, we have
focused on comparing the accuracy of explanation of Logistic Regression and
Naive Bayes and the comparison of explanations of other class of black-box
models are omitted. In addition, we have focused on the binary classification
tasks in our empirical study.

One obvious future direction is identifying procedures in which ground truth
can be extracted from other classes of machine learning models. The challenge
is to derive the ground truth in the models where intrinsic additive
structures are not present such as in decision trees.

On a final note, it should be emphasized that the proposed framework should be
seen as complementary to human-grounded evaluations. The functionally-grounded
evaluation methods could be seen as means to make the development of new
candidate techniques more efficient; they allow for rejecting some candidates
early on, but the finally produced candidates will most likely still need to
be qualitatively evaluated in a user-centered context in the end.

## References

  * Doshi-Velez and Kim (2017) Finale Doshi-Velez and Been Kim.  Towards a rigorous science of interpretable machine learning.  _arXiv preprint arXiv:1702.08608_ , 2017.
  * Hooker et al. (2018) Sara Hooker, Dumitru Erhan, Pieter-Jan Kindermans, and Been Kim.  A benchmark for interpretability methods in deep neural networks.  _arXiv preprint arXiv:1806.10758_ , 2018.
  * Ribeiro et al. (2016) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin.  " why should i trust you?" explaining the predictions of any classifier.  In _Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining_ , pages 1135–1144, 2016.
  * Lundberg and Lee (2017) Scott Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  _arXiv preprint arXiv:1705.07874_ , 2017.
  * Casalicchio et al. (2018) Giuseppe Casalicchio, Christoph Molnar, and Bernd Bischl.  Visualizing the feature importance for black box models.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , pages 655–670. Springer, 2018.
  * Alvarez-Melis and Jaakkola (2018) David Alvarez-Melis and Tommi S Jaakkola.  On the robustness of interpretability methods.  _arXiv preprint arXiv:1806.08049_ , 2018.
  * Camburu et al. (2019) Oana-Maria Camburu, Eleonora Giunchiglia, Jakob Foerster, Thomas Lukasiewicz, and Phil Blunsom.  Can i trust the explainer? verifying post-hoc explanatory methods.  _arXiv preprint arXiv:1910.02065_ , 2019.
  * Adebayo et al. (2018) Julius Adebayo, Justin Gilmer, Michael Muelly, Ian Goodfellow, Moritz Hardt, and Been Kim.  Sanity checks for saliency maps.  _arXiv preprint arXiv:1810.03292_ , 2018.
  * Smilkov et al. (2017) Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda Viégas, and Martin Wattenberg.  Smoothgrad: removing noise by adding noise.  _arXiv preprint arXiv:1706.03825_ , 2017.
  * Chen et al. (2018) Jianbo Chen, Le Song, Martin Wainwright, and Michael Jordan.  Learning to explain: An information-theoretic perspective on model interpretation.  In _International Conference on Machine Learning_ , pages 883–892. PMLR, 2018.
  * Lakkaraju et al. (2020) Himabindu Lakkaraju, Nino Arsov, and Osbert Bastani.  Robust and stable black box explanations.  In _International Conference on Machine Learning_ , pages 5628–5638. PMLR, 2020.
  * Rahnama and Boström (2019) Amir Hossein Akhavan Rahnama and Henrik Boström.  A study of data and label shift in the lime framework.  _arXiv preprint arXiv:1910.14421_ , 2019.
  * Montavon et al. (2018) Grégoire Montavon, Wojciech Samek, and Klaus-Robert Müller.  Methods for interpreting and understanding deep neural networks.  _Digital Signal Processing_ , 73:1–15, 2018.
  * Yang and Kim (2019) Mengjiao Yang and Been Kim.  Benchmarking attribution methods with relative feature importance.  _arXiv preprint arXiv:1907.09701_ , 2019.
  * Zar (2005) Jerrold H Zar.  Spearman rank correlation.  _Encyclopedia of biostatistics_ , 7, 2005.

## Checklist

  1. 1.

For all authors…

    1. (a)

Do the main claims made in the abstract and introduction accurately reflect
the paper’s contributions and scope? [Yes]

    2. (b)

Did you describe the limitations of your work? [Yes]

    3. (c)

Did you discuss any potential negative societal impacts of your work? [N/A]

    4. (d)

Have you read the ethics review guidelines and ensured that your paper
conforms to them? [Yes]

  2. 2.

If you are including theoretical results…

    1. (a)

Did you state the full set of assumptions of all theoretical results? [Yes]

    2. (b)

Did you include complete proofs of all theoretical results? [N/A]

  3. 3.

If you ran experiments…

    1. (a)

Did you include the code, data, and instructions needed to reproduce the main
experimental results (either in the supplemental material or as a URL)? [Yes]

    2. (b)

Did you specify all the training details (e.g., data splits, hyper-parameters,
how they were chosen)? [Yes]

    3. (c)

Did you report error bars (e.g., with respect to the random seed after running
experiments multiple times)? [N/A]

    4. (d)

Did you include the total amount of compute and the type of resources used
(e.g., type of GPUs, internal cluster, or cloud provider)? [Yes]

  4. 4.

If you are using existing assets (e.g., code, data, models) or
curating/releasing new assets…

    1. (a)

If your work uses existing assets, did you cite the creators? [Yes]

    2. (b)

Did you mention the license of the assets? [Yes]

    3. (c)

Did you include any new assets either in the supplemental material or as a
URL? [Yes]

    4. (d)

Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [Yes]

    5. (e)

Did you discuss whether the data you are using/curating contains personally
identifiable information or offensive content? [Yes]

  5. 5.

If you used crowdsourcing or conducted research with human subjects…

    1. (a)

Did you include the full text of instructions given to participants and
screenshots, if applicable? [N/A]

    2. (b)

Did you describe any potential participant risks, with links to Institutional
Review Board (IRB) approvals, if applicable? [N/A]

    3. (c)

Did you include the estimated hourly wage paid to participants and the total
amount spent on participant compensation? [N/A]

##  Appendix A Appendix

###  A.1 datasets

The URL to all datasets used throughout this study are presented in Table 4.

Table 4: List of datasets used in this study dataset | URL
---|---
Adult | https://archive.ics.uci.edu/ml/datasets/adult
Attrition | https://www.kaggle.com/philschmidt/employee-attrition-eda
Audit | https://www.kaggle.com/sid321axn/audit-data
Banking | https://www.kaggle.com/rashmiranu/banking-dataset-classification
Banknote | https://archive.ics.uci.edu/ml/datasets/banknote+authentication
Breast Cancer | https://archive.ics.uci.edu/ml/datasets/breast+cancer+wisconsin+(diagnostic)
Churn | https://www.kaggle.com/surendharanp/personal-loan
Donors | https://www.kaggle.com/momohmustapha/donorsprediction
Pima Indians | https://www.kaggle.com/uciml/pima-indians-diabetes-database
Haberman | https://archive.ics.uci.edu/ml/datasets/haberman’s+survival
Hattrick | https://www.kaggle.com/juandelacalle/hattirckorg-matches-dataset
Heart Disease | https://archive.ics.uci.edu/ml/datasets/heart+disease
HR | https://www.kaggle.com/sid321axn/audit-data
Insurance | https://www.kaggle.com/mhdzahier/travel-insurance
Iris | https://archive.ics.uci.edu/ml/datasets/iris
Loan | https://www.kaggle.com/teertha/personal-loan-modeling
Seismic | https://archive.ics.uci.edu/ml/datasets/seismic-bumps
Spambase | https://archive.ics.uci.edu/ml/datasets/spambase
Thera | https://www.kaggle.com/surendharanp/personal-loan
Titanic | https://www.kaggle.com/c/titanic

###  A.2 Experiments

The accuracy of the Logistic Regression and Naive Bayes models are shown in
Table 5. The rank of each explanation techniques according to median
Spearman’s rank correlation in a descending order across different datasets
and underlying models are depicted in Table 6.

Dataset LR NB Dataset LR NB Dataset LR NB Adult 0.841 0.805 Donor 1 0.99 Iris
1 1 Attrition 1 1 Pima Indians 0.802 0.766 Loan 1 1 Audit 0.99 0.974 Haberman
0.662 0.649 Seismic 0.946 0.837 Banking 0.906 0.879 Hattrick 1 0.939 Spambase
0.921 0.808 Banknote 0.98 0.854 Heart Disease 0.829 0.829 Thera 1 1 Breast
Cancer 0.965 0.916 HR 0.77 0.749 Titanic 0.794 0.78 Churn 0.805 0.826
Insurance 0.987 0.953

Table 5: Accuracy of Logistic Regression (LR) and Naive Bayes (NB) models with normalization on all datasets Table 6: Rank of explanation techniques based on median Spearman’s rank correlation within each dataset for Logistic Regression and Naive Bayes models | Logistic Regression | Naive Bayes
---|---|---
| LIME | SHAP | LPI | LIME | SHAP | LPI
Adult | 1 | 3 | 2 | 1 | 3 | 2
Attrition | 1 | 3 | 2 | 2 | 3 | 1
Audit | 2 | 3 | 1 | 2 | 3 | 1
Banking | 2 | 1 | 3 | 2 | 3 | 1
Banknote | 3 | 1 | 2 | 3 | 2 | 1
Breast Cancer | 3 | 1 | 2 | 3 | 2 | 1
Churn | 1 | 2 | 3 | 1.5 | 3 | 1.5
Donors | 3 | 2 | 1 | 2 | 3 | 1
HR | 2 | 3 | 1 | 1 | 3 | 2
Haberman | 3 | 1.5 | 1.5 | 3 | 2 | 1
Hattrick | 3 | 1 | 2 | 1 | 2 | 3
Heart Disease | 3 | 1 | 2 | 1 | 3 | 2
Insurance | 1 | 3 | 2 | 1.5 | 3 | 1.5
Iris | 3 | 2 | 1 | 3 | 2 | 1
Loan | 1 | 2 | 3 | 2 | 3 | 1
Pima Indians | 3 | 1 | 2 | 3 | 2 | 1
Seismic | 1 | 3 | 2 | 2 | 3 | 1
Spambase | 3 | 2 | 1 | 2 | 3 | 1
Thera | 1 | 3 | 2 | 2 | 3 | 1
Titanic | 2 | 3 | 1 | 2 | 3 | 1
Average (All data sets) | 2.1 | 2.08 | 1.82 | 2 | 2.7 | 1.3
Standard Deviation (All data sets) | 0.89 | 0.84 | 0.66 | 0.69 | 0.46 | 0.53

###  A.3 Pre-processing Effect

####  A.3.1 Min-max scaling

In Table 7, we present the rank of explanation techniques based on median
Spearman’s correlation for each dataset when min-max scaling is performed in
data sets.

Table 7: Rank of explanation techniques based on median Spearman rank correlation when decomposing log odds ratio under the effect of minmax pre-processing | Logistic Regression | Naive Bayes
---|---|---
| LIME | SHAP | LPI | LIME | SHAP | LPI
Adult | 1 | 3 | 2 | 1 | 3 | 2
Attrition | 2 | 1 | 3 | 2 | 3 | 1
Audit | 2 | 3 | 1 | 1 | 3 | 2
Banking | 2 | 1 | 3 | 2 | 3 | 1
Banknote | 1 | 2.5 | 2.5 | 3 | 2 | 1
Breast Cancer | 1 | 2 | 3 | 3 | 2 | 1
Churn | 1 | 2 | 3 | 1 | 3 | 2
Donors | 1 | 2 | 3 | 2 | 3 | 1
HR | 2 | 3 | 1 | 1 | 3 | 2
Haberman | 1.5 | 1.5 | 3 | 3 | 2 | 1
Hattrick | 2 | 1 | 3 | 1 | 2 | 3
Heart Disease | 1 | 2 | 3 | 1 | 3 | 2
Insurance | 1 | 3 | 2 | 1.5 | 3 | 1.5
Iris | 1 | 2.5 | 2.5 | 3 | 2 | 1
Loan | 2 | 1 | 3 | 2 | 3 | 1
Pima Indians | 1 | 2 | 3 | 3 | 2 | 1
Seismic | 2 | 1 | 3 | 2 | 3 | 1
Spambase | 2 | 1 | 3 | 3 | 2 | 1
Thera | 2 | 1 | 3 | 2 | 3 | 1
Titanic | 2 | 3 | 1 | 1 | 3 | 2
Average (All data sets) | 1.52 | 1.92 | 2.55 | 1.92 | 2.65 | 1.42
Standard Deviation (All data sets) | 0.49 | 0.79 | 0.72 | 0.81 | 0.48 | 0.58

####  A.3.2 Inter-quartile normalization

In Table 7, we present the rank of explanation techniques based on median
Spearman’s correlation for each dataset when inter-quartile pre-processing is
performed in data sets.

Table 8: Rank of explanation techniques under the inter-quartile normalization | Logistic Regression | Naive Bayes
---|---|---
| LIME | SHAP | LPI | LIME | SHAP | LPI
Adult | 1 | 2 | 3 | 1 | 3 | 2
Attrition | 1 | 2 | 3 | 2 | 3 | 1
Audit | 2 | 1 | 3 | 1 | 3 | 2
Banking | 2 | 1 | 3 | 2 | 3 | 1
Banknote | 3 | 1 | 2 | 3 | 2 | 1
Breast Cancer | 3 | 1 | 2 | 3 | 2 | 1
Churn | 2 | 1 | 3 | 1 | 3 | 2
Donors | 3 | 2 | 1 | 2 | 3 | 1
HR | 3 | 2 | 1 | 1 | 3 | 2
Haberman | 2.5 | 1 | 2.5 | 3 | 2 | 1
Hattrick | 2 | 1 | 3 | 1 | 2 | 3
Heart Disease | 3 | 1 | 2 | 1 | 3 | 2
Insurance | 1 | 3 | 2 | 1.5 | 3 | 1.5
Iris | 3 | 1.5 | 1.5 | 3 | 2 | 1
Loan | 2 | 1 | 3 | 2 | 3 | 1
Pima Indians | 3 | 1 | 2 | 3 | 2 | 1
Seismic | 2 | 1 | 3 | 2 | 3 | 1
Spambase | 2 | 1 | 3 | 3 | 2 | 1
Thera | 2 | 1 | 3 | 2 | 3 | 1
Titanic | 3 | 2 | 1 | 2 | 3 | 1
Average (All data sets) | 2.28 | 1.38 | 2.35 | 1.98 | 2.65 | 1.38
Standard Deviation (All data sets) | 0.7 | 0.57 | 0.74 | 0.78 | 0.48 | 0.57


Also at ]Laboratory for Manufacturing Systems and Automation, Department of
Mechanical Engineering and Aeronautics, University of Patras, 26110 Patra,
Greece

# Entanglement Entropy as a Probe Beyond the Horizon

Konstantinos Boutivas<EMAIL_ADDRESS>National and Kapodistrian
University of Athens, Department of Physics, 15784 Zografou, Attiki, Greece
Dimitrios Katsinis<EMAIL_ADDRESS>National and Kapodistrian University
of Athens, Department of Physics, 15784 Zografou, Attiki, Greece  Georgios
Pastras  [<EMAIL_ADDRESS>National and Kapodistrian University
of Athens, Department of Physics, 15784 Zografou, Attiki, Greece  Nikolaos
Tetradis<EMAIL_ADDRESS>National and Kapodistrian University of Athens,
Department of Physics, 15784 Zografou, Attiki, Greece

###### Abstract

The entanglement entropy of a free field in de Sitter space is enhanced by the
squeezing of its modes. We show analytically that the expansion induces a term
in the entanglement entropy that depends logarithmically on the size of the
overall system, which may extend beyond the horizon. In cosmology the size of
the system can be identified with the size of a spatially finite universe, or
with the wavelength of the first mode that exited the horizon in the beginning
of inflation.

##  I Introduction

According to the standard inflationary paradigm, the quantum fluctuations of
massless fields become classical upon horizon exit. Upon re-entry after the
end of inflation, they result in density perturbations that generate the
observed large-scale structure of the universe [1, 2, 3, 4, 5]. This quantum-
to-classical transition can be attributed to the fact that the field mode
functions consist of a growing and a decaying term, with the growing one
becoming dominant after horizon exit. If the decaying term is neglected, the
field commutes with its conjugate momentum. Therefore, it can be treated as a
classical stochastic field. This simple line of reasoning does not take into
account the fact that, no matter how suppressed the decaying mode is, the
field and its conjugate momentum always obey the canonical commutation
relation. Strictly speaking, the field never loses its quantum nature, even
though measuring quantum effects may be extremely difficult [6, 7, 8, 9, 10,
11].

It is known that the field modes evolve from a simple harmonic oscillator
ground state to an increasingly squeezed one, which ceases to be a state of
minimal uncertainty [6, 12]. From this point of view, it seems possible that
the quantum nature of the fluctuations may be discernible, as long as one
looks at the appropriate quantity. There is a connection between the squeezing
of the field modes in a cosmological setting and the production of entropy
[13, 14, 15, 16, 17, 18]. We focus here on entanglement entropy, which is a
quantity that reflects directly the quantum nature of the fluctuations. We
work in coordinate space, in an expanding background that includes a
cosmological horizon. We compute the entanglement entropy of a subsystem
enclosed by an entangling surface. We show that the expansion enhances the
effects of entanglement, leading to an increase in entropy, and also inducing
its dependence on global properties of space. For a de Sitter background, the
leading correction to the entanglement entropy beyond its flat-space value
depends on the size of the _entire system_ , which may extend beyond the
horizon.

The archetypical calculations in $(3+1)$-dimensional flat space [19, 20] led
to the unexpected finding that the leading contribution to the entanglement
entropy of a free massless scalar field is proportional to the area of the
surface enclosing the subsystem under consideration. This result is attributed
to the strong entanglement of short-distance modes on either side of the
entangling surface, which also makes the entropy UV-divergent. The numerical
calculation of the coefficient of the area-law term is a relatively easy task
nowadays, given the available computing power. However, the computation of the
coefficients of the subleading terms, such as the universal logarithmic term,
requires an elaborate numerical study [21]. Analytical calculations, even for
a free field theory, are more difficult [22]. There are only a few analytical
results, mainly in lower dimensions [23, 24, 25, 26, 27]. The effects of the
expansion have been analyzed even less [28, 29, 30, 31, 32, 33, 34].

We base our analysis on the approach of Srednicki in its seminal work [19]. We
calculate the entanglement entropy of a massless scalar field in the Bunch-
Davies vacuum of $(3+1)$-dimensional de Sitter space. The necessary formalism
has been developed in our recent work [33, 34] that uses the generalization of
Srednicki’s method for harmonic systems at squeezed states [35, 36]. We
consider a spherical entangling surface, expand the field in spherical
moments, and discretize the radial direction. The theory is decomposed into an
infinite number of independent sectors, each of them corresponding to a fixed
value of angular momentum $\ell$. Each sector amounts to a discretized
effective $(1+1)$-dimensional system, whose size is set by the number of
points $N$ of the radial lattice. The entanglement entropy is calculated as
the sum of the contributions of all sectors, namely

$S=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}S_{\ell
m}=\sum_{\ell=0}^{\infty}\left(2\ell+1\right)S_{\ell},$ (1)

taking into account that the dynamics of each sector depends only on $\ell$.

A crucial step in the calculation is taking the limit $N\sim
L/\epsilon\to\infty$, where $L$ is the size of the system and $\epsilon$ the
UV cutoff, identified with the lattice spacing. The main obstacle is that
subleading terms in the entropy must be disentangled from finite-size effects
that inevitably appear in the calculation. This problem is ubiquitous in
numerical calculations of the entropy. A method for dealing with the finite-
size effects in flat space has been developed in [21]. An important
observation is that these effects are suppressed for large angular momentum.
The higher the angular momentum, the more disentangled the local degrees of
freedom of the theory are. This happens because the angular momentum acts as a
position-dependent mass, making the normal modes more and more localized
around certain points as its value increases. In $(3+1)$-dimensional flat
space it turns out that the finite-size effects scale as [21]

$S^{\textrm{flat}}_{\textrm{fs}}\propto N^{-2\left(\ell+1\right)}.$ (2)

It is evident that these terms vanish in the infinite-size limit.

In an accompanying publication [37], we perform a thorough numerical study of
entanglement entropy in de Sitter space, modifying appropriately the
methodology developed for flat space in [21]. We consider a field in the
Bunch-Davies vacuum and observe how entanglement is enhanced through the
squeezing of the state. It turns out that the dominant finite-size effects now
scale as

$S^{\textrm{dS}}_{\textrm{fs}}\propto\begin{cases}\ln N,&\ell=0,\\\
N^{-2\ell},&\ell\neq 0.\end{cases}$ (3)

Rather surprisingly, the contribution of the $\ell=0$ sector does not vanish
in the infinite-size limit. The numerical calculation indicates that the
resulting contribution to the entropy is of the novel form

$S^{\textrm{dS}}_{\textrm{IR}}\propto\frac{R^{2}}{\tau^{2}}\ln\frac{L}{\epsilon},$
(4)

where $L$ is the radius of the _overall system_ , $R$ the radius of the
subsystem, $\tau$ the conformal time and $\epsilon$ the UV cutoff. The
coefficient of this term is found numerically to be approximately equal to
$1/3$. In this work we verify this fact with an analytical calculation.

The means to perform this calculation are provided by our recent work [38], in
which the $(1+1)$-dimensional discretized theory is studied analytically in
order to obtain the continuum limit of the entanglement spectrum, the
entanglement entropy and the corresponding modular Hamiltonian. The main idea
is to take the continuum limit of the eigenvalue equations of the discrete
problem, in which matrices become kernels and matrix multiplication turns into
integration. As we discussed above, the $\ln N$ term stems solely from the
sector of the theory with vanishing angular momentum. This sector is identical
to the $(1+1)$-dimensional free scalar field theory, so that the methodology
and some results of [38] are directly applicable to our calculation here.

##  II The Discretized Theory

We consider a free massless scalar field $\phi\left(\tau,\mathbf{x}\right)$ in
$(3+1)$-dimensional de Sitter space. Performing the field redefinition
$\phi\left(\tau,\mathbf{x}\right)=f\left(\tau,\mathbf{x}\right)/a\left(\tau\right)$,
where the scale factor is $a(\tau)=-1/(H\tau)$ and $H$ is the Hubble constant,
yields the action

$\mathcal{S}=\frac{1}{2}\int
d\tau\,d^{3}\mathbf{x}\,\left(\dot{f}^{2}-\left(\nabla
f\right)^{2}+\frac{\ddot{a}}{a}f^{2}\right),$ (5)

with a canonical kinetic term. The field $f\left(\tau,\mathbf{x}\right)$
obtains a time-dependent squared mass $-\ddot{a}/a=-2/\tau^{2}$. We consider a
spherical region of finite radius $L$ and expand the field in its spherical
moments $f_{\ell m}(r)$. The Hamiltonian corresponding to each
$(\ell,m)$-sector of the theory reads [39]

$H_{\ell m}=\frac{1}{2}\int_{0}^{L}dr\bigg{[}\pi_{\ell
m}^{2}+\left(\partial_{r}f_{\ell m}\right)^{2}\\\
+\left(\frac{\ell\left(\ell+1\right)}{r^{2}}-\frac{2}{\tau^{2}}\right)f_{\ell
m}^{2}\bigg{]}-\frac{1}{2}\frac{f_{\ell m}^{2}}{r}\bigg{|}^{L}_{0}.$ (6)

Considering Dirichlet boundary conditions we can drop the boundary term. Then
we discretize the radial coordinate as $r=j\epsilon$, with $L=(N+1)\epsilon$,
obtaining

$H_{\ell m}=\frac{1}{2\epsilon}\sum\limits_{j=1}^{N}\bigg{[}\pi_{\ell
m,j}^{2}+\left(f_{\ell m,j+1}-f_{\ell m,j}\right)^{2}\\\
+\left(\frac{\ell\left({\ell+1}\right)}{j^{2}}-\frac{2\epsilon^{2}}{\tau^{2}}\right)f_{\ell
m,j}^{2}\bigg{]}.$ (7)

The discretization we used here is slightly different, but completely
equivalent, to the one used by Srednicki [19]. The conformal time $\tau$ is
naturally measured in units of $\epsilon$, which we set equal to 1.

The system lies in the Bunch-Davies vacuum, which implies that the wave
functions of its modes for $\tau\to-\infty$ must reduce to those in flat
space. As a result, each mode of the system is described by a wave function
[33, 34]

$F(\tau,f)=\frac{1}{\sqrt{b(\tau)}}\,\exp\left({\frac{i}{2}\frac{b^{\prime}(\tau)}{b(\tau)}f^{2}}\right)\\\
\times F^{0}\left(\int\frac{d\tau}{b^{2}(\tau)},\frac{f}{b(\tau)}\right)$ (8)

where $F^{0}(\tau,f)$ is a solution of the standard simple harmonic oscillator
with constant frequency $\omega_{0}$. The function
$b(\tau)=\sqrt{1+1/(\omega^{2}_{0}\tau^{2})}$ quantifies the effect of
squeezing. It appoaches 1 for $\tau\to-\infty$, while it diverges for $\tau\to
0^{-}$. For the ground state we obtain

$F\left(\tau,x\right)=\sqrt[4]{\frac{\omega_{0}}{\pi}\frac{1}{1+\frac{1}{\omega_{0}^{2}\tau^{2}}}}\exp\left[-\frac{\omega_{0}}{2}\frac{\omega_{0}^{2}\tau^{2}+\frac{i}{\omega_{0}\tau}}{1+\omega_{0}^{2}\tau^{2}}x^{2}\right]\\\
\times\exp\left[-\frac{i}{2}\left(\omega_{0}\tau-\arctan\left(\omega_{0}\tau\right)-\frac{\pi}{2}\right)\right].$
(9)

We focus on the $\ell=0$ sector of the theory, because the effect we are
interested in originates in this sector only. The wave function describing the
degrees of freedom of this sector reads

$\Psi\left(\mathbf{x}\right)\propto\exp\left(-\frac{1}{2}x^{T}Wx\right),$ (10)

where the matrix $W$ has the form

$W=\Omega\left[I-\left(I-\frac{i}{\tau}\Omega^{-1}\right)\left(I+\Omega^{2}\tau^{2}\right)^{-1}\right].$
(11)

The matrix $\Omega$ is the positive square root of the flat-space coupling
matrix $K$, which can be read directly from the Hamiltonian (7). It is given
by

$K_{ij}=2\delta_{i,j}-\delta_{i,j+1}-\delta_{i+1,j}.$ (12)

We consider the subsystem $A$ consisting of the first $n$ oscillators, and we
trace out the remaining $N-n$ ones that constitute subsystem $A^{C}$.
Throughout the paper we denote the blocks of $N\times N$ matrices as
$\Omega=\begin{pmatrix}\Omega_{A}&\Omega_{B}\\\
\Omega_{B}^{T}&\Omega_{C}\end{pmatrix}$. In particular, the matrix $M_{A}$ is
$n\times n$, whereas the matrix $M_{C}$ is $(N-n)\times(N-n)$. The
entanglement entropy of Gaussian states, like the one in Eq. (10), can be
calculated by making use of correlation functions [40]. We define the
$2N\times 2N$ matrix $\mathcal{M}=2iJ\mathrm{Re}\,\left(M\right)$, where

$M=\begin{pmatrix}\left\langle x_{i}x_{j}\right\rangle&\left\langle
x_{i}\pi_{j}\right\rangle\\\ \left\langle
x_{i}\pi_{j}\right\rangle^{T}&\left\langle\pi_{i}\pi_{j}\right\rangle\end{pmatrix},\qquad
J=\begin{pmatrix}0&I\\\ -I&0\end{pmatrix}.$ (13)

The $2n\times 2n$ matrix $\mathcal{M}$ for subsystem $A$ has a relatively
simple form in terms of the blocks of half-integer powers of the couplings
matrix $K$ [37], namely

$\mathcal{M}_{A}=i\begin{pmatrix}\frac{1}{\tau^{3}}\left(\Omega^{-3}\right)_{A}&\left(\Omega\right)_{A}-\frac{1}{\tau^{2}}\left(\Omega^{-1}\right)_{A}+\frac{1}{\tau^{4}}\left(\Omega^{-3}\right)_{A}\\\
-\left(\Omega^{-1}\right)_{A}-\frac{1}{\tau^{2}}\left(\Omega^{-3}\right)_{A}&-\frac{1}{\tau^{3}}\left(\Omega^{-3}\right)_{A}\end{pmatrix}.$
(14)

Its eigenvalues come in pairs $\pm\lambda_{i}$ and satisfy
$\left|\lambda_{i}\right|\geq 1$. The entanglement entropy can be computed as

$S=\sum_{i=0}^{n}\left(\frac{\lambda_{i}+1}{2}\ln\frac{\lambda_{i}+1}{2}-\frac{\lambda_{i}-1}{2}\ln\frac{\lambda_{i}-1}{2}\right),$
(15)

where the sum is performed only over the positive eigenvalues.

Calculating analytically the eigenvalues of (14) for an arbitrary value of
$\tau$ is extremely hard. However, it is possible to implement perturbation
theory for large $|\tau|$ in order to calculate the leading corrections to the
entanglement entropy of flat space. The details of the perturbative expansion
are given in [37]. The calculation is facilitated if we consider the square of
the matrix $\mathcal{M}_{A}$, along with a similarity transformation,
obtaining

$\tilde{\mathcal{M}}^{2}_{A}=\begin{pmatrix}\frac{\mathcal{M}^{(0)T}}{\tau^{0}}+\frac{\mathcal{M}^{(-2)T}}{\tau^{2}}+\frac{\mathcal{M}^{(-4)}}{\tau^{4}}&\frac{\mathcal{M}^{(-2)}-\mathcal{M}^{(-2)T}}{\tau^{2}}-\frac{\mathcal{M}^{(-4)}}{\tau^{4}}\\\
\frac{\mathcal{M}^{(-4)}}{\tau^{4}}&\frac{\mathcal{M}^{(0)}}{\tau^{0}}+\frac{\mathcal{M}^{(-2)}}{\tau^{2}}-\frac{\mathcal{M}^{(-4)}}{\tau^{4}}\end{pmatrix},$
(16)

where

$\displaystyle\mathcal{M}^{(0)}\thickspace$
$\displaystyle=\left(\Omega^{-1}\right)_{A}\big{(}\Omega\big{)}_{A},$ (17)
$\displaystyle\mathcal{M}^{(-2)}$
$\displaystyle=\left(\Omega^{-3}\right)_{A}\big{(}\Omega\big{)}_{A}-\left(\Omega^{-1}\right)_{A}^{2},$
(18) $\displaystyle\mathcal{M}^{(-4)}$
$\displaystyle=\left(\Omega^{-3}\right)_{A}\left(\Omega^{-1}\right)_{A}-\left(\Omega^{-1}\right)_{A}\left(\Omega^{-3}\right)_{A}.$
(19)

The eigenvalues of $\tilde{\mathcal{M}}^{2}_{A}$, which are the squares of the
eigenvalues of $\mathcal{M}_{A}$, depend only on even powers of $1/|\tau|$.
The eigenvalues of $\mathcal{M}_{A}$ admit a perturbative expansion of the
form

$\lambda_{i}=\sum_{k=0}^{\infty}\frac{\lambda^{(2k)}_{i}}{\tau^{2k}}.$ (20)

In the following we use the same notation for similar expansions.

In the flat-space limit we have the eigenvalue problem

$\mathcal{M}^{(0)T}v^{(0)}_{i}=\Lambda_{i}^{(0)}v^{(0)}_{i},\qquad\mathcal{M}^{(0)}w^{(0)}_{i}=\Lambda_{i}^{(0)}w^{(0)}_{i},$
(21)

where $\Lambda_{i}^{(0)}=\left(\lambda^{(0)}_{i}\right)^{2}$. The vectors
$v^{(0)}_{i}$ and $w^{(0)}_{i}$ are the left and right eigenvectors of
$\mathcal{M}^{(0)}$. The leading corrections to the flat-space eigenvalues
$\lambda^{(0)}_{i}$ read

$\lambda^{(2)}_{i}=\frac{1}{2\lambda^{(0)}_{i}}\frac{w^{(0)T}_{i}\mathcal{M}^{(-2)T}v^{(0)}_{i}}{w^{(0)T}_{i}v^{(0)}_{i}}.$
(22)

For the whole system we have
$\Omega\,\Omega^{-3}=\Omega^{-2}=\Omega^{-1}\Omega^{-1}$ along with
$\Omega^{-1}\Omega=I_{N}$. This allows us to write

$\displaystyle\mathcal{M}^{(0)}$
$\displaystyle=I_{n}-\left(\Omega^{-1}\right)_{B}\big{(}\Omega\big{)}_{B}^{T},$
(23) $\displaystyle\mathcal{M}^{(-2)T}$
$\displaystyle=\left(\Omega^{-1}\right)_{B}\left(\Omega^{-1}\right)_{B}^{T}-\big{(}\Omega\big{)}_{B}\left(\Omega^{-3}\right)_{B}^{T}.$
(24)

The form of Eqs. (23) and (24) is advantageous because we avoid UV
singularities that arise in the continuum limit (when, for instance, we
integrate over $y$ in an expression such as Eq. (40) and $y\rightarrow x$ or
$y\rightarrow x^{\prime}$).

##  III The Continuum Limit

###  III.1 The Kernels

We are interested in the eigenvalues $\lambda_{i}^{(0)}$ and the corrections
$\lambda_{i}^{(2)}$ in the continuum limit. Restoring momentarily the cutoff
$\epsilon$ implies that the terms $\lambda_{i}^{(2k)}$ of Eq. (20) are
actually $\epsilon^{2k}\lambda_{i}^{(2k)}$. As a result, the kernels
$\Omega^{m}\left(x,x^{\prime}\right)$ corresponding to the matrices
$\Omega^{m}$ are defined as

$\Omega^{m}\left(x,x^{\prime}\right)=\lim_{\epsilon\rightarrow
0}\frac{\Omega^{m}_{ij}}{\epsilon^{m}},$ (25)

where $x=i\epsilon$ and $x^{\prime}=j\epsilon$. The $A$-block of the matrix
$\Omega^{m}$ corresponds to the kernel $\Omega^{m}\left(x,x^{\prime}\right)$
with the restriction $x\leq R$ and $x^{\prime}\leq R$.

The spectral decomposition of the coupling matrix of Eq. (12) is known and one
can take the continuum limit to obtain the eigenfunctions and the
corresponding eigenvalues [38]. We find

$\Omega^{m}\left(x,x^{\prime}\right)=\frac{2}{L}\sum_{k=0}^{\infty}\left(\frac{k\pi}{L}\right)^{m}\sin\frac{k\pi
x}{L}\sin\frac{k\pi x^{\prime}}{L},$ (26)

where $x$ and $x^{\prime}$ are valued in $\left[0,L\right]$.

The kernels that are relevant for our calculation can be written in closed
form as

$\displaystyle\Omega(x,x^{\prime})$
$\displaystyle=\frac{\pi}{4L^{2}}\left(\frac{1}{\sin^{2}\frac{\pi(x+x^{\prime})}{2L}}-\frac{1}{\sin^{2}\frac{\pi(x-x^{\prime})}{2L}}\right),$
(27) $\displaystyle\Omega^{-1}(x,x^{\prime})$
$\displaystyle=-\frac{1}{\pi}\ln\left|\frac{\sin\frac{\pi\left(x-x^{\prime}\right)}{2L}}{\sin\frac{\pi\left(x+x^{\prime}\right)}{2L}}\right|,$
(28)
$\displaystyle\begin{split}\Omega^{-3}(x,x^{\prime})&=\frac{L^{2}}{2\pi^{3}}\bigg{[}\mathrm{Li}_{3}\left(e^{-i\frac{\pi}{L}\left(x-x^{\prime}\right)}\right)\\\
&\qquad\qquad-\mathrm{Li}_{3}\left(e^{-i\frac{\pi}{L}\left(x+x^{\prime}\right)}\right)\bigg{]}+\textrm{c.c.}\thickspace.\end{split}$
(29)

The regime of interest is $L\gg R$, where $R$ is the radius of the entangling
surface. In this regime the kernels assume the form

$\displaystyle\Omega(x,x^{\prime})$
$\displaystyle=\frac{1}{\pi}\left(\frac{1}{(x+x^{\prime})^{2}}-\frac{1}{(x-x^{\prime})^{2}}\right),$
(30) $\displaystyle\Omega^{-1}(x,x^{\prime})$
$\displaystyle=-\frac{1}{\pi}\ln\left|\frac{x-x^{\prime}}{x+x^{\prime}}\right|,$
(31)
$\displaystyle\begin{split}\Omega^{-3}(x,x^{\prime})&=\frac{1}{\pi}\Bigg{[}xx^{\prime}\left(3+2\ln\frac{L}{\pi
R}\right)\\\
&+\frac{\left(x-x^{\prime}\right)^{2}}{2}\ln\frac{\left|x-x^{\prime}\right|}{R}-\frac{\left(x+x^{\prime}\right)^{2}}{2}\ln\frac{x+x^{\prime}}{R}\Bigg{]}.\end{split}$
(32)

The solution of the eigenvalue problem (21) is known [23, 38]. The right
eigenfunctions $f(x;\omega)$, the left eigenfunctions $g(x;\omega)$, and the
eigenvalues $\Lambda^{(0)}(\omega)$ of the kernel $\mathcal{M}^{(0)}$ read:

$\displaystyle f(x;\omega)$ $\displaystyle=\sin\left(\omega u(x)\right),$ (33)
$\displaystyle g(x;\omega)$
$\displaystyle=\frac{1}{R}\cosh^{2}\frac{u}{2}\sin\left(\omega u(x)\right),$
(34) $\displaystyle\Lambda^{(0)}(\omega)$
$\displaystyle=\coth^{2}\left(\pi\omega\right),$ (35)

where $x$ is related to the variable $u$ as

$x(u)=R\tanh\frac{u}{2},\qquad u(x)=\ln\frac{R+x}{R-x}.$ (36)

We need the inner product $v^{(0)T}_{i}w^{(0)}_{i}$, which is divergent in the
continuum limit. We impose a cutoff at $x=R-\epsilon$, where $\epsilon\ll R$.
Thus, we regularize this divergence and at the same time we discretize the
spectrum as

$\omega_{k}=\frac{k\pi}{u_{\textrm{max}}},\qquad
u_{\textrm{max}}=\ln\frac{2R}{\epsilon},\qquad k\in\mathbb{N}^{*}.$ (37)

It then follows that the continuum limit of the inner product
$w^{(0)T}_{i}v^{(0)}_{i}$ is given by

$w^{(0)T}_{i}v^{(0)}_{i}\rightarrow\frac{1}{2}\int_{0}^{u_{\textrm{max}}}du\sin^{2}\left(\omega_{i}u\right)=\frac{1}{4}u_{\textrm{max}}.$
(38)

###  III.2 The Correction to the Eigenvalues

We next proceed with the main part of our derivation, which is the calculation
of

$\left\langle M^{(-2)T}\right\rangle\\\
=\int_{0}^{R}dx\int_{0}^{R}dx^{\prime}f(x;\omega)M^{(-2)T}\left(x,x^{\prime}\right)g(x^{\prime};\omega).$
(39)

The kernel $M^{(-2)T}\left(x,x^{\prime}\right)$, which is the continuum limit
of the matrix $M^{(-2)T}$ in the form of Eq. (24), reads

$M^{(-2)T}\left(x,x^{\prime}\right)=\int^{L}_{\ell}dy\left[\Omega^{-1}(x,y)\Omega^{-1}(y,x^{\prime})\right.\\\
\left.-\Omega(x,y)\Omega^{-3}(y,x^{\prime})\right].$ (40)

Taking into account that

$\Omega^{m}(x,y)=-\frac{d^{2}}{dy^{2}}\Omega^{m-2}(x,y),$ (41)

which is a direct consequence of Eq. (26), we obtain

$M^{(-2)T}\left(x,x^{\prime}\right)=\Omega^{-1}(x,R)\left[\frac{\partial}{\partial
R}\Omega^{-3}(R,x^{\prime})\right]\\\ -\left[\frac{\partial}{\partial
R}\Omega^{-1}(x,R)\right]\Omega^{-3}(R,x^{\prime}).$ (42)

Notice that $M^{(-2)T}\left(x,x^{\prime}\right)$ is a sum of products of
functions of either $x$ or $x^{\prime}$. As a result, it is finite when
$x^{\prime}\rightarrow x$. Divergences occur only when $x\rightarrow R$ or
$x^{\prime}\rightarrow R$. For $L\gg R$, using Eqs. (31) and (32), we obtain

$M^{(-2)T}\left(x,x^{\prime}\right)=\frac{1}{\pi^{2}}\left[x^{\prime}+\frac{R^{2}-x^{\prime
2}}{2R}\ln\frac{R+x^{\prime}}{R-x^{\prime}}\right]\frac{2Rx}{R^{2}-x^{2}}\\\
+\frac{1}{\pi^{2}}\left[\frac{2Rx}{R^{2}-x^{2}}+\ln\frac{R+x}{R-x}\right]\left[2x^{\prime}\left(\ln\frac{L}{\pi
R}+1\right)\right.\\\
\left.+\left(R-x^{\prime}\right)\ln\frac{R-x^{\prime}}{R}-\left(R+x^{\prime}\right)\ln\frac{R+x^{\prime}}{R}\right].$
(43)

Substituting the eigenfunctions, given by Eqs. (33) and (34), and performing
the change of variable (36), one obtains $\left\langle
M^{(-2)T}\right\rangle$. Making use of the normalization (38), the corrections
of the eigenvalues $\lambda^{(2)}_{i}$, defined in Eq. (22), assume the form

$\lambda^{(2)}\left(\omega_{i}\right)=R^{2}\frac{\left[I_{1}\left(\omega_{i}\right)\left(3+2\ln\frac{L}{\pi
R}\right)+I_{2}\left(\omega_{i}\right)\right]}{u_{\textrm{max}}\coth\left(\pi\omega_{i}\right)}\\\
\times\left[I_{1}\left(\omega_{i}\right)+\frac{1}{2}I_{3}\left(\omega_{i}\right)\right],$
(44)

where

$\displaystyle I_{1}\left(\omega\right)$
$\displaystyle=\frac{1}{\pi}\int_{0}^{u_{\textrm{max}}}du\sin\left(\omega
u\right)\tanh\frac{u}{2},$ (45)
$\displaystyle\begin{split}I_{2}\left(\omega\right)&=\frac{1}{\pi}\int_{0}^{u_{\textrm{max}}}du\sin\left(\omega
u\right)\bigg{(}\frac{e^{-u/2}}{\cosh\frac{u}{2}}\ln\frac{e^{-u/2}}{\cosh\frac{u}{2}}\\\
&\qquad\qquad\qquad\qquad\quad\thickspace-\frac{e^{u/2}}{\cosh\frac{u}{2}}\ln\frac{e^{u/2}}{\cosh\frac{u}{2}}\bigg{)},\end{split}$
(46) $\displaystyle I_{3}\left(\omega\right)$
$\displaystyle=\frac{1}{\pi}\int_{0}^{u_{\textrm{max}}}du\,\frac{u\sin\left(\omega
u\right)}{\cosh^{2}\frac{u}{2}}.$ (47)

We are interested in the $u_{\textrm{max}}\rightarrow\infty$ limit of these
integrals, which are calculated as contour integrals. We obtain

$\displaystyle I_{1}\left(\omega\right)$
$\displaystyle=\frac{1}{\sinh\left(\pi\omega\right)},$ (48) $\displaystyle
I_{1}\left(\omega\right)+\frac{1}{2}I_{3}\left(\omega\right)$
$\displaystyle=\pi\omega\frac{\coth\left(\pi\omega\right)}{\sinh\left(\pi\omega\right)}.$
(49)

We do not make use of the explicit expression for $I_{2}\left(\omega\right)$
in the following, thus we do not present it.

###  III.3 The Correction to the Entanglement Entropy

Following the notation of the expansion (20), the first non-vanishing
correction to the entanglement entropy of flat space $S^{(0)}$ is given by

$S^{(2)}=\sum_{i}\lambda^{(2)}_{i}\mathrm{arccoth}\,\lambda^{(0)}_{i}=\sum_{i}\pi\omega_{i}\lambda^{(2)}_{i}.$
(50)

Using Eqs. (44), (48) and (49), we obtain

$S^{(2)}=R^{2}\frac{\pi}{u_{\textrm{max}}}\sum_{i}\left[\left(3+2\ln\frac{L}{\pi
R}\right)\frac{\pi\omega_{i}^{2}}{\sinh^{2}\left(\pi\omega_{i}\right)}\right.\\\
\left.+\frac{\pi\omega_{i}I_{2}\left(\omega_{i}\right)}{\sinh\left(\pi\omega_{i}\right)}\right].$
(51)

Since $R\gg\epsilon$ the eigenvalues $\omega_{i}$ become dense and the sum can
be approximated by an integral over $\omega$ [23] using Eq. (37), i.e.

$S^{(2)}=R^{2}\int_{0}^{\infty}d\omega\left[\left(3+2\ln\frac{L}{\pi
R}\right)\frac{\pi\omega^{2}}{\sinh^{2}\left(\pi\omega\right)}\right.\\\
\left.+\frac{\pi\omega
I_{2}\left(\omega\right)}{\sinh\left(\pi\omega\right)}\right].$ (52)

The integral of the first term can be calculated straightforwardly. Regarding
the second term, we perform first the integral over $\omega$, using

$\int_{0}^{\infty}d\omega\frac{\omega\sin\left(\omega
u\right)}{\sinh\left(\pi\omega\right)}=\frac{\sinh\frac{u}{2}}{4\cosh^{3}\frac{u}{2}}$
(53)

and then we integrate over $u$, obtaining

$\int_{0}^{\infty}du\frac{\sinh\frac{u}{2}}{4\cosh^{3}\frac{u}{2}}\left(\frac{e^{-u/2}}{\cosh\frac{u}{2}}\ln\frac{e^{-u/2}}{\cosh\frac{u}{2}}\right.\\\
\left.-\frac{e^{u/2}}{\cosh\frac{u}{2}}\ln\frac{e^{u/2}}{\cosh\frac{u}{2}}\right)=-\frac{1}{18}-\frac{1}{3}\ln
2.$ (54)

Putting everything together, the correction to the flat-space entanglement
entropy, arising from the $\ell=0$ sector at order $1/\tau^{2}$, reads

$S^{(2)}=\left(\frac{1}{3}\ln\frac{L}{2\pi R}+\frac{4}{9}\right)R^{2}.$ (55)

The higher-$\ell$ sectors give contributions independent of $L$, but with a
logarithmic dependence on $\epsilon$ [37]. The result (55) determines the
$L$-dependence of the total entropy. This means that the precise form of Eq.
(4), including the numerical coefficient, is

$S^{\textrm{dS}}_{\textrm{IR}}=\frac{1}{3}\frac{R^{2}}{\tau^{2}}\ln\frac{L}{\epsilon}.$
(56)

##  IV Discussion

In an accompanying publication [37] we calculate numerically the entanglement
entropy for a free field in $(3+1)$-dimensional de Sitter space. The numerical
study confirms the appearance of the term (56) as the leading correction to
the flat-space entanglement entropy. The numerical coefficient equal to $1/3$
is reproduced with very high accuracy. It is worth mentioning that the
analysis in [37] employs the original discretization scheme of Srednicki [19],
which differs slightly from the one we employed here. This confirms that the
continuum limit is independent of the details of the discretization.

The logarithmic term (56) originates in the sector with vanishing angular
momentum, which is equivalent to a $(1+1)$-dimensional theory. Moreover, it
results from the lowest frequency modes, whose wavelengths are comparable to
the size $L$ of the _overall system_. In flat space, their contribution is
suppressed. However, the cosmological expansion induces the strong squeezing
of these modes, enhancing their effect on entanglement. This is clearly seen
through the qualitative change of the finite-size corrections from equation
(2) to (3).

It must be pointed out that the consistency of our perturbative calculation
requires the term (56) to be subleading to the flat-space area-law term. We
have used an expansion in $1/|\tau|$, but the various orders include powers of
the combination $R^{2}/\tau^{2}$. In terms of the physical quantities
$R_{p}=aR$, $\epsilon_{p}=a\epsilon$ and $H=-1/(a\tau)$, this limits us to
subhorizon entangling radii $\epsilon_{p}\ll R_{p}\ll 1/H$. However, the
overall system size $L_{p}=aL$ can exceed the Hubble radius, as long as it is
not exponentially large to invalidate the expansion.

The presence of a novel term that depends logarithmically on the size of the
_overall system_ is intriguing. This term is not dominant within the region of
validity of the perturbative expansion we employed, and thus it does not
determine the scaling properties of the entanglement entropy. It should be
noted, however, that a presence of a term with a logarithmic dependence on $L$
is likely even in late-time or superhorizon regimes over which we do not have
analytical control at present. The reason is that such a term is incorporated
in the kernel $\Omega^{-3}$ given by (32), which appears in the full matrix
(14).

Even though it is unclear if the term (56) can become dominant in some regime,
its importance cannot be overstated: it represents the leading effect of the
cosmological expansion for subhorizon entangling radii, and also provides a
probe of regions beyond the horizon, which are otherwise inaccessible. A
relevant question concerns the identification of the parameter $L$. For a
finite universe, this quantity is fixed by its total size, which can be so
large that its presence cannot be deduced through curvature effects at
subhorizon distances. But even for a spatially flat infinite universe, there
exists a bound on $L$, set by the duration of inflation. This can be seen by
considering the function $b(\tau)$ that quantifies the effect of squeezing and
determines the wave function (8) of the field modes. The assumption of a
Bunch-Davies vacuum in exact de Sitter space fixes this function to have the
form $b(\tau)=\sqrt{1+1/(\omega^{2}_{0}\tau^{2})}$, so that $b(\tau)\to 1$ and
$b^{\prime}(\tau)\to 0$ for $\tau\to-\infty$. If, however, we assume that
inflation started at some finite time $\tau_{0}$, with no prior squeezing of
the wave function, these conditions must be imposed at $\tau_{0}$. It can be
checked that $b(\tau)$ is not modified considerably for modes with
$\omega_{0}\tau_{0}\gg 1$. Such modes start subhorizon and accumulate all the
squeezing induced by the subsequent expansion. However, for modes with
$\omega_{0}\tau_{0}\ll 1$, the function $b(\tau)$ is suppressed by
$\sim\omega_{0}\tau_{0}$ at late times, relatively to the exact de Sitter
scenario. Such modes are already superhorizon in the beginning of inflation,
so that their total squeezing is reduced. It can then be checked that these
modes do not contribute to the kernel $\Omega^{-3}$ that induces the
logarithmic term. Thus $L$ can be identified with the wavelength of the first
mode that exited the horizon in the beginning of inflation, which may extend
far beyond the present-day horizon.

The dependence of the entanglement entropy on the size of the overall system
implies that the degrees of freedom that lie beyond the cosmological horizon
affect its interior. Classically, this is not permitted. However, the nature
of quantum entanglement allows for this possibility. At the conceptual level
this effect is interesting because of its strong similarity to the EPR
paradox. A possibility to test the effect is provided by analogue gravity
experiments with ultra-cold atoms, which can simulate field evolution in an
expanding background [41, 42, 43, 44]. Similar experiments have measured
entanglement entropy, or related quantities, in a static setting [45, 46, 47].
An experiment that combines elements from both these lines of research could
provide a quantitative check of the dependence of the entropy on the size of
the overall system.

## References

  * Mukhanov and Chibisov [1981] V. F. Mukhanov and G. V. Chibisov, Quantum Fluctuations and a Nonsingular Universe, JETP Lett. 33, 532 (1981).
  * Hawking [1982] S. W. Hawking, The Development of Irregularities in a Single Bubble Inflationary Universe, Phys. Lett. B 115, 295 (1982).
  * Guth and Pi [1982] A. H. Guth and S. Y. Pi, Fluctuations in the New Inflationary Universe, Phys. Rev. Lett. 49, 1110 (1982).
  * Starobinsky [1982] A. A. Starobinsky, Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations, Phys. Lett. B 117, 175 (1982).
  * Bardeen  _et al._ [1983] J. M. Bardeen, P. J. Steinhardt, and M. S. Turner, Spontaneous Creation of Almost Scale - Free Density Perturbations in an Inflationary Universe, Phys. Rev. D 28, 679 (1983).
  * Albrecht  _et al._ [1994] A. Albrecht, P. Ferreira, M. Joyce, and T. Prokopec, Inflation and squeezed quantum states, Phys. Rev. D 50, 4807 (1994), arXiv:astro-ph/9303001 .
  * Polarski and Starobinsky [1996] D. Polarski and A. A. Starobinsky, Semiclassicality and decoherence of cosmological perturbations, Class. Quant. Grav. 13, 377 (1996), arXiv:gr-qc/9504030 .
  * Kiefer  _et al._ [1998a] C. Kiefer, D. Polarski, and A. A. Starobinsky, Quantum to classical transition for fluctuations in the early universe, Int. J. Mod. Phys. D 7, 455 (1998a), arXiv:gr-qc/9802003 .
  * Kiefer  _et al._ [1998b] C. Kiefer, J. Lesgourgues, D. Polarski, and A. A. Starobinsky, The Coherence of primordial fluctuations produced during inflation, Class. Quant. Grav. 15, L67 (1998b), arXiv:gr-qc/9806066 .
  * Allen  _et al._ [2000] B. Allen, E. E. Flanagan, and M. A. Papa, Is the squeezing of relic gravitational waves produced by inflation detectable?, Phys. Rev. D 61, 024024 (2000), arXiv:gr-qc/9906054 .
  * Maldacena [2016] J. Maldacena, A model with cosmological Bell inequalities, Fortsch. Phys. 64, 10 (2016), arXiv:1508.01082 [hep-th] .
  * Grishchuk and Sidorov [1990] L. P. Grishchuk and Y. V. Sidorov, Squeezed quantum states of relic gravitons and primordial density fluctuations, Phys. Rev. D 42, 3413 (1990).
  * Brandenberger  _et al._ [1993] R. H. Brandenberger, T. Prokopec, and V. F. Mukhanov, The Entropy of the gravitational field, Phys. Rev. D 48, 2443 (1993), arXiv:gr-qc/9208009 .
  * Brandenberger  _et al._ [1992] R. H. Brandenberger, V. F. Mukhanov, and T. Prokopec, Entropy of a classical stochastic field and cosmological perturbations, Phys. Rev. Lett. 69, 3606 (1992), arXiv:astro-ph/9206005 .
  * Prokopec [1993] T. Prokopec, Entropy of the squeezed vacuum, Class. Quant. Grav. 10, 2295 (1993).
  * Matacz [1994] A. L. Matacz, The Coherent state representation of quantum fluctuations in the early universe, Phys. Rev. D 49, 788 (1994), arXiv:gr-qc/9212008 .
  * Gasperini and Giovannini [1993] M. Gasperini and M. Giovannini, Entropy production in the cosmological amplification of the vacuum fluctuations, Phys. Lett. B 301, 334 (1993), arXiv:gr-qc/9301010 .
  * Kiefer  _et al._ [2000] C. Kiefer, D. Polarski, and A. A. Starobinsky, Entropy of gravitons produced in the early universe, Phys. Rev. D 62, 043518 (2000), arXiv:gr-qc/9910065 .
  * Srednicki [1993] M. Srednicki, Entropy and area, Phys. Rev. Lett. 71, 666 (1993), arXiv:hep-th/9303048 .
  * Bombelli  _et al._ [1986] L. Bombelli, R. K. Koul, J. Lee, and R. D. Sorkin, A Quantum Source of Entropy for Black Holes, Phys. Rev. D 34, 373 (1986).
  * Lohmayer  _et al._ [2010] R. Lohmayer, H. Neuberger, A. Schwimmer, and S. Theisen, Numerical determination of entanglement entropy for a sphere, Phys. Lett. B 685, 222 (2010), arXiv:0911.4283 [hep-lat] .
  * Casini and Huerta [2009] H. Casini and M. Huerta, Entanglement entropy in free quantum field theory, J. Phys. A 42, 504007 (2009), arXiv:0905.2562 [hep-th] .
  * Callan and Wilczek [1994] C. G. Callan, Jr. and F. Wilczek, On geometric entropy, Phys. Lett. B 333, 55 (1994), arXiv:hep-th/9401072 .
  * Holzhey  _et al._ [1994] C. Holzhey, F. Larsen, and F. Wilczek, Geometric and renormalized entropy in conformal field theory, Nucl. Phys. B 424, 443 (1994), arXiv:hep-th/9403108 .
  * Korepin [2004] V. Korepin, Universality of Entropy Scaling in One Dimensional Gapless Models, Phys. Rev. Lett. 92, 096402 (2004), arXiv:cond-mat/0311056 .
  * Calabrese and Cardy [2004] P. Calabrese and J. L. Cardy, Entanglement entropy and quantum field theory, J. Stat. Mech. 0406, P06002 (2004), arXiv:hep-th/0405152 .
  * Calabrese and Cardy [2009] P. Calabrese and J. Cardy, Entanglement entropy and conformal field theory, J. Phys. A 42, 504005 (2009), arXiv:0905.4013 [cond-mat.stat-mech] .
  * Maldacena and Pimentel [2013] J. Maldacena and G. L. Pimentel, Entanglement entropy in de Sitter space, JHEP 02, 038, arXiv:1210.7244 [hep-th] .
  * Kanno  _et al._ [2014] S. Kanno, J. Murugan, J. P. Shock, and J. Soda, Entanglement entropy of $\alpha$-vacua in de Sitter space, JHEP 07, 072, arXiv:1404.6815 [hep-th] .
  * Iizuka  _et al._ [2016] N. Iizuka, T. Noumi, and N. Ogawa, Entanglement entropy of de Sitter space $\alpha$-vacua, Nucl. Phys. B 910, 23 (2016), arXiv:1404.7487 [hep-th] .
  * Kanno  _et al._ [2017] S. Kanno, M. Sasaki, and T. Tanaka, Vacuum State of the Dirac Field in de Sitter Space and Entanglement Entropy, JHEP 03, 068, arXiv:1612.08954 [hep-th] .
  * Berges  _et al._ [2018] J. Berges, S. Floerchinger, and R. Venugopalan, Dynamics of entanglement in expanding quantum fields, JHEP 04, 145, arXiv:1712.09362 [hep-th] .
  * Boutivas  _et al._ [2023] K. Boutivas, G. Pastras, and N. Tetradis, Entanglement and expansion, JHEP 05, 199, arXiv:2302.14666 [hep-th] .
  * Boutivas  _et al._ [2024] K. Boutivas, D. Katsinis, G. Pastras, and N. Tetradis, Entanglement in cosmology, JCAP 04, 017, arXiv:2310.17208 [gr-qc] .
  * Katsinis  _et al._ [2023] D. Katsinis, G. Pastras, and N. Tetradis, Entanglement of harmonic systems in squeezed states, JHEP 10, 039, arXiv:2304.04241 [hep-th] .
  * Katsinis  _et al._ [2024] D. Katsinis, G. Pastras, and N. Tetradis, Entanglement Entropy of a Scalar Field in a Squeezed State, (2024), arXiv:2403.03136 [hep-th] .
  * [37] K. Boutivas, D. Katsinis, G. Pastras, and N. Tetradis, A Numerical Calculation of Entanglement Entropy in de Sitter Space, to appear .
  * Katsinis and Pastras [2024] D. Katsinis and G. Pastras, Entanglement in $(1+1)$-dimensional Free Scalar Field Theory: Tiptoeing between Continuum and Discrete, (2024), arXiv:2406.11031 [hep-th] .
  * Huerta and van der Velde [2023] M. Huerta and G. van der Velde, Modular Hamiltonian of the scalar in the semi infinite line: dimensional reduction for spherically symmetric regions, JHEP 06, 097, arXiv:2301.00294 [hep-th] .
  * Sorkin [2014] R. D. Sorkin, Expressing entropy globally in terms of (4D) field-correlations, J. Phys. Conf. Ser. 484, 012004 (2014), arXiv:1205.2953 [hep-th] .
  * Eckel  _et al._ [2018] S. Eckel, A. Kumar, T. Jacobson, I. B. Spielman, and G. K. Campbell, A rapidly expanding Bose-Einstein condensate: an expanding universe in the lab, Phys. Rev. X 8, 021021 (2018), arXiv:1710.05800 [cond-mat.quant-gas] .
  * Wittemer  _et al._ [2019] M. Wittemer, F. Hakelberg, P. Kiefer, J.-P. Schröder, C. Fey, R. Schützhold, U. Warring, and T. Schaetz, Phonon Pair Creation by Inflating Quantum Fluctuations in an Ion Trap, Phys. Rev. Lett. 123, 180502 (2019), arXiv:1903.05523 [quant-ph] .
  * Viermann  _et al._ [2022] C. Viermann _et al._ , Quantum field simulator for dynamics in curved spacetime, Nature 611, 260 (2022), arXiv:2202.10399 [cond-mat.quant-gas] .
  * Tolosa-Simeón  _et al._ [2022] M. Tolosa-Simeón _et al._ , Curved and expanding spacetime geometries in Bose-Einstein condensates, Phys. Rev. A 106, 033313 (2022), arXiv:2202.10441 [cond-mat.quant-gas] .
  * Islam  _et al._ [2015] R. Islam, R. Ma, P. M. Preiss, M. E. Tai, A. Lukin, M. Rispoli, and M. Greiner, Measuring entanglement entropy through the interference of quantum many-body twins 10.1038/nature15750 (2015), arXiv:1509.01160 [cond-mat.quant-gas] .
  * 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, aau4963 (2019).
  * Tajik  _et al._ [2023] M. Tajik _et al._ , Verification of the area law of mutual information in a quantum field simulator, Nature Phys. 19, 1022 (2023), arXiv:2206.10563 [cond-mat.quant-gas] .