\section{Introduction}
\label{sec:introduction}

\blfootnote{\textsuperscript{\textdagger}Work done primarily while an intern at ASAPP, Inc.}
Predicting a response given conversational context is a critical task for building open-domain chatbots and dialogue systems.
Recently developed conversational systems typically use either a generative or a retrieval approach for producing responses~\cite{wang2013,ji2014,vinyals2015,serban2015,li2016,xing2016,deb2019diversifying}.
While both of these approaches have demonstrated strong performance in the literature, retrieval methods often enjoy better control over response quality than generative approaches.
In particular, such methods select outputs from a \textit{whitelist} of candidate responses, which can be pre-screened and revised for desired qualities such as sentence fluency and diversity.

Most previous work on retrieval models has concentrated on designing neural architectures to improve response selection.
For instance, several works have improved model performance by encoding multi-turn conversation context instead of single-turn context~\cite{serban2015,zhou2016,wu2017}.
More recent efforts~\cite{zhou2018,zhang2018} have explored using more advanced architectures, such as the Transformer~\cite{vaswani2017}, to better learn the mapping between the context and the candidate responses.

Relatively little effort, however, has been devoted to the practical considerations of using such models in a real-world production setting.
For example, one critical consideration rarely discussed in the literature is the inference speed of the deployed model.
While recent methods introduce rich computation, such as cross-attention \cite{zhou2018}, to improve the modeling between the conversational context and candidate response, the model outputs must be re-computed for every pair of context and response.
As a consequence, these models are not well-suited to a production setting where the size of the response whitelist can easily extend into the thousands.

Another critical concern is the whitelist selection process and the associated retrieval evaluation.
Most prior work have reported Recall@$k$ on a small set of randomly selected responses which include the true response sent by the agent ~\cite{lowe2015,zhou2016,zhou2018,wu2017,zhang2018}.
However, this over-simplified evaluation may not provide a useful indication of performance in production, where the whitelist is not randomly selected, is significantly larger, and may not contain the target response.

In this paper, we explore and evaluate model and whitelist design choices for building retrieval-based conversation systems in production.
We present a dual encoder architecture that is optimized to select among as many as 10,000 responses within a couple tens of milliseconds.
The model makes use of a fast recurrent network implementation~\cite{lei2018} and multi-headed attention~\cite{lin2017} and achieves over a 4.1x inference speedup over traditional encoders such as LSTM~\cite{hochretier1997}.
The independent dual encoding allows pre-computing the embeddings of candidate responses, thereby making the approach highly scalable with the size of the whitelist.
In addition, we compare two approaches for generating the response candidates, and we conduct a comprehensive analysis of our model and whitelists on a large, real-world help desk dataset, using human evaluation and metrics that are more relevant to use in a production setting.

\section{Related Work}
\label{sec:related_work}

This paper extends the line of work on conversational retrieval models for multi-turn response selection \cite{lowe2015,alrfou2016,zhou2016,zhou2018,wu2016,wu2017,yan2016,lu2017,zhang2018,shalyminov2018,deb2019diversifying,yang2019hybrid}. Our model is most similar to \citet{lowe2015}, who construct the context of the conversation by concatenating all previous utterances. They use an RNN to separately encode the context and each candidate response, and they then compute a matching score between the context and response representations to determine the best response for that context.

Other recent work has explored more complex methods of incorporating information from the context of a conversation. \citet{serban2015} and \citet{zhou2016} employ a hierarchical architecture in which they encode the context using RNNs at both the word level and the utterance level. In contrast to these models, which generate a single context encoding, \citet{wu2017} designed a network that matches a response to each utterance in the context individually.

While many of the models cited above implement their RNNs with an LSTM \cite{hochretier1997}, we instead use an SRU \cite{lei2018}. SRU uses light recurrence, which makes it highly parallelizable, and \citet{lei2018} showed that it trains 5-9x faster than cuDNN LSTM. SRU also exhibits a significant speedup in inference time compared to LSTM (by a factor of 4.1x in our experiments), which is particularly relevant in a production setting. Furthermore, \citet{lei2018} showed that SRU matches or exceeds the performance of models using LSTMs or the Transformer architecture \cite{vaswani2017} on a number of NLP tasks, meaning significant speed gains can be achieved without a drop in performance.

Despite the abundance of prior work on retrieval models for dialogue, whitelist selection has received relatively little attention. Since practical use of conversational models has typically not been addressed, most models are evaluated on their ability to select the correct response from a small list of randomly sampled responses \cite{lowe2015}. Another option, from \citet{wu2017}, is to use \textit{Apache Lucene}\footnote{\url{http://lucene.apache.org/}} to select a list of response candidates relevant to each context. However, neither method produces a single whitelist that can be used for every context and reviewed for quality. The closest work to ours is \citet{lu2017}, who build a whitelist using a $k$-means clustering of responses. We extend this work by doing a more comprehensive analysis of different whitelist selection methods, and we further analyze the effect of whitelist size on performance.

\section{Model Architecture}
\label{sec:model_architecture}

\begin{figure*}
\centering
\begin{subfigure}{.45\textwidth}
  \centering
  \includegraphics[width=\linewidth]{model.pdf}
\end{subfigure}
\begin{subfigure}{.45\textwidth}
  \centering
  \includegraphics[width=\linewidth]{loss.pdf}
\end{subfigure}
\caption{(Left) The dual encoder architecture, which takes as input a context $c$ and a response $r$ and computes the score $s(c,r)$. (Right) Computing the model's loss based on the scores between a context $c$ and each response in $R = \{ r^+, r_1^-, \dots, r_k^- \}$, where $r^+$ is the actual agent response and  $r_1^-, \dots, r_k^-$ are randomly sampled agent responses.}
\label{fig:model}
\end{figure*}

Next we describe the architecture of our retrieval model. The two inputs to the model are a context $c$, which is a concatenation of all utterances in the conversation, and a candidate response $r$. In the context, we use special tokens to indicate whether each utterance comes from the customer or the agent. The model outputs a score $s(c,r)$ indicating the relevance of the response to the context. The model architecture is described in detail below and is illustrated in Figure \ref{fig:model}.

\subsection{Dual Encoders}
\label{ssec:dual_encoders}

At the core of our model are two neural encoders $f_c$ and $f_r$ to encode the context and the response, respectively. These encoders have identical architectures but learn separate weights.

Each encoder takes a sequence of tokens $w = \{w_1, w_2, \dots, w_n\}$ as input, which is either a context or a response.
Due to the prevalence of typos in both user and agent utterances in real chats, we use {\tt fastText} \cite{bojanowski2016} as the word embedding method. {\tt fastText} learns both word-level and character-level features and is therefore more robust to misspellings. We pretrained {\tt fastText}\footnote{\url{https://github.com/facebookresearch/fastText}} embeddings on a corpus of 15M utterances from help desk conversations and then fixed the embeddings while training the neural encoders.

Each encoder consists of a recurrent neural network followed by a multi-headed attention layer~\cite{lin2017} to perform pooling.
We use multi-layer, bidirectional SRUs as the recurrent network.
Each layer of the SRU involves the following computation:
\begin{equation}
    \begin{split}
        \textbf{f}_t &= \sigma(\textbf{W}_f \textbf{x}_t + \textbf{v}_f \odot \textbf{c}_{t-1} + \textbf{b}_f) \\
        \textbf{c}_t &= \textbf{f}_t \odot \textbf{c}_{t-1} + (1 - \textbf{f}_t) \odot (\textbf{W} \textbf{x}_t) \\
        \textbf{r}_t &= \sigma(\textbf{W}_r \textbf{x}_t + \textbf{v}_r \odot \textbf{c}_{t-1} + \textbf{b}_r) \\
        \textbf{h}_t &= \textbf{r}_t \odot \textbf{c}_t + (1 - \textbf{r}_t) \odot \textbf{x}_t
    \end{split}
\end{equation}
where $\sigma$ is the sigmoid activation function, $\textbf{W}, \textbf{W}_f, \textbf{W}_r \in \mathbb{R}^{d_h \times d_e}$ are learned parameter matrices, and $\textbf{v}_f, \textbf{v}_r, \textbf{b}_f, \textbf{b}_v \in \mathbb{R}^{d_h}$ are learned parameter vectors.

The multi-headed attention layer compresses the encoded sequence $\textbf{h} = \{\textbf{h}_1, \textbf{h}_2, \dots, \textbf{h}_n\}$ into a single vector.
For each attention head $i$, attention weights are generating with the following computation:
\begin{equation}
    \bm{\alpha}^{(i)} = {\tt softmax}(\sigma(\textbf{h}^T \textbf{W}^{(i)}_a) \textbf{v}^{(i)}_a)
\end{equation}
where $\sigma$ is a non-linear activation function, $\textbf{W}^{(i)}_a \in \mathbb{R}^{d_h \times d_a}$ is a learned parameter matrix, and $\textbf{v}^{(i)}_a \in \mathbb{R}^{d_a}$ is a learned parameter vector.

The encoded sequence representation is then pooled to a single vector for each attention head $i$ by summing the attended representations:
\begin{equation}
    \widetilde{\textbf{h}}^{(i)} = \sum_{j=1}^n \bm{\alpha}^{(i)}_j \textbf{h}_j\ .
\end{equation}

Finally, the pooled encodings are averaged across the $n_h$ attention heads:
\begin{equation}
    \widetilde{\textbf{h}} = \frac{1}{n_h} \sum_{i=1}^{n_h} \widetilde{\textbf{h}}^{(i)}\ .
\end{equation}
The output of the encoder is the vector $f(w) = \widetilde{\textbf{h}}$.

\subsection{Scoring}
\label{ssec:scoring}

To determine the relevance of a response $r$ to a context $c$, our model computes a matching score between the context encoding $f_c(c)$ and the response encoding $f_r(r)$. This score is simply the dot product of the encodings:
\begin{equation}
    s(c, r) = f_c(c) \cdot f_r(r)\ .
\end{equation}

\subsection{Training}
\label{ssec:training}

We optimize the model to maximize the score between the context $c$ and the response $r^+$ actually sent by the agent while minimizing the score between the context and each of $k$ random (``negative'') responses $r_1^-, \dots, r_k^-$. This is accomplished by training the model to minimize the cross-entropy loss:
\begin{equation}
    \mathcal{L} = -s(c, r^+) + \log \sum_{r \in R} e^{s(c, r)}
    \label{eq:loss}
\end{equation}
where $R = \{r^+, r_1^-, \dots, r_k^-\}$.

Although negative responses could be sampled separately for each context-response pair, we instead use a method inspired by \citet{logeswaran2018} and share a set of negative responses across all examples in a batch. Specifically, for each batch, we sample $k$ responses from the set of all agent responses (weighted according to response frequency), and we use those $k$ responses as the negative responses for every context-response pair in the batch. This has the benefit of reducing the number of responses that need to be encoded in each batch of size $b$ from $\mathcal{O}(bk)$ to $\mathcal{O}(b + k)$, thereby significantly accelerating training.

\begin{table*}[!h]
\centering
\begin{tabular}{|c|c|}
    \hline
    Conversations & 594,555 \\
    \hline
    Utterances & 15,217,773 \\
    \hline
    Customer utterances & 6,943,940 \\
    \hline
    Agent utterances & 8,273,833 \\
    \hline
    Mean conversation length (\# utterances) & 25.60 \\
    \hline
    Mean utterance length (\# tokens) & 12.70 \\
    \hline
    Mean customer utterance length (\# tokens) & 7.53 \\
    \hline
    Mean agent utterance length (\# tokens) & 17.15 \\
    \hline
\end{tabular}
\caption{Summary statistics for the propriety help desk dataset.}
\label{tab:data_stats}
\end{table*}

\subsection{Whitelist Generation}
\label{ssec:whitelist_generation}

After training, we experimented with two methods of creating the whitelist from which our model selects responses at inference time. For each method, we created both a 1,000 response whitelist and a 10,000 response whitelist. Having a whitelist with any more than 10,000 responses would likely make a manual review infeasible.

\paragraph{Frequency-Based Method.}

Responses that are sent frequently are more likely to be relevant in multiple conversations and are less likely to contain errors. Therefore, one method of building a high-quality whitelist is simply to collect messages that are sent often. We created frequency-based whitelists by selecting the 1,000 or 10,000 most common agent responses, after accounting for minor variations in capitalization, punctuation, and whitespace.

\paragraph{Clustering-Based Method.}

Although selecting responses based on frequency may help guarantee quality, manual examination of the frequency whitelists showed that they contained many redundant responses. Therefore, we experimented with a clustering-based whitelist selection method in the hope of reducing redundancy and increasing response diversity. Specifically, we encoded all agent responses using our response encoder $f_r$ and then used $k$-means clustering with $k=1,000$ or $k=10,000$ to cluster the responses. We then selected the most common response from each cluster to create the whitelists.

\section{Experiments and Results}
\label{sec:experiments_and_results}

We evaluated our model and whitelists on a large, proprietary help desk chat dataset using several offline metrics and a human evaluation. We particularly emphasize metrics relevant to production, such as inference speed and Recall@$k$ from a large candidate set. The human evaluation illustrates how our model and whitelists compare to each other and to the responses sent by a real human agent.

\subsection{Data}
\label{ssec:data}

\begin{table}
\centering
\begin{tabular}{|p{7.5cm}|}
    \hline
    \multicolumn{1}{|c|}{\textbf{Example Conversation}} \\
    \Xhline{3\arrayrulewidth}
    \textbf{Customer:} I would like to pay my bill can you help me \\
    \textbf{Agent:} I can definitely help you to pay your bill. Are we going to work with the account logged in now? \\
    \textbf{Customer:} Yes it still says there is no money on my account \\
    \textbf{Agent:} I understand that. I have reviewed your account and its shows here that the payment has been posted and you're all good until next month service. \\
    \textbf{Customer:} Oh ok thank you for all your help \\
    \textbf{Agent:} You're welcome. Anything for a valued customer like you! \\
    \hline
\end{tabular}
\caption{A sample conversation from the propriety help desk chat dataset. The sample has been lightly edited to remove proprietary information.}
\label{tab:chat_example}
\end{table}

The help desk chat dataset used in our experiments consists of 15M utterances from 595K conversations. We randomly split the conversations into train, validation, and test sets with 80\%, 10\%, and 10\% of the conversations, respectively. Since each conversation includes several agent responses, each of which produces a context-response example, our dataset consists of 6.6M training examples, 828K validation examples, and 828K test examples. Additional dataset statistics are provided in Table \ref{tab:data_stats}. An example chat conversation can be seen in Table \ref{tab:chat_example}.

\subsection{Model Details}
\label{ssec:model_details}

We implemented the dual encoder model using PyTorch~\cite{paszke2017}.
We use pre-trained {\tt fastText} embeddings of dimension $d_e=300$, a 4-layer bidirectional SRU\footnote{SRU code available at \url{https://github.com/taolei87/sru/tree/master/sru}} with hidden size $d_h=300$, and multi-headed attention with 16 heads and a hidden size of $d_a=64$. The batch size was 200 and we used $k=200$ negative responses for each positive response. To ensure quick encoding even for long inputs, contexts were restricted to the 500 most recent tokens and responses were restricted to the 100 most recent tokens\footnote{A context with 500 tokens contains 39 utterances on average, which is typically more than enough to understand the topic of conversation. Almost all responses are shorter than 100 tokens.}. The model was optimized using Adam~\cite{kingma2014} with the Noam learning rate schedule from ~\citet{vaswani2017}. The model was trained for 30 epochs, with each epoch limited to 10,000 training batches (2M training examples). Training took about 32 hours on a single Tesla V100 GPU.

\subsection{Results and Analysis}
\label{ssec:performance_and_analysis}

\paragraph{AUC and AUC@$p$.}

To determine the model's ability to use context to distinguish between true responses and negative responses, we use the metrics AUC and AUC@$p$. AUC is the area under the receiver operating characteristic curve when using the score $s(c,r)$ to determine whether each response is the true response or a negative response. AUC@$p$ is the area under the portion of the ROC curve where the false positive rate is $\leq p$, renormalized so that the maximum AUC@$p$ is 1.

The performance of our model according to these AUC metrics can be seen in Table~\ref{tab:auc}. The high AUC indicates that our model can easily distinguish between the true response and negative responses. Furthermore, the AUC@$p$ numbers show that the model has a relatively high true positive rate even under the difficult requirement of a low false positive rate.

\begin{table}
\centering
\begin{tabular}{|c|c|c|}
    \hline
    \textbf{Metric} & \textbf{Validation} & \textbf{Test} \\
    \Xhline{3\arrayrulewidth}
    AUC & 0.991 & 0.977 \\
    AUC@0.1 & 0.925 & 0.885 \\
    AUC@0.05 & 0.871 & 0.816 \\
    AUC@0.01 & 0.677 & 0.630 \\
    \hline
\end{tabular}
\caption{AUC and AUC@$p$ of our model on the propriety help desk dataset.}
\label{tab:auc}
\end{table}

\begin{table}
\centering
\begin{tabular}{|c|c|c|c|c|}
    \hline
    \textbf{Candidates} & \textbf{R@1} & \textbf{R@3} & \textbf{R@5} & \textbf{R@10} \\
    \Xhline{3\arrayrulewidth}
    10 & 0.892 & 0.979 & 0.987 & 1 \\
    100 & 0.686 & 0.842 & 0.894 & 0.948 \\
    1,000 & 0.449 & 0.611 & 0.677 & 0.760 \\
    10,000 & 0.234 & 0.360 & 0.421 & 0.505 \\
    \hline
\end{tabular}
\caption{Recall@$k$ from $n$ response candidates for different values of $n$ using random whitelists. Each random whitelist includes the correct response along with $n-1$ randomly selected responses.}
\label{tab:random_recall}
\end{table}

\paragraph{Recall and Whitelist Size.}

In order to determine our model's ability to select the correct response from a whitelist, we use recall at $k$ from $n$ (R\textsubscript{$n$}@$k$), which is the proportion of times that the true response is ranked as one of the top $k$ responses in a whitelist containing $n$ candidate responses.

Table \ref{tab:random_recall} shows R\textsubscript{$n$}@$k$ on the test set for different values of $n$ and $k$ when using a random whitelist, meaning a whitelist which contains the true response and $n-1$ randomly sampled responses\footnote{To be precise, we sampled responses without replacement weighted according to the frequency with which the response was sent by agents.}. As discussed in the introduction, most prior work evaluate their models using a random whitelist with $n=10$ candidates. However, a production whitelist needs to contain hundreds or thousands of response candidates in order to provide relevant responses in a variety of contexts. Therefore, a more meaningful metric for production purposes is R\textsubscript{$n$}@$k$ for $n \geq 100$. Table \ref{tab:random_recall} shows that recall drops significantly as $n$ grows, meaning that the R\textsubscript{10}@$k$ evaluation performed by prior work may significantly overstate model performance in a production setting.

\begin{table}
\centering
\resizebox{0.49\textwidth}{!}{
    \begin{tabular}{|c|c|c|c|c|c|}
        \hline
        \textbf{Whitelist} & \textbf{R@1} & \textbf{R@3} & \textbf{R@5} & \textbf{R@10} & \textbf{BLEU} \\
        \Xhline{3\arrayrulewidth}

        Random 10K+ & 0.252 & 0.400 & 0.472 & 0.560 &  37.71 \\
        Frequency 10K+ & 0.257 & 0.389 & 0.455 & 0.544 & 41.34 \\
        Clustering 10K+ & 0.230 & 0.376 & 0.447 & 0.541 &  37.59\\
        \hline

        Random 1K+ & 0.496 & 0.663 & 0.728 & 0.805 & 59.28\\
        Frequency 1K+ & 0.513 & 0.666 & 0.726 & 0.794 & 67.05\\
        Clustering 1K+ & 0.481 & 0.667 & 0.745 & 0.835 & 61.88 \\
        \hline
        \hline
        Frequency 10K & 0.136 & 0.261 & 0.327 & 0.420 & 30.46\\
        Clustering 10K & 0.164 & 0.292 & 0.360 & 0.457 &  31.47\\
        \hline
        Frequency 1K & 0.273 & 0.465 & 0.550 & 0.658 &  47.13 \\
        Clustering 1K & 0.331 & 0.542 & 0.650 & 0.782 & 49.26\\
        \hline
    \end{tabular}
}
\caption{Recall@$k$ for random, frequency, and clustering whitelists of different sizes. The ``+'' indicates that the true response is added to the whitelist.}
\label{tab:random_vs_frequency_vs_clustering}
\end{table}

\begin{table}
\centering
\begin{tabular}{|c|c|c|c|c|}
    \hline
    \textbf{Whitelist} & \textbf{R@1} & \textbf{Coverage} \\
    \Xhline{3\arrayrulewidth}

    Frequency 10K & 0.136 & 45.04\% \\
    Clustering 10K & 0.164 & 38.38\% \\
    \hline

    Frequency 1K & 0.273 & 33.38\% \\
    Clustering 1K & 0.331 & 23.28\% \\
    \hline
\end{tabular}
\caption{Recall@1 versus coverage for frequency and clustering whitelists.}
\label{tab:frequency_vs_clustering}
\end{table}

\paragraph{Comparison Between Whitelists.}
An interesting question we would like to address is whether a random whitelist serves as a good proxy for whitelists generated using other methods.
To this end, we also evaluate recall on the frequency and clustering whitelists from Section \ref{ssec:whitelist_generation}.

First, we compute recall when the true response is added to the whitelist, as in the case of the random whitelists described above. Second, we compute recall only on the subset of examples for which the true response is already contained in the whitelist. The latter recall measure is more relevant to a production setting since the true response is not known at inference time and therefore cannot be artificially added to the whitelist.

The results in Table \ref{tab:random_vs_frequency_vs_clustering} show that the three types of whitelists perform comparably to each other when the true response is added.
However, in the more realistic second case, when recall is only computed on examples with a response already in the whitelist, performance on the frequency and clustering whitelists drops significantly.

Additionally, we compute the BLEU ~\cite{papineni2002bleu, ward2002corpus} scores between the true responses and the best suggested responses. The BLEU score allows us to measure the semantic similarity when the true and suggested responses are not exactly matched. The BLEU scores computed with the frequency and clustering whitelists are slightly higher than those computed with random whitelists.

\paragraph{Recall versus Coverage.}

Although recall is a good measure of performance, recall alone is not a sufficient criterion for whitelist selection. The recall results in Table \ref{tab:random_vs_frequency_vs_clustering} seem to indicate that the clustering-based whitelists are strictly superior to the frequency-based whitelists in the realistic case when we only consider responses that are already contained in the whitelist, but this analysis fails to account for the frequency with which this is the case. For instance, a whitelist may have very high recall but may only include responses that were sent infrequently by agents, meaning the whitelist will perform well for a handful of conversations but will be irrelevant in most other cases.

To quantify this effect, we introduce the notion of \textit{coverage}, which is the percent of all context-response pairs where the agent response appears in the whitelist, after accounting for minor deviations in capitalization, punctuation, and whitespace. A whitelist that contains responses that are sent more frequently by agents will therefore have a higher coverage.

Table \ref{tab:frequency_vs_clustering} shows R@1 and coverage for the frequency and clustering whitelists. While the clustering whitelists have higher recall, the frequency whitelists have higher coverage. This is to be expected since the frequency whitelists were specifically chosen to maximize the frequency of the included responses. Since both recall and coverage are necessary to provide good responses for a wide range of conversations, these results indicate the importance of considering the trade-off between recall and coverage inherent in a given whitelist selection method.

It may be interesting in future work to further investigate these trade-offs in order to identify a whitelist selection method that can simultaneously optimize recall and coverage.

\begin{table}
\centering
\resizebox{0.49\textwidth}{!}{
    \begin{tabular}{|c|c|c|c?{0.5mm}c|}
        \hline
        \textbf{Whitelist} & \textbf{Great} & \textbf{Good} & \textbf{Bad} & \textbf{Accept} \\
        \Xhline{3\arrayrulewidth}
        Freq. 1K & 54\% & 26\% & 20\% & 80\% \\
        Cluster. 1K & 55\% & 21\% & 23\% & 77\% \\
        Freq. 10K & 56\% & 24\% & 21\% & 80\% \\
        Cluster. 10K & 57\% & 23\% & 20\% & 80\% \\
        \hline
        Real response & 60\% & 24\% & 16\% & 84\% \\
        \hline
    \end{tabular}
}
\caption{Results of the human evaluation of the responses produced by our model. A response is acceptable if it is either good or great. Note: Numbers may not add up to 100\% due to rounding.}
\label{tab:human}
\end{table}

\paragraph{Human Evaluation.}
While offline metrics are indicative of model performance, the best measure of performance is a human evaluation of the model's predictions. Therefore, we performed a small-scale human evaluation of our model and whitelists. We selected 322 contexts from the test set and used our model to generate responses from the Frequency 10K, Frequency 1K, Clustering 10K, and Clustering 1K whitelists. Three human annotators were shown each context followed by five responses: one from each of the four whitelists and the true response sent by the agent. The annotators were blinded to the source of each response. The annotators were asked to rate each response according to the following categories:

\vspace{3mm}

\noindent
\textbf{Bad:} The response is not relevant to the context.

\noindent
\textbf{Good:} The response is relevant to the context but is vague or generic.

\noindent
\textbf{Great:} The response is relevant to the context and directly addresses the issue at hand.

\vspace{3mm}

\noindent
For example, three such responses for the context ``My phone is broken'' would be:

\vspace{3mm}

\noindent
\textbf{Bad response:} Goodbye!

\noindent
\textbf{Good response:} I'm sorry to hear that.

\noindent
\textbf{Great response:} I'm sorry to hear that your phone is broken.

\vspace{3mm}

The results of the human evaluation are in Table \ref{tab:human}. Our proposed system works well, selecting acceptable (i.e. good or great) responses about 80\% of the time and selecting great responses more than 50\% of the time.

Interestingly, the size and type of whitelist seem to have little effect on performance, indicating that all the whitelists contain responses appropriate to a variety of conversational contexts. Since the frequency whitelists are simpler to generate than the clustering whitelists and since the 1K whitelists contain fewer responses to manually review than the 10K whitelists, the Frequency 1K whitelist is the preferred whitelist for our production system.

\begin{table}
\centering
\begin{tabular}{|c|c|c|c|}
    \hline
    \textbf{Encoder} & \textbf{Layer} & \textbf{Params} & \textbf{Time} \\
    \Xhline{3\arrayrulewidth}
    SRU & 2 & 3.7M & 14.7 \\
    SRU & 4 & 8.0M & 21.9 \\
    LSTM & 2 & 7.3M & 90.9 \\
    LSTM & 4 & 15.9M & 174.8  \\

    \hline
    $ $+rank response& - & - & 0.9 \\
    \hline
\end{tabular}

\caption{Inference time (milliseconds) of our model to encode a context using an SRU or an LSTM encoder on a single CPU core. The last row shows the extra time needed to compare the response encoding to 10,000 cached candidate response encodings in order to find the best response.}
\label{tab:inference_speed}
\end{table}
\paragraph{Inference Speed.}
A major constraint in a production system is the speed with which the system can respond to users.
To demonstrate the benefit of using an SRU encoder instead of an LSTM encoder in production, we compared the speed with which they encode a random conversation context at inference time, averaged over 1,000 samples.
We used a single core of Intel Core i9 2.9 GHz CPU.
As seen in Table \ref{tab:inference_speed}, an SRU encoder is over 4x faster than an LSTM encoder with a similar number of parameters, making it more suitable for production use.

Table~\ref{tab:inference_speed} also highlights the scalability of using a dual encoder architecture.
Since the embeddings of the candidate responses are independent from the conversation context, the embeddings of the whitelist responses can be pre-computed and stored as a matrix.
Retrieving the best candidate once the context is encoded takes a negligible amount of time compared to the time to encode the context.

\begin{table*}[!t]
\centering
\begin{tabular}{|c|c|c|c|c|}
    \hline
    \textbf{Model} & \textbf{Parameters} & \textbf{Validation AUC@0.05} & \textbf{Test AUC@0.05} \\
    \Xhline{2\arrayrulewidth}
    Base & 8.0M & \textbf{0.871} & 0.816 \\
    \hline
    4L SRU $\rightarrow$ 2L LSTM & 7.3M & 0.864 & \textbf{0.829} \\
    \hline
    4L SRU $\rightarrow$ 2L SRU & 7.8M & 0.856 & \textbf{0.829} \\
    \hline
    Flat $\rightarrow$ hierarchical & 12.4M & 0.825 & 0.559 \\
    \hline
    Cross entropy $\rightarrow$ hinge loss & 8.0M & 0.765 & 0.693 \\
    \hline
    6.6M $\rightarrow$ 1M examples & 8.0M & 0.835 & 0.694 \\
    \hline
    6.6M $\rightarrow$ 100K examples & 8.0M & 0.565 & 0.417 \\
    \hline
    200 $\rightarrow$ 100 negatives & 8.0M & 0.864 & 0.647 \\
    \hline
    200 $\rightarrow$ 10 negatives & 8.0M & 0.720 & 0.412 \\
    \hline
\end{tabular}
\caption{An ablation study showing the effect of different model architectures and training regimes on performance on the proprietary help desk dataset.}
\label{tab:ablation}
\end{table*}

\paragraph{Ablation analysis.}
Finally, we performed an ablation analysis to identify the effect of different aspects of the model architecture and training regime. The results are shown in Table \ref{tab:ablation}, and details of the model variants are available in the Appendix.

As Table \ref{tab:ablation} shows, the training set size and the number of negative responses for each positive response are the most important factors in model performance. The model performs significantly worse when trained with hinge loss instead of cross-entropy loss, indicating the importance of the loss function. We also experimented with a hierarchical encoder, where two different recurrent neural networks are used to encode contexts, one at the word level and one at the utterance level. We observed no advantage to using a hierachical encoder, despite its complexity and popularity for encoding conversations \cite{serban2015,zhou2016}. Finally, we see that a 2 layer LSTM performs similarly to either a 4 layer or a 2 layer SRU with a comparable number of parameters. Since the SRU is more than 4x faster at inference time with the same level of performance, it is the preferred encoder architecture.

\section{Conclusion}
\label{sec:conclusion}

In this paper, we present a fast dual encoder neural model for retrieval-based human-computer conversations.
We address technical considerations specific to the production setting, and we evaluate our model and two whitelist generation methods on a large help desk chat dataset. We observe that traditional offline evaluation metrics significantly overestimate model performance, indicating the importance of using evaluation metrics more relevant to a production setting. Furthermore, we find that our proposed model performs well, both on offline metrics and on a human evaluation. Due to its strong performance and its speed at inference time, we conclude that our proposed model is suitable for use in a production conversational system.

One important direction for future work is a deeper analysis of the whitelist selection process. Although our analysis found similar performance across whitelists according to a human evaluation, our offline metrics indicate underlying trade-offs between different characteristics of the whitelists such as recall and coverage. A better understanding the implications of these trade-offs may lead to improved whitelist generation methods, thereby further improving the performance of retrieval-based models.

\section*{Acknowledgments}
\label{sec:acknowledgments}

We would like to thank Howard Chen for the invaluable conversations we had with him during the development of our model. We would also like to thank Anna Folinsky and the ASAPP annotation team for their help performing the human evaluation, and Hugh Perkins for his support on the experimental environment setup. Thank you as well to Ethan Elenberg, Kevin Yang, and Adam Yala for reviewing early drafts of this paper and providing valuable feedback. Finally, thank you to the anonymous reviewers for their constructive feedback.

\section{Appendix}
\label{sec:appendix}

\subsection{Ablation Study}
\label{ssec:ablation_study}

Table \ref{tab:ablation} shows the results of an ablation study we performed to identify the most important components of our model architecture and training regime. Each variant is described below.

\vspace{3mm}

\noindent
\textbf{Base.} This is the model architecture described in Section \ref{sec:model_architecture}.

\vspace{3mm}

\noindent
\textbf{4L SRU $\rightarrow$ 2L LSTM.} We replace the 4 layer SRU encoder with a 2 layer LSTM encoder, which has a comparable number of parameters when using the same hidden sizes.

\vspace{3mm}

\noindent
\textbf{4L SRU $\rightarrow$ 2L SRU.} We use an SRU with 2 layers instead of 4 layers. In order to match parameters, we use a hidden size of $d_h=475$ instead of $d_h=300$ in the model with 2 layers.

\vspace{3mm}

\noindent
\textbf{Flat $\rightarrow$ hierarchical.} We replace the SRU encoder with two SRU encoders, one which operates at the word level and one which operates at the utterance level, following the architectures of \citet{serban2015,wu2017}.

\vspace{3mm}

\noindent
\textbf{Cross entropy $\rightarrow$ hinge loss.} Instead of using the cross-entropy loss defined in Equation \ref{eq:loss}, we use the hinge loss, which is defined as:
\begin{equation}
    \mathcal{L} = \sum_{i=1}^k |s(c,r^+) - s(c,r_i^-) + m|\
\end{equation}
where the margin $m=0.25$ encourages separation between the score of the correct response and the score of each negative response.

\vspace{3mm}

\noindent
\textbf{6.6M $\rightarrow$ 1M examples.} We train on a dataset with 1 million examples instead of the full 6.6 million training examples.

\vspace{3mm}

\noindent
\textbf{6.6M $\rightarrow$ 100K examples.} We trained on a dataset with 100 thousand examples instead of the full 6.6 million training examples.

\vspace{3mm}

\noindent
\textbf{200 $\rightarrow$ 100 negatives.} During training, we sample 100 negatives for each context-response pair instead of 200 negatives.

\vspace{3mm}

\noindent
\textbf{200 $\rightarrow$ 10 negatives.} During training, we sample 10 negatives for each context-response pair instead of 200 negatives.


\section{Introduction}

The fact that one can incorporate thermodynamical attributes to a black hole was
first introduced in the seminal work of Bekenstein \cite{Bekenstein:1972tm,
Bekenstein:1973ur}. In these articles and others \cite{Bekenstein:1974ax,
Kallosh:1992wa, hawking1975, Bekenstein:1975tw, Hawking:1976de, Gibbons:1976pt,
Davies:1978mf, Israel:2003kv, Ashtekar:1997yu, Mukohyama:1997aq, Bardeen:1973gs}
the authors demonstrated that intrinsic entropy $\mathcal{S}_{BH}$ of a black
hole should be proportional to the area $\mathcal{A}_{h}$ of its event horizon
$\mathcal{S}_{BH} = \frac{1}{4}M_{Pl}^2 \mathcal{A}_{h}$, where $M_{Pl}$ is the
Planck mass. Then the natural question appeared is how to connect the concept of
quantum states to this entropy of event horizon \cite{Frolov:1993ym,
Barvinsky:1994jca, Kiefer:1998mi, Strominger:1996sh} as horizon is not different
than any other classical surface with no special local dynamics. To answer this
question and to provide a more general realization of the entropy associated to
a black hole the authors in \cite{Bombelli:1986rw, Srednicki:1993im} presented
the idea in terms of entanglement entropy. Here it is shown that entanglement
entropy of a free scalar field in a certain spatial region is proportional to
its area. In these articles the reduced density matrix, essential for estimating
the entanglement entropy, is obtained by tracing over the spatial degrees of
freedom of the ground state density matrix residing inside the considered
region.

In the regular formulation of entanglement entropy estimation
\cite{Srednicki:1993im} firstly the scalar field is partially Fourier
transformed with respect to the angular coordinates. The resulting Fourier field
Hamiltonian is still dependent on the radial coordinate and it is discretized by
assuming a lattice of finite size and inter-atomic spacing. This discretization
transforms the Fourier Hamiltonian to be a collection of coupled harmonic
oscillators. The ground state wave-function for these coupled harmonic
oscillators then provide the corresponding ground state density matrix for the
field. Subsequently using the Gaussian nature of this ground state wave-function
the reduced density matrix and its eigenvalues are obtained which would produce
the entanglement entropy. However this Gaussian nature is a feature specific to
the Schrodinger quantization. In polymer quantization \cite{Ashtekar:2002sn,
Halvorson-2004-35, Hossain:2010eb}, the quantization method used in loop quantum
gravity \cite{Ashtekar:2004eh, Rovelli2004quantum, Thiemann2007modern}, the
ground state wave-functions are expressed in terms of Mathieu functions. Using
these polymer wave-functions it is still unknown how to obtain the analytic form
of reduced density matrix.

In this article we consider a perturbative approach to circumvent these
difficulties and obtain the entanglement entropy for free scalar field using
Fock and polymer quantization. We treat the interaction between coupled harmonic
oscillators in perturbative manner to get the related ground state and
eigen-values of the reduced density matrix. Firstly we use this procedure to
evaluate the entanglement entropy for two coupled harmonic oscillators in
Schrodinger and polymer quantization. Then by considering free scalar field we
obtain the area law in Fock quantization. Furthermore we apply polymer
quantization in this formulation and verify that the field theoretic
entanglement entropy obeys a similar area law.

In section \ref{EE_usual} we briefly review the procedures to derive the
entanglement entropy in usual formulation. In this section the detailed
description of the considered system is given. Following in section
\ref{EE_perturbed} we recall the perturbative formulation and construct the
framework to estimate the entanglement entropy utilizing this technique. In the
subsequent sections we use this formulation to obtain entanglement entropy for
two coupled harmonic oscillators in Schrodinger and polymer quantization.
Following parts include the realization of the area law of entanglement entropy
in Fock and polymer quantization utilizing perturbative formulation. We argue on
the implications of the obtained results and conclude with the discussion.

\section{Entanglement entropy and the area law}\label{EE_usual}

In the standard derivations of entanglement entropy \cite{Bombelli:1986rw,
Srednicki:1993im, Das:2008sy, Das:2005ah, Ahmadi:2005mh, Das:2007mj, Das:2007sj,
Muller:1995mz, Das:2006vx, Jonker2016EntanglementEO} one considers a system of
coupled harmonic oscillators as a basis. In particular the eigen-values of the
reduced density matrix for two coupled oscillators give the entanglement entropy
corresponding to a single oscillator. These eigenvalues are used for a set of
coupled harmonic oscillators, which are obtained from the discretized
Hamiltonian of a free scalar field, to get the area law of entanglement entropy.
In this section we briefly review the key aspects of these procedures and the
considered systems, which will also be useful to construct the perturbative
formulation.

\subsection{Entanglement entropy for two coupled harmonic oscillators}

In order to understand entropy from entanglement at first a system of two
coupled harmonic oscillators \cite{Srednicki:1993im, Das:2008sy} is considered.
The two unit mass oscillators are denoted by their position and momentum
$(x_1,p_1)$ and $(x_2,p_2)$. The total system can be described by the
Hamiltonian
\begin{eqnarray}
  H &=& \frac{1}{2}\left[p_1^2 + p_2^2+\omega_0^2 (x_1^2+
x_2^2)+k_1^2(x_1-x_2)^2 \right] \label{eq:hamiltonian-full1} \nonumber \\
&=& \frac{1}{2}\left[p_+^2 + \omega_+^2 x_+^2\right] + \frac{1}{2}\left[p_-^2 +
\omega_-^2 x_-^2\right] ~,
\end{eqnarray}
where the normal coordinates $x_\pm=(x_1\pm x_2)/\sqrt{2}$, $p_\pm=(p_1\pm
p_2)/\sqrt{2}$ and normal frequencies $\omega_+=\omega_0$,
$\omega_-=(\omega_0^2+2k_1^2)^{1/2}$ are defined to make the Hamiltonian
decoupled. In decoupled form the ground state wave-function becomes simplified
and can be expressed in terms of the normal coordinates as
\begin{equation}
 \psi_0(x_1,x_2) = \left( \frac{\omega_+ \omega_-}{\pi^2}
\right)^{\frac{1}{4}}
\exp\left[ -\frac{\omega_+ x_+^2 + \omega_- x_-^2}{2}
\right]~.\label{eq:GS-wavefn-normal-freq}
\end{equation}
From the expression (\ref{eq:GS-wavefn-normal-freq}) one can find out the ground
state density matrix to be $\rho(x_1,x_2;x_1',x_2') =
\psi_0(x_1,x_2)~\psi_0^{*}(x_1',x_2')$. To discuss about the entanglement
entropy corresponding to a single oscillator  one needs to find its associated
reduced density matrix. The reduced density matrix is obtained by tracing out
the density matrix with respect to the position degree of freedom of a single
oscillator, expressed as
\begin{equation}
\rho^{r}(x_2,x_2') = \int_{-\infty}^{\infty}dx_1 ~
\psi_0(x_1,x_2)~\psi_0^{*}(x_1,x_2')~.\label{eq:reduced-general}
\end{equation}
The reduced density matrix describes whether the system is in mixed or pure
state and the corresponding entanglement entropy is defined as $S_{E} =
-\mathrm{Tr}\hspace{1pt}[\rho^{r} \ln\rho^{r}]$. In a suitable basis one can evaluate the
entanglement entropy by obtaining the eigenvalues of the reduced density matrix.
In particular for two coupled harmonic oscillators the resulting reduced density
matrix from equation (\ref{eq:reduced-general}) has eigenvalues
\begin{equation}
\lambda_n = (1-\xi^2)\;\xi^{2n}~~\textup{where}~~\xi=\frac{\sqrt{\omega_-}
-\sqrt{\omega_+} } { \sqrt{\omega_-}
+\sqrt{\omega_+}}~.\label{eq:eigenvalue-gen}
\end{equation}
Then the corresponding entanglement entropy \cite{Srednicki:1993im,
Das:2008sy, Chandran:2018wwc} becomes
\begin{eqnarray}
 S_{E}(\xi) = -\sum_n\lambda_n\ln{\lambda_n}=
-\ln{(1-\xi^2)}-\frac{\xi^2}{1-\xi^2}\ln{\xi^2}.\label{eq:ent-entropy-2HO-gen}
\end{eqnarray}

\subsection{Entanglement entropy for N-coupled harmonic oscillators}

Now it is important to understand the entanglement entropy corresponding to
$N-$coupled harmonic oscillators to get the area law of entanglement entropy for
free scalar field. The general Hamiltonian for $N$ coupled harmonic oscillator
is
\begin{equation}
 H = \frac{1}{2} \sum_{j=1}^{N} p_{j}^2 + \frac{1}{2} \sum_{j,k=1}^{N}
x_{j}K_{jk}x_{k}~,\label{eq:general-Ncoupled-HO-hamiltonian}
\end{equation}
where the matrix $K$ describes the potential and interaction. The diagonal
elements of $K$ give the frequency square of individual oscillator and symmetric
off diagonal elements provide the interaction between two adjacent oscillators.
With the help of a suitably chosen orthogonal matrix $U$ this interaction matrix
is diagonalized to $K_D$ as $K=U^T K_D U$. The ground state wave-function of
this $N-$coupled harmonic oscillator (\ref{eq:general-Ncoupled-HO-hamiltonian})
can be expressed as
\begin{equation}
 \psi_{0}(x_1,...,x_N) = \left(\frac{Det.\Omega}{\pi^N}\right)^{\frac{1}{4}}
\exp{\left[ -\frac{x.\Omega.x}{2} \right]}~,\label{eq:general-Ncoupled-HO-GS}
\end{equation}
where $\Omega=U^T K^{1/2}_D U$. From this wave-function one can obtain the
reduced density matrix when first $\mathbbm{n}$ of the total $N$ oscillators are
traced out \cite{Srednicki:1993im}. The reduced density matrix is further
evaluated using a general form of the matrix $\Omega$,
\begin{equation}
 \Omega = \left( {\begin{array}{ccccc}
   A & B\\
   B^T & C\\
  \end{array} } \right)~,
\end{equation}
where $A$ is a $\mathbbm{n}\times \mathbbm{n}$ matrix corresponding to the first
$\mathbbm{n}$ oscillators, $C$ is a $(N-\mathbbm{n})\times(N-\mathbbm{n})$
matrix and $B$ is a $\mathbbm{n}\times(N-\mathbbm{n})$ matrix. In terms of few
newly defined quantities $\beta=(1/2)B^T A^{-1}B$ and $\gamma= C-\beta$ the
reduced density matrix becomes
\begin{equation}
 \rho_{out}(x,x')\sim \exp{\left[ -(x.\gamma.x+x'.\gamma.x')/2+x.\beta.x'
\right]}~,\label{eq:Rrho_Nc}
\end{equation}
where $x$ and $x'$ consist of the $(N-\mathbbm{n})$ oscillators after the
integration over the first $\mathbbm{n}$ degrees of freedom.
$x=V^T\gamma^{-1/2}_{D}y$ is defined, where $\gamma=V^T\gamma_{D}V$ such that
$\gamma_{D}$ is diagonal and $V$ is orthogonal. Then one shall get
$\rho_{out}(x,x')\sim \exp{\left[ -(y.y+y'.y')/2+y.\beta'.y' \right]}$, where
$\beta'=\gamma^{-1/2}_{D} V\beta V^T \gamma^{-1/2}_{D}$. Now moving to the basis
$z=W^{T}y$, such that $\beta'$ is diagonalized as  $\beta'_{D}=W^T\beta'W$, one
gets
\begin{equation}\label{eq:rho_nreduced}
 \rho_{out}(z,z')\sim \prod_{j=\mathbbm{n}+1}^{N}\exp{\left[
-(z_{j}^2+{z'}_{j}^{2})/2+\beta'_{j} z_{j}z'_{j} \right]}~,
\end{equation}
where $\beta'_{j}$ are the eigenvalues of $\beta'$. Then the entanglement
entropy \cite{Srednicki:1993im} corresponding to $(N-\mathbbm{n})$ oscillators
turns out to be $S=\sum_jS(\xi_{j})$, with $S(\xi)$ given by equation
(\ref{eq:ent-entropy-2HO-gen}) and $\xi_{j}^{2} =
\beta'_{j}/[1+(1-{\beta'}_{j}^2)^{1/2}]$.

\subsection{Entanglement entropy for free scalar field and area law}

In order to discuss about the area law for entanglement entropy, a free massive
scalar field $\Phi(\vv{x})$ is considered with mass $\mu$ and conjugate momenta
$\Pi(\vv{x})$. In Minkowski spacetime the Hamiltonian \cite{book:Peskin,
book:ADas, padmanabhan2016quantum} corresponding to the scalar field is
\begin{equation}\label{eq:Field-hamiltonian-massive}
 H = \frac{1}{2}\int d^3x~\left[
\Pi^2(\vv{x})+|\nabla\Phi(\vv{x})|^2 + \mu^2 \Phi^2(\vv{x}) \right]~.
\end{equation}
In terms of partial Fourier decomposition the field $\Phi(\vv{x})$ and the
conjugate momentum $\Pi(\vv{x})$ are transformed with respect to the angular
coordinates as
\begin{eqnarray}\label{eq:Field-fourier-transform}
 \Phi(\vv{x}) &=& \sum_{l,m}
Z_{lm}(\theta,\phi)~\frac{\varphi_{lm}(r)}{r}~,\nonumber\\
\Pi(\vv{x}) &=& \sum_{l,m}
Z_{lm}(\theta,\phi)~\frac{\pi_{lm}(r)}{r}~,
\end{eqnarray}
where $Z_{lm}(\theta,\phi)$ denotes real spherical harmonics and the Fourier
field modes satisfy a commutation relation $[\varphi_{lm}(r),\pi_{l'm'}(r')] =
i\delta_{ll'}\delta_{mm'}\delta(r-r')$ among themselves. With this definition of
field decomposition from Eqn. (\ref{eq:Field-fourier-transform}) the
Hamiltonian now becomes $H=\sum_{lm}H_{lm}$, where
\begin{eqnarray}\label{eq:Hamiltonian-massive-FT}
H_{lm} &=&  \frac{1}{2}\int_{0}^{\infty} dr \left\{
\pi^2_{lm}(r)+r^2\left[\frac{\partial}{\partial
r}\left(\frac{\varphi_{lm}(r)}{r}\right)\right]^2
\right.\nonumber\\
~&& ~~~~~~~~~~~~~\left. +\left(\frac{l(l+1)}{r^2}+\mu^2\right) \varphi^2_ {
lm}(r) \right\}~.
\end{eqnarray}
Next the radial coordinate $r$ is discretized forming a lattice with
inter-atomic spacing $a$ and size $L=(N+1)a$. The inverse of the spacing
$a^{-1}$ signifies the ultraviolet cutoff while the inverse system size $L^{-1}$
denotes infrared cutoff. This discretization makes the Hamiltonian look like a
set of coupled harmonic oscillators
 \begin{eqnarray}
H_{lm} &=&  \frac{1}{2a} \sum_{j=1}^{N} \left[
\pi^2_{lm,j}+\left(j+\frac{1}{2}\right)^2\left(\frac{\varphi_{lm,j}}{j}
-\frac{\varphi_ {lm,j+1}}{j+1}
\right)^2 \right.\nonumber\\
~&& ~~~~~~~~~~~~~~~~\left.
+\left\{\frac{l(l+1)}{j^2}+\mu^2a^2\right\}\varphi^2_{
lm,j} \right]~,\label{eq:Hamiltonian-massive-FT-discrete1}
\end{eqnarray}
such that $\varphi_{lm,N+1}=0$ and $[\varphi_{lm,j},\pi_{l'm',j'}] =
i\delta_{ll'}\delta_{mm'}\delta_{jj'}$. Comparison of this Hamiltonian with the
Hamiltonian for $N-$coupled harmonic oscillators from equation
(\ref{eq:general-Ncoupled-HO-hamiltonian}) gives \cite{Das:2008sy}
\begin{eqnarray}
 K_{jk} &=& \frac{1}{j^2} \left[
l(l+1)\delta_{jk}+\frac{9}{4}\delta_{j1}\delta_{k1}+\left(N-\frac{1}{2}
\right)^2\delta_{jN}\delta_{kN} \right.\nonumber\\
&+& \left. (\mu a)^2j^2\delta_{jk} +
\left\{\left(j+\frac{1}{2}\right)^2+\left(j-\frac{1}{2}\right)^2\right\}\delta_{
jk(j\ne 1, N)} \right]\nonumber\\
&-&\left[\frac{\left(k+\frac{1}{2}\right)^2}{k\left(k+1\right)}\right]\delta_{j,
 k+ 1}
-\left[\frac{\left(j+\frac{1}{2}\right)^2}{j\left(j+1\right)}\right]\delta_{j,
k- 1}~.\label{eq:Interaction_matrix_discrete_hamiltonian}
\end{eqnarray}
The discretization of the radial coordinate enables one to get a finite
expression of matrix $K$ denoting the potential energy and interaction. This in
turn would enable one to obtain the entanglement entropy when a finite number
$\mathbbm{n}$ of spatial points are traced out in total $N+1$ points. Then as
one plots the entanglement entropy with respect to $(\mathbbm{n}+1/2)^2$, one
gets a straight line which represents the celebrated area law for entanglement
entropy \cite{Srednicki:1993im}. We shall present the area curve of entanglement
entropy coming from perturbative formulation along with the curve obtained from
this usual formulation together in the next section.

\section{Entanglement entropy in perturbative approach}\label{EE_perturbed}

As discussed in previous section using Gaussian ground state-wave function of
coupled harmonic oscillators from Schrodinger quantization Eqn.
(\ref{eq:GS-wavefn-normal-freq}), one can easily evaluate the exact form of
reduced density matrix in Eqn. (\ref{eq:reduced-general}). However, in polymer
quantization, a quantization method used in \emph{loop quantum gravity}, the
ground state wave function is obtained in terms of Mathieu functions. To the
best of our knowledge evaluation of exact analytical form of reduced density
matrix is not possible even for two coupled oscillators using these polymer
wave-functions. This constraint further debars one to obtain the eigen-values
for $N-$coupled oscillators in polymer quantization and motivates us to take
help of the perturbation technique.

In this section we are going to apply perturbation to describe entanglement
entropy of coupled harmonic oscillators. We express the Hamiltonian
corresponding to the coupled oscillators in terms of a non interacting free
Hamiltonian $H_{0}$ and a net interaction term $\lambda H_{int}$ as
\begin{equation}\label{eq:Hamiltonian_interaction}
 H = H_0 +\lambda H_{int}~,
\end{equation}
When interaction strength is smaller than the strength of free Hamiltonian,
which is obtained for small $\lambda$, one can express the ground state
$|\Omega\rangle$ corresponding to the whole system in a perturbative manner as
\begin{equation}
 |\Omega\rangle=|0\rangle + \lambda |0^{1}\rangle + \lambda ^2 |0^{2}\rangle +
\dots~,\label{eq:GS-general-perturbation}
\end{equation}
where $|0\rangle$ denotes the ground state corresponding to the non-interacting
Hamiltonian $H_0$. On the other hand $|0^{1}\rangle$ and $|0^{2}\rangle$ denotes
the first and second order perturbative corrections to the non-interacting ground
state. From the time independent perturbation theory \cite{sakurai1985modern,
griffiths2018introduction, cohen1991quantum} one obtains the first order
correction to the ground state as
\begin{equation}
 |0^{1}\rangle = \sum_{n\neq 0}  \frac{\langle
n|\hat{H}_{int}|0\rangle}{E_0-E_n} |n\rangle ~,
\label{eq:perturbative-gen-wavefn1}
\end{equation}
where $E_n$ denotes the energy of the $n^{th}$ excited state corresponding to
the non interacting Hamiltonian. The second order correction to the ground state
is expressed as
\begin{eqnarray}
|0^{2}\rangle &=&- \sum_{n\neq 0}  \frac{\langle
0|\hat{H}_{int}|0\rangle\langle
n|\hat{H}_{int}|0\rangle}{(E_0-E_n)^2} |n\rangle \nonumber\\
&& +~\sum_{m\neq 0} \sum_{n\neq 0}  \frac{\langle
m|\hat{H}_{int}|n\rangle\langle
n|\hat{H}_{int}|0\rangle}{(E_0-E_m)(E_0-E_n)} |m\rangle ~.
\label{eq:perturbative-gen-wavefn2}
\end{eqnarray}
We shall use these perturbative corrections to obtain the actual ground state
upto certain perturbative order in the system of coupled harmonic oscillators.
We mention that while discussing polymer quantization \cite{Hossain:2010eb} we
shall consider only the $\pi-$periodic sector for our calculations. In
$\pi-$periodic sector, except for the ground state, the even and odd energies
become degenerate in high energy regimes. Now as we are interested in the ground
state density matrix it is convenient for us to consider the non-degenerate
perturbation theory.

\subsection{Entanglement entropy for two coupled harmonic oscillators}

We begin with a system of two coupled harmonic oscillators. We recall the
Hamiltonian from equation (\ref{eq:hamiltonian-full1}) and observe that it can
be expressed in form of Eqn. (\ref{eq:Hamiltonian_interaction}) with $H_0= H_1 +
H_2$ and $\lambda H_{int} = -k_1^2 x_1x_2$, where $H_j = [p_j^2+\omega^2
x_j^2]/2$  and $\omega=(\omega_{0}^2+k_1^2)^{1/2}$. Perturbative methods can be
applied when $k_1^2$ is smaller than $\omega^2$, which is always true for
nonzero $\omega_{0}$. Then in this system of two coupled oscillators the
correction to the ground state wave-function due to  first order perturbation
would be

 \begin{eqnarray}
\lambda |0^{1}\rangle = A_{nn}~ |n\rangle_{1} \otimes |n\rangle_{2} = A_{nn}~
|n,n\rangle ~,\label{eq:perturbative-wavefn1}
 \end{eqnarray}
where in the second compact notation of the wave-function the first index
corresponds to first oscillator and the second one corresponds to second
oscillator. Here the operation of $\hat{x}_j$ on the corresponding ground state
is given by
\begin{equation}
 \hat{x}_j |0\rangle_{j} = \sum_nC^j_{0n}|n\rangle_{j}~,\label{eq:operator_x}
\end{equation}
where in general the most dominating term comes from a single excitation
$|n\rangle_{j}$. Then we get for two coupled oscillators
\begin{equation}\label{eq:perturbative_coefficient}
 A_{nn} = \frac{k_1^2}{E_{n,n}-E_{0,0}}C^1_{0n}C^2_{0n}~.
\end{equation}
Considering up to the $1^{st}$ order perturbation, the normalized ground
state will be $ |\Omega\rangle=N_{nn}^{1}\left[|0,0\rangle + A_{nn}~
|n,n\rangle\right]$, where the normalization factor
$N_{nn}^{1}=(1+A_{nn}^2)^{-1/2}$. The corresponding reduced density matrix for
the first oscillator would be
\begin{equation}
 \hat{\rho}_{1} = \mathrm{Tr}\hspace{1pt}_{2}(|\Omega\rangle \langle\Omega|) =
(N_{nn}^{1})^2\left[
|0\rangle \langle 0| + A_{nn}^2 |n\rangle \langle n| \right]~,
\end{equation}
where the states now correspond to the first oscillator. This reduced density
matrix has eigen-values $(N_{nn}^{1})^2$ and $(N_{nn}^{1}A_{nn})^2$, and it
would give the entanglement entropy
\begin{equation}
 S_{E}^1 = -(N_{nn}^{1})^2\left[ \ln{\left(N_{nn}^{1}\right)^2} + A_{nn}^2
\ln{\left(A_{nn}N_{nn}^{1}\right)^2}\right]~.\label{eq:ent-entropy-1st-pert1}
\end{equation}
Now we consider $2^{nd}$ order perturbation and from Eqn.
(\ref{eq:perturbative-gen-wavefn2}) we observe that the first quantity would
vanish as $\langle 0|\hat{x}_1\hat{x}_{2}|0\rangle=0$, when discussing two
coupled oscillators. Then second order correction to ground state can be
expressed as
 \begin{eqnarray}
\lambda^2|0^{2}\rangle = A_{nn}k^2_{1}\sum_{m\neq0}
\frac{\langle m|\hat{x}_1\hat{x}_{2}|n,n\rangle}{E_{m}-E_{0,0}}~|m\rangle
~.\label{eq:perturbative-wavefn2}
 \end{eqnarray}
We shall evaluate this quantity explicitly in Schrodinger quantization and
compare the qualitative difference of resulting entanglement entropy with the
result obtained from first order perturbation.

\subsubsection{Schrodinger quantization}

\begin{figure}
  \includegraphics[width=0.8\linewidth]{com_EC_PD.pdf}
  \caption{(a) We have plotted the entanglement entropy of two coupled harmonic
oscillators with respect to varying $k_{1}$ in unit of $\omega_{0}$, from both
actual and perturbative formulation. (b) The percentage difference of
perturbative entanglement entropies from actual entropy is plotted, with respect
to varying $k_{1}$ in unit of $\omega_{0}$.}
  \label{fig:Perturbative-EE-comparison}
\end{figure}

In Schrodinger quantization $C^j_{0n}=\delta_{1,n}/\sqrt{2\omega}$ and
$E_n=(n+1/2)\omega$, then we have $\hat{x}_j|0\rangle_{j} =
1/\sqrt{2\omega}~|1\rangle_{j}$ and $E_{11}-E_{00}=2\omega$. Then
$A_{11}=A=k_1^2/4\omega^2$ and $\lambda |0^{1}\rangle = k_1^2/4\omega^2~
|1,1\rangle $. The entanglement entropy would be given by Eqn.
(\ref{eq:ent-entropy-1st-pert1}). When the interaction $k_1^2$ is very small
compared to the frequency square $\omega^2(\approx\omega_0^2, ~\textup{in~
that~ case})$ of the individual oscillator, the expression of entanglement
entropy from Eqn. (\ref{eq:ent-entropy-1st-pert1}) can be simplified to
\begin{equation}
 S_{E}^1 \approx - \ln{\left[1-\left(\frac{k_1}{2\omega_0}\right)^4\right]} -
\left(\frac{k_1}{2\omega_0}\right)^4
\ln{\left[\left(\frac{k_1}{2\omega_0}\right)^4\right]}~.\label{
eq:ent-entropy-1st-pert2 }
\end{equation}
One can observe that in similar conditions an exactly same expression for
entanglement entropy is obtained from equation (\ref{eq:ent-entropy-2HO-gen})
using (\ref{eq:eigenvalue-gen}). Thus at least for two coupled oscillators, when
the interaction is comparatively much lower than the frequency square, the first
order perturbation produces reasonable results in accordance with the results
from actual formulation. This fact can also be verified from FIG.
\ref{fig:Perturbative-EE-comparison}$:(a)$. Similarly in Schrodinger
quantization the second order correction to the wave-function from Eqn.
(\ref{eq:perturbative-wavefn2}) becomes
 \begin{eqnarray}
\lambda^2|0^{2}\rangle = A^2\left[|2,2\rangle +
\sqrt{2} \left(|0,2\rangle + |2,0\rangle \right)\right]
~.\label{eq:perturbative-wavefn22}
 \end{eqnarray}
Then the normalized ground state wave-function would be
\begin{small}
\begin{eqnarray}
 |\Omega\rangle &=& N_{2}\left[|0,0\rangle + A~ |1,1\rangle \right.\nonumber\\
 ~&+& \left. A^2\left\{|2,2\rangle + \sqrt{2} \left(|0,2\rangle + |2,0\rangle
\right)\right\}\right]~,
\end{eqnarray}
\end{small}
where $N_{2}=(1+A^2+5A^4)^{-1/2}$ is the normalization constant. One obtains
the reduced density matrix corresponding to the first oscillator as
\begin{eqnarray}
 \hat{\rho}_{1} &=& N_{2}^2\left[ (1+2A^4)|0\rangle \langle 0| + A^2 |1\rangle
\langle 1| + 3A^4 |2\rangle \langle 2| \right.\nonumber\\
&&~+ ~ \left. \sqrt{2}A^2(1+A^2) \left\{|2\rangle \langle 0| + |0\rangle \langle
2|\right\} \right]~.\label{eq:reduced-perturb1}
\end{eqnarray}
This reduced density matrix has eigenvalues
\begin{eqnarray}
\lambda_{1} &=& N_{2}^2A^2~\nonumber\\
\textup{and}\nonumber\\
\lambda_{2,3} &=& \frac{N_{2}^2}{2}\left[
1+5A^4\pm\sqrt{(1+A^2)^2(9A^4-2A^2+1)}
\right]~,\nonumber\\
\end{eqnarray}

which would give the entanglement entropy to be $S_{E} = -\sum_{s=1}^3
\lambda_s\ln{\lambda_s}$ . This entanglement entropy and the entanglement
entropy obtained from first order perturbation are plotted with the actual
entropy in FIG. \ref{fig:Perturbative-EE-comparison}$:(a)$ with respect to a
varying coupling between the two oscillators. The percentage difference of the
obtained result using perturbative techniques from the actual entanglement
entropy is plotted in FIG. \ref{fig:Perturbative-EE-comparison}$:(b)$. From
these figures we observe that when the coupling is small compared to the
individual frequency square of the oscillators, perturbation method is quite
elegant to study entanglement entropy of coupled harmonic oscillators.
Furthermore from these figures we also observe that the results from second
order perturbation does not drastically improve compared to first order
perturbation. On the other hand as our main objective is to understand the
qualitative nature of entanglement entropy from perturbation, it is expected
that first order perturbation would be good enough to satisfy our requirement.

\subsubsection{Polymer quantization}

In this part we are going to discuss about entanglement entropy for the system
of two coupled harmonic oscillators in polymer quantization. Perturbation
techniques will be used to obtain the entanglement entropy as the wave-functions
arising from polymer quantization cannot be handled analytically like the
Gaussian wave-functions. Here we start with a brief overview of the technical
aspects of polymer quantization.

Polymer quantization \cite{Hossain:2010eb} is a background independent
quantization procedure arising from loop quantum gravity (LQG). In polymer
quantization apart form the Planck constant $\hbar$ a new dimension-full
parameter $\lambda$ is introduced. Here the elementary operators are
configuration operator $\hat{x}$ and translation operator $\hat{U}_{\lambda}
\equiv \hat{e^{i\lambda p}}$ and their actions are defined as
\begin{equation}\label{eq:poly_operator}
 \hat{x} e^{i p x_{k}} = i\frac{\partial}{\partial p} e^{i p x_{k}}~,~~~~
 \hat{U}_{\lambda} e^{i p x_{k}} =  e^{i p (x_{k}+\lambda)}~.
\end{equation}
These operators satisfy the basic commutator $[\hat{x}, \hat{U}_{\lambda}] =
-\lambda \hat{U}_{\lambda}$. Now with the definition of translation operator
from Eqn. (\ref{eq:poly_operator}) and the inner product
\begin{equation}\label{eq:poly_innerP}
 \langle x|x' \rangle = \lim_{T\to\infty}\frac{1}{2T} \int_{-T}^{T} dp e^{-i p
x}e^{i p x'} = \delta_{x,x'}~,
\end{equation}
it is observed that a momentum operator cannot be defined as the translation
operator is not continuous in its parameter. However to describe the kinetic
energy part of the Hamiltonian one must have a suitable expression of the
momentum operator. In this case the momentum operator should be $\lambda$
dependent and to be given in terms of the translation operator. One simple
definition of the momentum operator as considered in \cite{Hossain:2010eb} is
\begin{equation}\label{eq:poly_pdef}
 \hat{p}_{\lambda} =
\frac{1}{2i\lambda}(\hat{U}_{\lambda}-\hat{U}^{\dagger}_{\lambda})~.
\end{equation}
One can then express the eigen-value equation $\hat{H}\psi=E\psi$, where
$H=p^2/2m+m\omega^2x^2/2$ represents the Hamiltonian corresponding to a simple
harmonic oscillator with mass $m$, as
\begin{equation}\label{eq:poly_eom}
 \psi''(u) + \left[ \sigma-\frac{1}{2g^2}\cos{(2u)} \right]\psi(u) = 0~,
\end{equation}
which represents a Mathieu equation \cite{Abramowitz1964handbook}. Here
$\lambda=\lambda_{\star}$, $u=\lambda_{\star}p+\pi/2$,
$g=m\omega\lambda_{\star}^2$ and $\sigma = 2E/g\omega-1/2g^2$. The above
differential equation has periodic solutions for $\sigma$ representing the
Mathieu characteristic value functions
\begin{eqnarray}\label{eq:poly_WF}
 \psi_{2n}(u) &=& \pi^{-1/2} ce_{n}(1/4g^2,u),
~~~\sigma=\mathcal{A}_{n}(g)\nonumber\\
 \psi_{2n+1}(u) &=& \pi^{-1/2} se_{n+1}(1/4g^2,u),
~~~\sigma=\mathcal{B}_{n}(g)~.\nonumber\\
\end{eqnarray}
For $(n=0,1,...)$, $ce_{n}$ and $se_{n}$ represent the elliptic cosine and
sine functions, where for even $n$ they are $\pi-$periodic and for odd $n$ they
are $\pi-$antiperiodic functions. The corresponding energy eigen values are
given by
\begin{eqnarray}\label{eq:poly_energy}
 E_{2n} &=& \frac{\omega}{4g}\left[2g^2\mathcal{A}_{n}(g)+1\right]\nonumber\\
 E_{2n+1} &=& \frac{\omega}{4g}\left[2g^2\mathcal{B}_{n+1}(g)+1\right]~.
\end{eqnarray}
Using the asymptotic expansions of the Mathieu characteristic value functions
$\mathcal{A}_{n}(g)$ and $\mathcal{B}_{n}(g)$ one can get in small $g$ limit,
i.e. when $g\ll1$
\begin{eqnarray}\label{eq:poly_energy_gL}
\frac{E_{2n}}{\omega}\approx\frac{E_{2n+1}}{\omega}=\left(n+\frac{1}{2}
\right)-\frac{(2n+1)^2+1}{16}g+\mathcal{O}(g^2)~.\nonumber\\
\end{eqnarray}
On the other hand in high $g$ regimes, i.e. when $g\gg1$, one gets
\begin{eqnarray}\label{eq:poly_energy_gH}
\frac{E_{0}}{\omega}=\frac{1}{4g}+\mathcal{O}\left(\frac{1}{g^3}\right),~\frac{
E_ { 2n-1}}{\omega} \approx\frac{E_{ 2n}}{\omega} =
\frac{n^2g}{2}+\mathcal{O}\left(\frac{1}{g}\right).\nonumber\\
\end{eqnarray}
From the asymptotic expression (\ref{eq:poly_energy_gL}) we observe that when
$g\ll1$ the energy levels corresponding to the $\pi-$periodic and
$\pi-$antiperiodic sectors becomes degenerate among themselves. On the other
hand from (\ref{eq:poly_energy_gH}) we observe that for $g\gg1$ the energy
levels within the separate $\pi-$periodic and $\pi-$antiperiodic sectors becomes
degenerate. We shall consider only the $\pi-$periodic sector of the
wave-functions, containing the non-degenerate ground state, to discuss about the
corresponding entanglement entropy. We want to mention that the asymptotic
expressions of the energy eigen-values from Eqn. (\ref{eq:poly_energy_gL},
\ref{eq:poly_energy_gH}) are also utilized in \cite{Hossain:2014fma,
Hossain:2015xqa, Barman:2017vqx} to observe the Unruh and Hawking effect for
polymer observer. One can also look into \cite{Kajuri:2015oza,
Stargen:2017xii, Louko:2017emx, Kajuri:2017zmr} where polymer quantization is
used in different systems to study particle creation.
\vspace{0.5cm}

\emph{Entanglement entropy in polymer quantization:}
In this part we evaluate the perturbative corrections to the ground state in
polymer quantization, which basically requires the estimation of $A_{nn}$. The
operation $\hat{x}_j|0\rangle_j$ is already discussed in Eqn.
(\ref{eq:operator_x}) and in polymer quantization $C^j_{0n}$ are given by
\begin{eqnarray}
  C^j_{0n} &=& {}_{j}\langle
n|\hat{x}_{j}|0\rangle_{j}=\lambda_{\star}\int_0^{2\pi}
\psi^j_n ~\frac{\partial}{\partial u^j}\psi^j_0~ du^j
~,\label{eq:Cn_general_poly}
\end{eqnarray}
where $u^j=\lambda_{\star}p^j+\pi/2$ and $\lambda_{\star}$ is the polymer length
scale, see \cite{Hossain:2010eb}. There are infinite number of non-zero
$C^j_{0n}$ in polymer quantization, where as in Schrodinger quantization there
is only one $C^j_{0n}=\delta_{1,n}/\sqrt{2\omega}$. In order to compute polymer
corrections we only consider the first and most dominating non-zero $C^j_{0n}$,
which is $C^j_{03}$. In small $g=\omega \lambda_{\star}^2$ limit, i.e. when
$g\ll1$, these coefficients are given by
\begin{equation}
C^j_{03}= C_{03}=\frac{i}{\sqrt{2\omega}}[1-\frac{3}{4} g+\mathcal{O}(g^2)]~,
\end{equation}
and the corresponding energy correction is given by
\begin{equation}
 \Delta
E_{30}=E_3-E_0=\omega[1-\frac{g}{2}+\mathcal{O}(g^2)]~.
\end{equation}
The expression of $A_{nn}$, now obtained as $A_{33}$, is changed and using Eqn.
(\ref{eq:perturbative_coefficient}) becomes
\begin{eqnarray}
 A_{33} = A&\approx& -\frac{k_1^2}{4\omega^2(1-g/2)}(1-\frac{3}{2}g)\nonumber\\
 &\approx&-\frac{k_1^2}{4\omega^2}(1-g)~.
\end{eqnarray}
One can observe from this expression of $A$ that as one takes $g\to 0$, one gets
back the result from Schrodinger quantization. We note that the sign of $A$ do
not affect the end result as the entanglement entropy is obtained using $A^2$.
In the ultraviolet limit when $g\gg1$, one has the expressions
\begin{equation}
 C^j_{03}=C_{03} = i\sqrt{\frac{g}{2\omega}}\left[\frac{1}{4g^2} +
\mathcal{O}\left(\frac{1}{g^6}\right)\right]~,
\end{equation}
and
\begin{equation}
 \Delta
E_{30}=\omega\left[2g+\mathcal{O}\left(\frac{1}{g^3}\right)\right]~.
\end{equation}
Then the expression of $A_{33}$ is
\begin{equation}
 A_{33}=A\approx -\left(\frac{k_1}{2\omega}\right)^2 \frac{1}{2(2g)^4}~.
\end{equation}
We know that the reduced density matrix in first order perturbation has
eigenvalues $\lambda_{1}=N_{1}^2=(1+A^2)^{-1}$ and $\lambda_{2} = (N_{1}A)^{2}$.
Then for fixed $(k_1/\omega)$ as we take $g\to\infty$, we observe that
$\lambda_{1}\to1$ and $\lambda_{2}\to0$, because in this limit $A\to0$. Then the
entanglement entropy evaluated from these eigenvalues would vanish, providing a
very new feature in ultraviolet regime of energy in polymer quantization.

\begin{figure}
  \includegraphics[width=0.8\linewidth]{Poly_EE_fixed.pdf}
  \caption{(a) The entanglement entropy of two coupled harmonic oscillators
using perturbative formulation in polymer quantization, with varying $g$ keeping
$\alpha=k_1/\omega$ fixed. (b) The $Log-Log$ plot of the entanglement entropy in
polymer quantization, with respect to varying $g$ keeping $\alpha=k_1/\omega$
fixed.}
\label{fig:polymer-EE}
\end{figure}

We now intent to express this result in a more general fashion without using
asymptotic forms of the Mathieu functions. For this we shall need some
numerical help. First the expression of general energy difference of our
concern in polymer quantization is
\begin{equation}
 E^{0}_{33}-E^{0}_{00} = \omega g[\mathcal{B}_{2}(g)-\mathcal{A}_{0}(g)]~,
\label{eq:dEn_general_poly}
\end{equation}
where $\mathcal{A}_{n}(g)$ and $\mathcal{B}_{n}(g)$ are the Mathieu
characteristic value functions corresponding to even and odd Mathieu functions.
The expressions of $C_{03}$ are obtained from Eqn. (\ref{eq:Cn_general_poly})
with wave-functions represented in terms of Mathieu functions. Then we have
\begin{equation}\label{eq:perturbative_coefficient_poly}
 A=A_{33} = \frac{k_1^2}{\omega
g[\mathcal{B}_{2}(g)-\mathcal{A}_{0}(g)]}C^1_{03}C^2_{03}~.
\end{equation}
and the corresponding eigen-values are
\begin{equation}
 \lambda_{1}=N_{1}^2=\left[1 + \left\{\frac{(\alpha
~C_{03}/\lambda_{\star})^2}{\mathcal{B}_{2}(g)-\mathcal{A}_{0}(g)}\right\}^2
\right]^{-1}~,
\end{equation}
and
\begin{equation}
 \lambda_{2} = \left\{N_{1}~\frac{(\alpha
~C_{03}/\lambda_{\star})^2}{\mathcal{B}_{2}(g)-\mathcal{A}_{0}(g)}\right\}^2 ~,
\end{equation}
where $\alpha =k_1/\omega$. We have plotted this entanglement entropy coming
from first order perturbation for different fixed values of $\alpha$ with
varying polymer parameter $g$, in FIG. \ref{fig:polymer-EE}. Then the change of
$g$ signifies the change in harmonic oscillator frequency $\omega$ for the fixed
ratio $k_1/\omega$ and fixed polymer length scale $\lambda_{\star}$. In these
plots the fixed ratio $(k_1/\omega)$ is considered to be less than one. It
implies that the interaction strength is less than one and permits the
application of perturbation theory even for high $g$ regimes. We observe that at
high frequency regime as $g$ increases the entanglement entropy decreases and
becomes very low at large $g$, see FIG. \ref{fig:polymer-EE}. This situation was
not present in Schrodinger quantization as there the entanglement entropy is a
function of $A=k_1^2/4\omega^2$ and its value is fixed for fixed $\alpha
=k_1/\omega$. We want to note that same phenomena can also be observed in
polymer quantization using second order perturbation with very little
quantitative difference. For a unit mass harmonic oscillator one can interpret
the inverse square-root of it's frequency to be a length scale characteristic of
the harmonic oscillator. Now as the frequency increases this length decreases
and even reaches polymer length scale $\lambda_{\star}$ when $g$ becomes very
high. One then interprets the above deviation of entanglement entropy in polymer
quantization from usual quantization, as a result of the physics in very high
energy or in a very small length scale addressed by polymer quantization.

\subsection{Area law for free scalar field}

In this part we are going to use perturbation technique to evaluate the
entanglement entropy corresponding to a massive free scalar field described by
the Hamiltonian (\ref{eq:Field-hamiltonian-massive}). We first provide a
prescription for the eigen-values of the reduced density matrix
corresponding to $N-$coupled oscillators. By considering only first order
perturbation, one can express the ground state wave-function of $N$ weakly
coupled harmonic oscillators as
\begin{eqnarray}
|\Omega\rangle_N &=& |000...\rangle+A_1|nn0...\rangle+
...+A_{N-1}|...0nn\rangle,\nonumber\\\label{eq:Gen_perturb_waveFn}
\end{eqnarray}
where $A_1,...,A_{N-1}$ are coefficients of the first order perturbative
correction to wave-function and they are functions of the individual frequency
and interaction between different oscillators. In both Fock and polymer
quantization the state $|00...n_{j}n_{j+1}...\rangle$ is obtained from
$\hat{x}_{j}\hat{x}_{j+1}|00...\rangle=C^j_{0n} C^{j+1}_{0n} |00...n_{j}n_{j+1}
...\rangle$, where we have omitted the sum on $n$ as the most dominating
contribution comes from a single term. We first include appropriate
normalization factor to the wave function from Eqn.
(\ref{eq:Gen_perturb_waveFn}) and calculate the corresponding density matrix and
its eigen-values for successively increasing number of coupled harmonic
oscillators. Then the eigen-values, corresponding to the reduced density matrix
after tracing over $\mathbbm{n}-$degrees of freedom out of total $N-$coupled
harmonic oscillators, can be found by guessing from these consecutive
eigen-value evaluation, as
\begin{eqnarray}
 \lambda^{\mathbbm{n}}_{s} &=& N^{2}_{c}A^2_{\mathbbm{n}},
~\frac{N^{2}_{c}}{2}\left[\sum_{\substack{j=0 \\
j\ne \mathbbm{n}}}^{N-1}A^2_{j}\pm \left\{\left(\sum_{\substack{j=0 \\ j\ne
\mathbbm{n}}}^{N-1}A^2_{j}\right)^2
\right.\right.\nonumber\\
~&& ~~~~~~~\left.\left. -~~4\left(\sum_{\substack{j=1 \\ j\ne
\mathbbm{n}}}^{\mathbbm{n}-1}A^2_{j}
\right)\left(\sum_{\substack{j=\mathbbm{n}+1 \\ j\ne \mathbbm{n}}}^{N-1}A^2_{j}
\right)\right\}^{\frac{1}{2}}
\right],\label{eq:eigenvalues-Ncoupled-perturbed-gen}
\end{eqnarray}
where $A_{0}=1$ and subscript $s$ denotes different eigenvalues which are three
in number for any particular reduction $\mathbbm{n}$. $N_{c} = (\sum_{j=0}^{N-1}
A^{2}_{j})^{-1/2}$ is the normalization factor corresponding to the perturbed
ground state. As discussed earlier the entanglement entropy corresponding to
these eigenvalues would be $S_{E}^{\mathbbm{n}} = -\sum_{s=1}^{3}
\lambda^{\mathbbm{n}}_s \ln{\lambda^{\mathbbm{n}}_s}$.

\subsubsection{Area law in Fock quantization}

\begin{figure}
  \includegraphics[width=0.8\linewidth]{com_Fock_Area.pdf}
  \caption{(a) The plot of entanglement entropy with respect to
$(\mathbbm{n}+0.5)^2$ for a free massless scalar field in a lattice of total
size $N=60$. Here $\mathbbm{n}$ is the number of position degrees of freedom
traced out of the ground state of the scalar field. This plot shows the area law
for entanglement entropy in both actual and perturbative formulation. Here the
ratio of the slopes from perturbative entropy and actual entropy is $\sim 1.85$.
(b) This figure provides the $Log-Log$ plot of entanglement entropy vs
$(\mathbbm{n}+0.5)$ in both actual and perturbative formulation. The slope of
the curve from actual formulation is about $\sim1.9$ where the slope from
perturbative formulation is $\sim2.1$. We have used first order perturbation to
evaluate the entanglement entropy.}
  \label{fig:area-law-comparison}
\end{figure}

In order to obtain the area law of entanglement entropy in Fock quantization we
first consider the discrete Hamiltonian, formed out of partially Fourier
transformed field Hamiltonian in a lattice of finite size, from Eqn.
(\ref{eq:Hamiltonian-massive-FT-discrete1}). With the help of Eqn.
(\ref{eq:Interaction_matrix_discrete_hamiltonian}) one can get frequency
$\omega_j$ of the $j^{th}$ oscillator and coupling $k^2_j$ between $j^{th}$ and
$(j+1)^{th}$ oscillator. Now according to Eqn.
(\ref{eq:perturbative_coefficient}) we want to find the expressions of the
coefficients of first order perturbation $A_{j} = k^2_{j}/
\left\{2\sqrt{\omega_{j} \omega_{j+1}} (\omega_{j} + \omega_{j+1}) \right\}$,
which are used in Eqn. (\ref{eq:eigenvalues-Ncoupled-perturbed-gen}) to get the
eigen-values. We note that the perturbative coefficients $A_{j}$ are in
principle functions of $l$ and $m$ and sums over these quantities are taken to
evaluate the entropy, but for brevity we omitted their index from the notation.
In Fock quantization they are evaluated using $C^j_{0n} =
\delta_{1,n}/\sqrt{2\omega_j}$ and $E^j_n=(n+1/2)\omega_j$. We have numerically
computed the entanglement entropy using the obtained eigen-values. In FIG.
\ref{fig:area-law-comparison} the entanglement entropy from perturbative and
actual formulation are presented for a massless free scalar field.
The entanglement entropy in actual formulation is obtained non perturbatively,
utilizing the Gaussian nature of the ground state wave-functions. We have used
the results from Eqn. (\ref{eq:rho_nreduced}) and the potential from Eqn.
(\ref{eq:Interaction_matrix_discrete_hamiltonian}) to evaluate the entanglement
entropy in actual formulation \cite{Srednicki:1993im}.
FIG. \ref{fig:area-law-comparison} shows that first order perturbation is
sufficient enough to provide an area law for entanglement entropy. We want to
note that the slope from this area curve($\sim 0.42$) is different than the one
obtained from actual formulation($\sim 0.29$).
One can notice from Eqn. (\ref{eq:Hamiltonian-massive-FT-discrete1},
\ref{eq:Interaction_matrix_discrete_hamiltonian}) that the ratio of the
frequency square to interaction strength decreases when $j$ increases as we
consider a larger system size. Thus perturbation theory becomes less effective
and the results obtained from first order perturbation deviates more from the
actual non perturbative result for large $j$.
We also want to note that as second order perturbation is employed in this
formulation the area curve is quantitatively very little improved and gets
closer to the area curve from actual formulation. Now as we consider a field
with increasing mass, the ratio of the frequency square and coupling strength
from Eqn. (\ref{eq:Interaction_matrix_discrete_hamiltonian}) increases and the
perturbative formulation becomes more effective in describing the original
system.

\begin{figure}
    \includegraphics[width=0.8\linewidth]{FArea_mass24.pdf}
    \caption{(a) The plot of entanglement entropy of a massive free scalar field
with respect to $(\mathbbm{n}+0.5)^2$ for $N=60$. Here $\mu a=2$ and the ratio
of the slopes from perturbative entropy and actual entropy is $\sim 1.22$. (b)
The entanglement entropy vs $(\mathbbm{n}+0.5)^2$ plot for massive free scalar
field with $N=60$. Here $\mu a=4$ and the ratio of the slopes from perturbative
entropy and actual entropy is $\sim 1.04$. In both cases only first order
perturbation is used.}
\label{fig:area-law-comparison-mu24}
\end{figure}

In addition to the massless case we have also considered massive scalar field
with discretized Hamiltonian from Eqn.
(\ref{eq:Hamiltonian-massive-FT-discrete1}). We have taken different values for
the parameter $\mu a$ and for each value obtained the area law in both the
actual and first order perturbation. Then the massless case becomes a special
case when this parameter $\mu a$ is taken to be zero. In particular we have
observed that as we increase the value of this parameter the perturbation
becomes stronger and the slopes of the two area curves corresponding to actual
and perturbative formulation get closer. The plots corresponding to this
feature are shown in FIG. \ref{fig:area-law-comparison-mu24}.

\subsubsection{Area law in polymer quantization}

We take the Fourier Hamiltonian density for a massive free scalar field as
according to equation (\ref{eq:Hamiltonian-massive-FT-discrete1}). We note that
the field modes $\varphi_{lm,j}$ and their conjugate momentums $\pi_{lm,j}$ are
dimensionless here. To consistently introduce polymer quantization in this
formulation we make the transformations $\tilde{\pi}_{lm,j} \equiv \pi_{lm,j}/\sqrt{a}$
and $\tilde{\varphi}_{lm,j} \equiv \sqrt{a}~\varphi_{lm,j}$ such that the new field
mode and the conjugate momentum become dimension full. Then the Fourier
Hamiltonian density from Eqn. (\ref{eq:Hamiltonian-massive-FT-discrete1}) will
take the form
 \begin{eqnarray}
H_{lm} &=&  \tfrac{1}{2} \sum_{j=1}^{N} \left[
\tilde{\pi}^2_{lm,j}+\left(\frac{j+\frac{1}{2}}{a}\right)^2\left(\frac{\tilde{\varphi}_{lm,j
} }
{j}
-\frac{\tilde{\varphi}_ {lm,j+1}}{j+1}
\right)^2 \right.\nonumber\\
~&& ~~~~~~~~~~~~~\left. +\left\{\frac{l(l+1)}{a^2j^2}+\mu^2\right\}\tilde{\varphi}^2_{
lm,j} \right]~.\label{eq:fourier-hamiltonian-field-discrete3}
\end{eqnarray}
This Hamiltonian also describes a system of $N$ coupled harmonic oscillators
given by Eqn. (\ref{eq:general-Ncoupled-HO-hamiltonian}). In polymer
quantization a new dimension full parameter $\lambda_{\star}$ is introduced with
dimension $(length)^{1/2}$, inverse of the dimension of momentum. Here the basic
variables are taken to be $\tilde{\varphi}_{lm,j}$ and $U^{\lambda_{\star}}_{lm,j} =
\exp{\left\{i\lambda_{\star} \tilde{\pi}_{lm,j}\right\}}$ with Poisson bracket
$\{\tilde{\varphi}_{lm,j},U^{\lambda_{\star}}_{lm,j}\} = i\lambda_{\star}
U^{\lambda_{\star}}_{lm,j}$. From the above system of coupled harmonic
oscillators we observe that for a general $j^{th}$ oscillator the frequency is
\begin{eqnarray}
 \omega_{lm,j} &=& \frac{1}{a~j}\left[l(l+1)+\left\{\left(j+\tfrac{1}{2}
\right)^2+\left(j-\tfrac{1}{2}\right)^2\right\}\right|_{j\neq1,N}
\nonumber\\
~&& \left.+\mu^2 a^2 j^2+\tfrac{9}{4}\delta_{j1} +
\left(N-\tfrac{1}{2}\right)^2\delta_{jN}\right]^{\frac{1}{ 2 } }\nonumber\\
~&=& \frac{\Omega_{lm,j}}{a}~.
\end{eqnarray}

\begin{figure}
  \includegraphics[width=0.8\linewidth]{com_Poly_Area.pdf}
  \caption{(a) The entanglement entropy of a massless free scalar field with
respect to $(\mathbbm{n}+0.5)^2$ using first order perturbation in polymer
quantization. This plot shows the area law in polymer quantization with slope
$6.12\times10^{-6}$. (b) The $Log-Log$ plot of the entanglement entropy vs
$(\mathbbm{n}+0.5)$ in perturbative formulation using polymer quantization. The
total lattice size is taken to be $N=60$ and the obtained slope is $2.002$.}
  \label{fig:Poly-area-law}
\end{figure}

Now we want to get the expressions of perturbative coefficients $A_{j}$  used in
Eqn. (\ref{eq:eigenvalues-Ncoupled-perturbed-gen}). They are constructed using
expressions of $C^j_{0n}$ and $E^j_n$ in polymer quantization from Eqn.
(\ref{eq:Cn_general_poly}) and (\ref{eq:dEn_general_poly}), which are further
given by the dimension less polymer parameter $g_{lm,j} = \omega_{lm,j}
\lambda^2_{\star} = \Omega_{lm,j}l_{\star}/a$. We also want to note that the
inter atomic distance $a$ and the polymer length scale
$l_{\star}=\lambda^2_{\star}$ both have same nature and should have same order
as they signify the ultraviolet cutoff. Then we take their ratio
$\gamma=l_{\star}/a$ to be unity, which adds further simplification to the
evaluation of entanglement entropy. We have plotted the entanglement entropy
from first order perturbative formulation in FIG. \ref{fig:Poly-area-law}
considering massless free scalar fields in polymer quantum field theory. In
these figures we observe that the area law is valid in polymer quantization too.
However the corresponding slope is now very low compared to results from Fock
quantization. One can also get the area law in polymer quantization for massive
free scalar field with a further decreased slope.

\emph{Implication of the result:}
From \cite{Das:2008sy} we get to understand that the slope of the area curve for
entanglement entropy can be different due to many reasons, such as due to
different discretization procedures, inclusion of mass or taking excited states
instead  of the ground state. We want to mention here a consistency check to
understand whether this result from polymer quantization is a plausible one or
not. It is noted in \cite{Hossain:2010eb} that in low energy regimes polymer
quantization reproduces the results from usual Fock quantization. Now in this
formulation of entanglement entropy evaluation we observe that one direct
influence of polymer quantization over Fock is dictated by the factor
$(l_{\star}/a)$. When this factor is unity the system is completely interpreted
in terms of polymer quantization. On the other hand when the value of this
quantity decreases the value of the dimensionless polymer frequency $g_{lm,j}$
decreases and the system becomes more and more Fock like as the lower energy
regimes of polymer quantization tends to contribute to the description of the
system. We have plotted the entanglement entropy for different values of this
factor and we observed that as the value decreases the area curve of
entanglement entropy from polymer quantization approaches the one from Fock
quantization, see FIG. \ref{fig:Poly_Slope_Change} and FIG.
\ref{fig:log_Poly_Slope_Change}. Thus the very low slope of the entanglement
entropy can be described as a feature coming from the  disentangling nature of
polymer quantization at high energy regimes. We want to note that massive scalar
fields also show disentangling nature and lowers the slope of the area curve
\cite{Katsinis:2017qzh, Riera:2006vj, Balasubramanian:2011wt}.

\begin{figure}
  \includegraphics[width=0.8\linewidth]{Poly_Slope_Change.pdf}
  \caption{The plot of entanglement entropy in polymer quantization for
different values of $\gamma=l_{\star}/a$.}
  \label{fig:Poly_Slope_Change}
\end{figure}

\begin{figure}
  \includegraphics[width=0.8\linewidth]{loglog_Poly_Slope_Change.pdf}
  \caption{The $Log-Log$ plot of entanglement entropy in polymer quantization
for different values of $\gamma=l_{\star}/a$.}
  \label{fig:log_Poly_Slope_Change}
\end{figure}

The entanglement entropy of free scalar field in polymer quantization gives rise
to another question, which relates to the corrections to the area law as
predicted by quantum gravity \cite{Pasqua:2014hea, Sen:2012cj,
El-Menoufi:2015cqw, Pathak:2016vfc}. In this manner we want to note that the
slope of the entanglement entropy vs $(\mathbbm{n}+0.5)$ curve in \emph{Log-Log}
plot is $2.002$, which automatically discards any possible departure from the
area law. This area dependence of entanglement entropy in polymer quantization
is enthralling in its own right since it validates the generality of the area
law in quantizations other than Fock.

\section{Discussion}

In usual formulation, procurement of the area curve for entanglement entropy
\cite{Bekenstein:1974ax, Kallosh:1992wa, hawking1975, Bekenstein:1975tw,
Hawking:1976de, Gibbons:1976pt, Davies:1978mf} is simplified using the
mathematical structure of Gaussian ground state wave-function from Schrodinger
quantization. However not all quantization procedures provide this Gaussian
nature of ground state and polymer quantization is one of them.
We note that though entanglement entropy for two coupled harmonic oscillators
are specifically evaluated for polymer quantization in \cite{Demarie:2012tz},
the framework to obtain entanglement entropy for large number of coupled
oscillators is not provided thus one can not obtain the area law. In this work
we have treated the interaction between coupled harmonic oscillators in
perturbative manner. Our procedure is different than the ones discussed in
\cite{Katsinis:2017qzh, Riera:2006vj, Balasubramanian:2011wt, Kumar:2017ctm},
where the eigen-values of the reduced density matrix and momentum space
entanglement are estimated using perturbation.
For two coupled harmonic oscillators we noticed disentangling nature from
polymer quantization at high frequency regime. We observed that in Schrodinger
quantization the entanglement entropy is unchanged while in polymer quantization
it decreases at high oscillator frequencies, keeping the ratio of interaction
strength to frequency square fixed. We showed that in our formulation, by
considering free scalar field, one obtains the area law of entanglement entropy
for Fock quantization. As the mass of the scalar field increases the individual
oscillator frequency increases, thus perturbation strength increases and the
obtained area curve approaches the area curve from usual formulation.
Furthermore we showed that in polymer quantization also this formulation
provides a similar area law, but with a very decreased slope. We inferred that
this decrease of slope is due to the disentangling nature of polymer
quantization at higher energies. We further noticed that as the effect of
polymer quantization becomes smaller, by lowering the value of the ratio of
polymer length scale $l_{\star}$ to inter-atomic distance $a$, the area curve
from polymer quantization using first order perturbation tends to approach the
area curve from Fock quantization. This phenomena is not quite surprising as in
the limit $l_{\star}/a \rightarrow 0$, the physical result from polymer
quantization would converge to the result obtained from the standard Fock
quantization.
The disentangling nature of polymer quantization is very intriguing in its own
right as it is known that usual quantization looses its predictability in
trans-Planckian energy regimes \cite{Cohen:1998zx, Carmona:2000gd}.
We note that this disentangling phenomena in polymer quantization is analogous
to the suppression of propagation at large energies.
We mention that there are other derivations to obtain the area law and harvest
entanglement entropy for scalar field \cite{Calabrese:2004eu, Callan:1994py,
Calabrese:2005zw, Solodukhin:2011gn, VanRaamsdonk:2009ar, Casini:2009sr,
Cramer:2005mx, Page:1983ug, Kabat:1994vj, Plenio:2004he, Solodukhin:1994yz,
Holzhey:1994we, Das:2001ic, Allouche:2018err, Pagani:2018mke, Henderson:2017yuv,
Reznik:2002fz} and it would be interesting to see whether an exact form of the
entanglement entropy can be found using these derivations in polymer
quantization. In conclusion we address that our formulation opens up an avenue
to understand entanglement entropy in terms of perturbative corrections.\\

\begin{acknowledgments}
We would like to thank Golam Mortuza Hossain, Narayan Banerjee, Ritesh K. Singh
and Ananda Dasgupta for discussions. We would also like to thank Sumanta
Chakraborty, Abhishek Majhi and Sudipta Saha for useful suggestions. S.B. would
like to thank IISER Kolkata for supporting this work through a doctoral
fellowship.
\end{acknowledgments}