\section{Introduction}
\label{1}
Since gaining popularity in the ImageNet Large Scale Visual Recognition Competition (ILSVRC) \cite{b1} of 2012, convolutional neural networks (CNNs) have received much attention from various research fields. A CNN consists of several elements, including convolutional layers, pooling layers, and a fully connected layer. Since CNN first became famous in ILSVRC 2012, pooling has been included in all popular networks: AlexNet \cite{b2}, GoogleNet \cite{b3}, Visual Geometry Group (VGG) \cite{b4}, Residual Neural Network (ResNet) \cite{b5}, and DenseNet \cite{b6}.

Pooling plays several important roles in CNNs. First, it reduces the size of the feature map, enabling training and testing with a limited amount of computation. It also reduces the computation of the network. Second, pooling increases the size of receptive field size in each neuron, enabling each neuron to "see" a larger portion of the input image. This functionality improves the recognition performance. Third, pooling reduces the negative effects of noise and slight distortion. Pooling admits only the representative values and removes the small noises. In essence, pooling simultaneously improves the recognition performance and reduces the computational demands.

The deep learning performance depends on the choice of the pooling method. The most popular pooling methods are average pooling, max pooling, and stride pooling. Overall, these simple methods are efficient because they easily compute the representative values. However, they do not consider the various input patterns and there is still a room to exploit for improving the pooling. Recently, several new pooling approaches have been reported, including a linear combination of simple existing poolings \cite{b15}, and an adaptive pooling method that magnifies the spatial changes and preserves the important structural details \cite{b16}. Most of the new pooling methods assume fixed weights for feature training, which do not account for the characteristics of individual images. Therefore, the fixed-weight constraint on the pooling might degrade the performance of the CNN.

 On the other hands, researches on attention module \cite{b7, b8} have been actively conducted to improve the performance of the network. The attention concept weights a major feature in the feature map. Several studies have shown performance improvements after weighting appropriate parts of the feature maps. Attenuation modules  have been investigated in the spatial domain, time domain, and channel direction.

This paper proposes a new pooling method named universal pooling. The proposed pooling is the most general pooling and includes the previous simple poolings such as max pooling, average pooling and stride pooling as special cases. The basic idea of universal pooling is to interpret pooling as attention and extend it to the general channel-wise local spatial attention. That is, the universal pooling selects pooling weights for each channel and they are trained together with other feature extraction parts. Because it embraces the existing simple poolings, universal pooling is expected to outperform the independent simple poolings.
The remainder of this paper is organized as follows. Section \ref{2} reviews some related works, and Section \ref{3} introduces the preliminary fundamentals of pooling. Section \ref{4} proposes the universal pooling concept and discusses its relationship with the simple pooling approaches. Section \ref{5} applies the proposed pooling to two benchmark datasets and presents the experimental results. Conclusions are drawn in Section \ref{6}.

\section{Related Work}
\label{2}
The commonest pooling methods are max and average pooling, but stride pooling in which the feature map is reduced using the stride in the convolution layer has also gained traction in recent years.

Max pooling divides the feature map into blocks and collects the maximum feature value in each block into a smaller output matrix. Max pooling is commonly used between convolution layers and is employed in AlexNet \cite{b2} and VGG \cite{b4}. Average pooling operates similarly to max pooling, but outputs the average of each block in the feature map. Global average pooling (GAP), which applies average pooling over the entire feature map, is commonly used in convolutional networks such as ResNet \cite{b5} and DenseNet \cite{b6}. In DenseNet, average pooling is applied between the convolution layers. Meanwhile, stride pooling is equivalent to importing values from a fixed position after convoluting across the entire area. This approach is adopted in ResNet.

All of these pooling methods are efficient but simple, and it seems that there is a room to improve the performance. S3pool \cite{b9} and stochastic pooling \cite{b19} adopt a probability-based pooling approach. $L_p$  pooling of various coefficients' norms was proposed in \cite{b10} and \cite{b11}, and a fractional version of max pooling was proposed in \cite{b12}. The spectral space was down-sampled through a filter in \cite{b13} and \cite{b14}. In \cite{b15}, the existing simple pooling methods were combined to improve the pooling performance. Detail-preserving pooling, which preserves the feature details by applying existing down-sampling techniques in the image processing area and by learning the parameters of the function over the network, was proposed in \cite{b16}.

\section{Preliminary Fundamentals}
\label{3}
The pooling operation, which down-samples the feature map, can be considered as a multiple input function
\begin{equation}
\rho :{\Re ^{{S^2}}} \to \Re
\label{eq1}
\end{equation}
which returns only a scalar value. $S^2$ is the area of the pooling block $\left( S \times S \right)$. For example, let us consider a pooling of size 2 with a stride 2, as shown in Figure \ref{fig1}.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure1.pdf}}
\caption{Pooling downsizes the feature map. Each red block denotes one pooling block.}
\label{fig1}
\end{figure}

In this case, the pooling can be represented by a multiple input function
\begin{equation}
{o_{0,0}} = \rho ({f_{0,0}},{f_{0,1}},{f_{1,0}},{f_{1,1}}) .
\label{eq2}
\end{equation}
As mentioned above, there are three popular pooling methods: max pooling, average pooling, and stride pooling. In the following, these methods will be explained on a pooling block of size $S=2$ and stride 2.

\subsection{Max pooling}
\label{3A}
Max pooling computes the output using the maximum function. The output from the max pooling is defined as
\begin{equation}
{o_{0,0}} = \max ({f_{0,0}},{f_{0,1}},{f_{1,0}},{f_{1,1}})
\label{eq3}
\end{equation}
and its example is given in Figure \ref{fig2}.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure2.pdf}}
\caption{Max pooling takes the maximum value within the pooling block. The orange areas denote the locations from which the output values are taken.}
\label{fig2}
\end{figure}

\subsection{Average pooling}
\label{3B}
Average pooling returns the average of each pooling block in the feature map. The averaging function is given by
\begin{equation}
{o_{0,0}} = \frac{1}{4}({f_{0,0}} + {f_{0,1}} + {f_{1,0}} + {f_{1,1}}).
\label{eq4}
\end{equation}

In average pooling, all inputs of the corresponding block in the feature map equally contribute to the output. A simple example of average pooling is given in Figure \ref{fig3}.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure3.pdf}}
\caption{Average pooling averages the feature-map entries in each pooling block.}
\label{fig3}
\end{figure}

\subsection{Stride pooling (stride convolution)}
\label{3C}
Stride convolution reduces the size of the feature map by applying the stride in the convolution layer, without invoking an extra function. The convolution layer with a stride performs the same function as the convolution layer without a stride, followed by fixed-position pooling. The stride pooling is represented by
\begin{equation}
{o_{0,0}} = {f_{0,0}} ,
\label{eq5}
\end{equation}
and is demonstrated by a simple example in Figure \ref{fig4}.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure4.pdf}}
\caption{Stride convolution can be considered as stride pooling without the intermediate operations. Stride pooling takes the values from certain regions of the input.}
\label{fig4}
\end{figure}

\section{Universal Pooling}
\label{4}
This section proposes our new pooling method named universal pooling, which automatically determines the pooling mechanism by training the pooling weights. Universal pooling can be considered as a channel-wise local spatial attention module. Here we show that the proposed method encompasses max pooling, average pooling and stride pooling. For simplicity, we assume a pooling size of 2 and a stride of 2, but the concept is easily extended to more general cases.

\subsection{Motivation}
\label{4A}
The popular standard pooling methods, demonstrated by simple examples in Section III, are represented by linear combinations of the weighted features within a pooling block in the feature map. More specifically, the simple poolings are expressed as

\begin{align}
{o_{0,0}} & = {\mathrm{universal\_pooling}} \left( {f_{0,0}},{f_{0,1}},{f_{1,0}},{f_{1,1}} \right) \nonumber \\
 & = {\pi _{0,0}}{f_{0,0}} + {\pi _{0,1}}{f_{0,1}} + {\pi _{1,0}}{f_{1,0}} + {\pi _{1,1}}{f_{1,1}} \label{eq6}\\
 & \buildrel \Delta \over = {\boldsymbol{\pi }} \otimes {\boldsymbol{f}} \nonumber
\end{align}
where ${\boldsymbol{\pi }}$ denote the pooling weights, ${\boldsymbol{f}}$ is the input feature map, and the operator $\otimes$ performs linear combinations within the pooling block. Equation \eqref{eq6} interprets the pooling operation as an element-wise multiplication and summation within each pooling-size block in the feature map, as shown in Figure \ref{fig5}. In universal pooling, the pooling weight ${\boldsymbol{\pi }}$ is not manually selected, but is automatically determined from the input image or feature map. The pooling weight is also trained by back-propagation, by a mechanism similar to an attention map.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure5.pdf}}
\caption{Standard pooling can be considered as a linear combination within each pooling block. The weights in this linear combination are called the pooling weights.}
\label{fig5}
\end{figure}

\subsection{Structure}
\label{4B}
The pooling operation computes a scalar value representing the pooling block in the feature map. In implementations, the universal pooling operation should be divided into two smaller operations. The first sub-operation determines the contribution of each element in each pooling block of the feature map, which is equivalent to determining the pooling weight ${\boldsymbol{\pi }}$. The second sub-operation applies the $\otimes$ operator to the pooling weights and feature map, and returns the representative value of an element in the next-layer feature map.

The characteristics of the proposed universal pooling are similar to those of the existing pooling methods, and are subjected to the following constraints: (C1) Pooling is conducted in a channel-wise manner and excludes the information from all other channels. The pooling weight ${\boldsymbol{\pi }}$ is thus trained separately from channel to channel. (C2) Only the feature map elements in a given pooling block are pooled; the elements outside of the block are excluded. When training the pooling weights, the receptive fields of the weights are thus restricted to their corresponding pooling block. (C3) The scale of the features is unchanged by the pooling. For this purpose, the pooling weights are normalized such that ${\pi _{0,0}} + {\pi _{0,1}} + {\pi _{1,0}} + {\pi _{1,1}} = 1$. Our universal pooling is designed to satisfy constraints C1, C2, and C3. As shown in Figure 6, the universal pooling module has two main blocks: Block B1, in which the pooling weights are computed, and Block B2, in which the pooling is performed using the pooling weights.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure6.pdf}}
\caption{Structure of the proposed universal pooling module. The universal pooling module has two main blocks. B1 computes the pooling weight, and B2 performs the pooling using the pooling weights.}
\label{fig6}
\end{figure}

The pooling-weight block consists of several layers-the convolution, fully connected, batch normalization, and activation layers-as in other CNN modules. Using the ${\mathrm{softmax}} \left( \cdot \right)$ function, the output of this block is then normalized over the corresponding pooling block to obtain ${\pi _{0,0}} + {\pi _{0,1}} + {\pi _{1,0}} + {\pi _{1,1}} = 1$. The second block applies the $\otimes$ operator to the pooling weights and feature map over the pooling block to obtain the pooling result. To simplify the description, we assume identical size and stride number in the pooling method of the proposed approach, but unequal size and stride number are equally valid.

To clarify the relationship between the input and output of universal pooling, we give the equations of the proposed pooling. Before normalization, the output of the first block (B1) is represented as
\begin{equation}
\begin{array}{l}
{{{\boldsymbol{\overline{f}}}}^{\ell ,c}} = \phi ({{\boldsymbol{f}}^{\ell ,c}};{\boldsymbol{W}})\\
{{\boldsymbol{f}}^{\ell ,c}} = \left( {f_{i,j}^{\ell ,c}} \right) \in {\Re ^{H \times W}},{{{\boldsymbol{\overline{f}}}}^{\ell ,c}} = \left( {\overline{f}_{i,j}^{\ell ,c}} \right) \in {\Re ^{H \times W}}
\end{array}
\label{eq7}
\end{equation}
where $\phi ( \cdot )$ denotes the CNN module receiving ${\boldsymbol{f}}^{\ell,c}$, the feature map of channel $c$ in the $l$-th layer, and returning the temporary feature map ${{\boldsymbol{\overline{f}}}^{\ell ,c}} = \phi ({\boldsymbol{f}}^{\ell ,c};{\boldsymbol{W}})$ of channel $c$ (which is the same size as ${\boldsymbol{f}}^{\ell,c}$). ${\boldsymbol{W}}$ is a trainable parameter inside the network $\phi ( \cdot )$. To satisfy constraint (C3), the corresponding pooling block is processed by the softmax function, which guarantees that the pooling weights within the pooling block sum to one, as shown in Figure \ref{fig7}. This operation is represented by
\begin{equation}
\begin{gathered}
  {{\boldsymbol{\pi }}^{\ell ,c}} = {\mathrm{softmax}}({{{\boldsymbol{\overline{f}}}}^{\ell ,c}}) \hfill \\
  \pi _{pS + m,qS + n}^{\ell ,c} \hfill \\
  \quad = \frac{1}{{\sum\limits_{i = 0}^{S - 1} {\sum\limits_{j = 0}^{S - 1} {\exp (\overline{f}_{pS + i,qS + j}^{\ell ,c})} } }}\exp (\overline{f}_{pS + m,qS + n}^{\ell ,c}) \hfill \\
  {{\boldsymbol{\pi }}^{\ell ,c}} = \left( {\pi _{i,j}^{\ell ,c}} \right) \in {\Re ^{H \times W}} \hfill\\
  p = 0,1, \cdots ,\left\lfloor {{\raise0.7ex\hbox{$H$} \!\mathord{\left/
 {\vphantom {H S}}\right.\kern-\nulldelimiterspace}
\!\lower0.7ex\hbox{$S$}}} \right\rfloor  - 1,q = 0,1, \cdots ,\left\lfloor {{\raise0.7ex\hbox{$W$} \!\mathord{\left/
 {\vphantom {W S}}\right.\kern-\nulldelimiterspace}
\!\lower0.7ex\hbox{$S$}}} \right\rfloor  - 1 \hfill \\
\end{gathered}
\label{eq8}
\end{equation}
where $S$ denotes the stride number and size of the pooling, ${\mathrm{softmax}} \left( \cdot \right)$ is applied to an $S \times S$ pooling block, and  denotes the learnable pooling weight. Here, ${{\boldsymbol{\pi }}^{\ell ,c}}$ and $\overline{f}_{i,j}$ are the same size.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure7.pdf}}
\caption{Application of softmax within the pooling block. Each red square delineates a pooling block. Note that the pooling weights within each pooling block sum to one.}
\label{fig7}
\end{figure}

Then, the output of the first block (B1) is related to its input as follows:
\begin{equation}
{{\boldsymbol{\pi }}^{\ell ,c}} = {\mathrm{softmax}}(\phi ({{\boldsymbol{f}}^{\ell ,c}};{\boldsymbol{W}})).
\label{eq9}
\end{equation}

In \eqref{eq9}, the learnable parameter $\boldsymbol{W}$ inside the module determines the pooling function. To satisfy constraint (C1), $\phi ( \cdot )$ is applied only to the corresponding channel. The universal pooling then learns the different pooling functions in the different channels.
In the second block (B2), pooling weights are assigned to the feature map, and pooling is performed through the $ \otimes $ operation. The output from (B2) is represented by
\begin{equation}
\begin{gathered}
  {{\boldsymbol{o}}^{\ell ,c}} = {{\boldsymbol{\pi }}^{\ell ,c}} \otimes {{\boldsymbol{f}}^{\ell ,c}} \hfill \\
  o_{p,q}^{\ell ,c} = \sum\limits_{m = 0}^{S - 1} {\sum\limits_{n = 0}^{S - 1} {\pi _{pS + m,qS + n}^{\ell ,c} \cdot f_{pS + m,qS + n}^{\ell ,c}} }  \hfill \\
  {{\boldsymbol{o}}^{\ell ,c}} = \left( {o_{p,q}^{\ell ,c}} \right) \in {\Re ^{\scriptsize{\left\lfloor {{\raise0.7ex\hbox{$H$} \!\mathord{\left/
 {\vphantom {H S}}\right.\kern-\nulldelimiterspace}
\!\lower0.7ex\hbox{$S$}}} \right\rfloor  \times \left\lfloor {{\raise0.7ex\hbox{$W$} \!\mathord{\left/
 {\vphantom {W S}}\right.\kern-\nulldelimiterspace}
\!\lower0.7ex\hbox{$S$}}} \right\rfloor }}} \hfill \\
  p = 0,1, \cdots ,\left\lfloor {{\raise0.7ex\hbox{$H$} \!\mathord{\left/
 {\vphantom {H S}}\right.\kern-\nulldelimiterspace}
\!\lower0.7ex\hbox{$S$}}} \right\rfloor  - 1,{\text{ }}q = 0,1, \cdots ,\left\lfloor {{\raise0.7ex\hbox{$W$} \!\mathord{\left/
 {\vphantom {W S}}\right.\kern-\nulldelimiterspace}
\!\lower0.7ex\hbox{$S$}}} \right\rfloor  - 1{\text{ }} \hfill \\
\end{gathered}
\label{eq10}
\end{equation}
where the output ${{\boldsymbol{o}}^{\ell ,c}}$ from (B2) defines a new feature map in the $\left(\ell+1\right)$-th layer. Obviously, the scale of the features is unchanged and constraint (C2) is satisfied because the pooling weights are already normalized by \eqref{eq9}. In \eqref{eq10}, the pooling weights   can be rewritten as
\begin{equation}
{{\boldsymbol{o}}^{\ell ,c}} = {{\boldsymbol{\pi }}^{\ell ,c}} \otimes {{\boldsymbol{f}}^{\ell ,c}} = {\mathrm{softmax}}(\phi ({{\boldsymbol{f}}^{\ell ,c}};{\boldsymbol{W}})) \otimes {{\boldsymbol{f}}^{\ell ,c}}.
\label{eq11}
\end{equation}
Because all functions in this equation are simple multiplications and additions, all the learnable parameters are easily updated by back-propagation.

\subsection{Relationship between Universal Pooling and the Existing Pooling Methods}
\label{4C}
As mentioned above, the proposed pooling module can be considered as a channel-wise local spatial attention module. This subsection demonstrates that the popular existing pooling methods are special cases of universal pooling.

\subsubsection{Average Pooling}
\label{4CA}
Average pooling averages the elements inside the pooling block of the feature map. The averaging is performed as follows:
\begin{equation}
o_{p,q}^{\ell ,c} = \frac{1}{{{S^2}}}\sum\limits_{m = 0}^{S - 1} {\sum\limits_{n = 0}^{S - 1} {f_{pS + m,qS + n}^{\ell ,c}} }
\label{eq12}
\end{equation}
where the pooling block is an $S \times S$ matrix. Using the pooling weights ${{\boldsymbol{\pi }}^{\ell ,c}}$, we can rewrite \eqref{eq12} as
\begin{equation}
\begin{gathered}
  o_{p,q}^{\ell ,c} = \sum\limits_{m = 0}^{S - 1} {\sum\limits_{n = 0}^{S - 1} {\pi _{pS + m,qS + n}^{\ell ,c} \cdot f_{pS + m,qS + n}^{\ell ,c}} } , \hfill \\
  \pi _{pS + m,qS + n}^{\ell ,c} = \frac{1}{{{S^2}}} . \hfill \\
\end{gathered}
\label{eq13}
\end{equation}
The proposed universal pooling behaves like average pooling if all pooling weights are equal and sum to unity. The pooling weights are constant for any input image. This implies that before ${\mathrm{softmax}}( \cdot )$ the function, the output $\overline{f}^{\ell ,c}$ of (B1.1) is the same over the whole pooling block. Thus, when universal pooling performs like average pooling, its form reduces to
\begin{equation}
\begin{gathered}
  o_{p,q}^{\ell ,c} = \sum\limits_{m = 0}^{S - 1} {\sum\limits_{n = 0}^{S - 1} {\pi _{pS + m,qS + n}^{\ell ,c} \cdot f_{pS + m,qS + n}^{\ell ,c}} } , \hfill \\
  \pi _{pS + m,qS + n}^{\ell ,c} \hfill\\
  \quad = \frac{1}{{\sum\limits_{i = 0}^{S - 1} {\sum\limits_{j = 0}^{S - 1} {\exp (\overline{f}_{pS + i,qS + j}^{\ell ,c})} } }}\exp (\overline{f}_{pS + m,qS + n}^{\ell ,c}) \hfill\\
  \quad = \frac{1}{{{S^2}}} , \hfill \\
  \overline{f}_{pS + i,qS + j}^{\ell ,c} = \overline{f}_{pS + x,qS + y}^{\ell ,c}{\text{        }}(x \ne i,y \ne j) . \hfill \\
\end{gathered}
\label{eq14}
\end{equation}
In the simplest example of average pooling by the proposed universal pooling, the learnable parameter $\boldsymbol{W}$ is set to 0. In this case, all outputs are zero and the ${\mathrm{softmax}}( \cdot )$ function generates equal pooling weights that sum to one.

\subsubsection{Max Pooling}
\label{4CB}
Max pooling outputs the maximum value of the elements inside a given pooling block. The max pooling function can be written as follows.
\begin{equation}
o_{p,q}^{\ell ,c} = \mathop {\max }\limits_{0 \leqslant m,n \leqslant S - 1} f_{pS + m,qS + n}^{\ell ,c} .
\label{eq15}
\end{equation}

In terms of the pooling weights ${{\boldsymbol{\pi }}^{\ell ,c}}$, \eqref{eq15} can be rewritten as
\begin{equation}
\begin{gathered}
  o_{p,q}^{\ell ,c} = \sum\limits_{m = 0}^{S - 1} {\sum\limits_{n = 0}^{S - 1} {\pi _{pS + m,qS + n}^{\ell ,c} \cdot f_{pS + m,qS + n}^{\ell ,c}} } , \hfill \\
  \pi _{pS + m,qS + n}^{\ell ,c} = \left\{ {\begin{array}{*{20}{c}}
  1&{m,n = \mathop {\arg \max }\limits_{0 \leqslant i,j \leqslant S - 1} f_{pS + i,qS + j}^{l,c}} \\
  0&{otherwise}
\end{array}} . \right. \hfill \\
\end{gathered}
\label{eq16}
\end{equation}

The pooling weight is set to one at the site containing the maximum value in the block, and zero at all other sites. If the output $\overline{f}^{\ell ,c}$ of (B1.1) equals the input feature map $f^{\ell ,c}$ before invoking the ${\mathrm{softmax}}( \cdot )$ function, the universal pooling becomes
\begin{equation}
\begin{gathered}
  o_{p,q}^{\ell ,c} = \sum\limits_{m = 0}^{S - 1} {\sum\limits_{n = 0}^{S - 1} {\pi _{pS + m,qS + n}^{\ell ,c} \cdot f_{pS + m,qS + n}^{\ell ,c}} } , \hfill \\
  \pi _{pS + m,qS + n}^{\ell ,c} = {\mathrm{softmax}}({{{\boldsymbol{\overline{f}}}}^{\ell ,c}}) = {\mathrm{softmax}}({{\boldsymbol{f}}^{\ell ,c}}) \hfill \\
  \qquad \qquad \; \approx \left\{ {\begin{array}{*{20}{c}}
  1&{m,n = \mathop {\arg \max }\limits_{0 \leqslant i,j \leqslant S - 1} f_{pS + i,qS + j}^{\ell,c}} \\
  0&{otherwise}
\end{array}} .\right. \hfill \\
\end{gathered}
\label{eq17}
\end{equation}

In this case, the proposed universal pooling becomes very similar to max pooling. In summary, the condition under which the proposed method behaves like max pooling is
\begin{equation}
f^{\ell ,c} = \overline{f}^{\ell ,c} = \phi (f^{\ell ,c}) .
\label{eq18}
\end{equation}
That is, universal pooling reduces to max pooling when (B1.1) is trained as an identity mapping.

\subsubsection{Stride Pooling}
\label{4CC}
Stride pooling obtains an output from a specific location inside an area. When the specific location is the top-left site of the pooling block, the stride pooling is described by
\begin{equation}
o_{p,q}^{\ell ,c} = f_{pS,qS}^{\ell ,c}.
\label{eq19}
\end{equation}

In terms of the pooling weights, \eqref{eq19} is rewritten as
\begin{equation}
\begin{gathered}
  o_{p,q}^{\ell ,c} = \sum\limits_{m = 0}^{S - 1} {\sum\limits_{n = 0}^{S - 1} {\pi _{pS + m,qS + n}^{\ell ,c} \cdot f_{pS + m,qS + n}^{\ell ,c}} } , \hfill \\
  \pi _{pS + m,qS + n}^{\ell ,c} = \left\{ {\begin{array}{*{20}{c}}
  1&{m,n = 0,0} \\
  0&{otherwise}
\end{array}} . \right. \hfill \\
\end{gathered}
\label{eq20}
\end{equation}

The proposed universal pooling behaves like stride pooling when the first block (B1) outputs unity at the top-left location of the pooling block, and zero at the remaining positions.

Above, we showed that the proposed universal pooling includes the popular existing pooling methods as special cases. This implies that the performance of the existing pooling methods can be improved by replacing them with universal pooling. In particular, by appropriately training $\boldsymbol{W}$ in universal pooling, we can expect to improve the pooling to an unprecedented level, with concomitant improvement in CNN performance.

\section{Experimentation}
\label{5}
In this section, the proposed method is experimentally tested on two benchmark classification datasets: CIFAR10 \cite{b17} and Places2 \cite{b18}. The VGG19 and ResNet18 datasets are employed as backbone networks.

\subsection{Implementation}
\label{5A}
Pooling can be divided into local pooling and global pooling. In local pooling, the pooling block is smaller than the feature map, whereas in global pooling, it covers the entire feature map. In the present experiment, local pooling was implemented by one or two fully connected layers, as shown in Figure \ref{fig8a}. Global pooling also can be implemented either by fully connected layers (Figure \ref{fig8b}), or by several convolution layers (Figure \ref{fig8c}).

\begin{figure*}[!t]
    \centering

    \subfloat[]{
    \raisebox{11mm}{
        \includegraphics[width=0.25\textwidth]{figures/figure8-a.pdf}
        \label{fig8a}
    }}

    \hfill
    \subfloat[]{
    \raisebox{11mm}{
        \includegraphics[width=0.25\textwidth]{figures/figure8-b.pdf}
        \label{fig8b}
    }}

    \hfill
    \subfloat[]{
        \includegraphics[width=0.28\textwidth]{figures/figure8-c.pdf}
        \label{fig8c}
    }
    \caption{Local and global pooling implemented by fully connected and convolutional layers. The red squares delineate the pooling blocks. Local pooling shares the weights in the fully connected layers of the same channel. (a) Local pooling with a pooling module of fully connected layers; (b) Global pooling with a pooling module of fully connected layers; (c) Global pooling with a pooling module of several convolution layers.}
    \label{fig8}
\end{figure*}

\subsection{CIFAR10}
\label{5B}
The Canadian Institute for Advanced Research 10 (CIFAR10) dataset is a classification dataset containing 60,000 images \cite{b17}, each of size $32 \times 32$ on the RGB channels. The dataset contains ten object classes: airplane, automobile, bird, cat, deer, dog, frog, horse, ship, and truck. Typical images contained in the dataset are shown in Figure \ref{fig9}. The CNNs were trained by the stochastic gradient descent (SGD) method. The initial learning rate, momentum, and weight decay in the SGD were set to 0.1, 0.9, and 0.0001, respectively. The whole training procedure was implemented over 450 epochs, reducing the learning rate by one tenth at 150-epoch intervals.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure9.pdf}}
\caption{Examples from the CIFAR10 dataset, which comprises small $32 \times 32$ RGB images.}
\label{fig9}
\end{figure}

The backbone networks were two pre-trained models: VGG19 and ResNet18. The original VGG19 network consists of 16 convolution layers, 3 fully connected layers, and 5 max pooling layers. Instead of global pooling, VGG vectorizes the last feature map and passes it to the fully connected layer. The existing VGG structure, which is designed for large-sized images, was here modified for processing small images by employing one fully connected layer. As the last feature map is already vectorized by the final pooling process, it requires no further vectorization, so the final pooling can be considered as a global pooling. This structure of the network is given in Table \ref{tab1}.

The original VGG network (here termed V1) performs five max poolings. When the input image is a $32 \times 32$ CIFAR 10 image, the last-pooled map and the input feature map are the same size. Thus, the final pooling is considered as a global pooling while the preceding poolings are local poolings. As competing methods for comparison, we also employed five variants of VGG19. In variant V2, all five max poolings were replaced with their average pooling procedures. In variant V3, the local max pooling was replaced with local stride pooling in the first four pooling layers, and GAP was applied in the last pooling layer. The structure of V3 matched that of ResNet. Variants V4-V6 were constructed by various combinations of max and average pooling \cite{b15}, as shown in Table \ref{tab3}. In this table, "Mixed" refers to a linear combination of max and average pooling with real-numbered weights. Meanwhile, "Gated-channel" and "Gated-layer" refer to gated max-average pooling in which the gating mask is a linear combination of the weighted channels and weighted layers, respectively. Four variants of the proposed method, labeled P1-P5, were also constructed for a comparison evaluation. In variants P1-P3, both local and global poolings were implemented by fully connected layers, while in variants P4 and P5, the first four local poolings were implemented by fully connected layers, and the final (global) pooling was implemented by a convolution layer (see Table \ref{tab3}). Interestingly, the average pooling V2 (rather than the original VGG V1) delivered the best performance among the standard pooling methods. The combined max and average pooling methods (V4-V6) outperformed the standard pooling methods (V1-V3). Meanwhile, among variants P1-P5 of the proposed universal pooling method, constructed as competitors against the standard methods, P4 and P3 delivered the best and second-best performances, respectively. Overall, increasing the number of layers improved the network performance. However, the performance was degraded in the deepest network (P5), possibly because the network overfitted the training data.

Second, let us consider ResNet. The original ResNet18 consists of 17 convolution layers, 1 fully connected layer, 5 local pooling layers, and 1 global pooling layer. In the original ResNet, the local and global poolings are performed by striding and averaging (V3), respectively. The existing ResNet structure, which is suitable for large-sized images, was here modified for small-image processing by employing three (rather than five) local pooling layers. The network structure is given in Table \ref{tab2}. As in VGG, two standard pooling methods (V1 with max pooling and V2 with average pooling) and three combinations of max and average poolings (V4-V6) were constructed as competing methods. All experiments were conducted ten times, and the performances are summarized in Table \ref{tab2} and Figure \ref{fig10}. As confirmed in the table, average pooling delivered the best performance among the standard pooling methods (V1-V3, where V3 is the original ResNet18). The performance of average pooling (V2) was significantly higher than that of max pooling (V1). Interestingly, the average pooling outperformed even the combination methods (V4-V6). This suggests that training all poolings to behave like average pooling is a difficult task. However, the proposed universal pooling outperformed the average pooling, as shown in Table \ref{tab2}. The overall best performance was achieved by a network with one fully connected layer for local pooling, and two fully connected layers for global pooling. As demonstrated in this experiment, the optimal performance is dictated by the structure of the entire network, and the proposed universal pooling outperformed all standard pooling methods and their combinations proposed in \cite{b15}.

\begin{table}
\centering
\caption{VGG network, modified to suit the CIFAR10 dataset.}
\setlength{\tabcolsep}{3pt}
\begin{tabular}{ccc}
\toprule
Layer Name      & Output Size & Size, Channel, Stride                          \\
\midrule
Convolution 1   & $32 \times 32$ & $\left[3 \times 3, 64, 1\right] \times 2$  \\ \hline
Local Pooling1  & $16 \times 16$ & $\left[2 \times 2, 64, 2\right]$          \\ \hline
Convolution 2   & $16 \times 16$ & $\left[3 \times 3, 128, 1\right] \times 2$ \\ \hline
Local Pooling2  & $8 \times 8$   & $\left[2 \times 2, 128, 2\right]$         \\ \hline
Convolution 3   & $8 \times 8$   & $\left[3 \times 3, 256, 1\right] \times 4$ \\ \hline
Local Pooling3  & $4 \times 4$   & $\left[2 \times 2, 256, 2\right]$         \\ \hline
Convolution 4   & $4 \times 4$   & $\left[3 \times 3, 512, 1\right] \times 4$ \\ \hline
Local Pooling4  & $2 \times 2$   & $\left[2 \times 2, 512, 2\right]$         \\ \hline
Convolution 5   & $2 \times 2$   & $\left[3 \times 3, 512, 1\right] \times 4$ \\ \hline
Global Pooling  & $1 \times 1$   & $\left[2 \times 2, 512, 2\right]$         \\ \hline
Fully Connected & $1 \times 1$   & $\left[1 \times 1, 10, 1\right]$           \\
\bottomrule
\end{tabular}
\label{tab1}
\end{table}

\begin{table}
\centering
\caption{ResNet network, modified to suit the CIFAR10 dataset.}
\setlength{\tabcolsep}{3pt}
\begin{tabular}{ccc}
\toprule
Layer Name                 & Output Size    & Size, Channel, Stride                      \\
\midrule
Conv 1              & $32 \times 32$ & $\left[3 \times 3, 64, 1\right]$           \\ \hline
Conv 2-1 + Shortcut & $32 \times 32$ & $\left[3 \times 3, 64, 1\right] \times 2$  \\ \hline
Conv 2-2 + Shortcut & $32 \times 32$ & $\left[3 \times 3, 64, 1\right] \times 2$  \\ \hline
Conv 3-1 + Shortcut & $32 \times 32$ & $\left[3 \times 3, 128, 1\right] \times 2$ \\ \hline
Local Pooling 1            & $16 \times 16$ & $\left[2 \times 2, 128, 2\right]$          \\ \hline
Conv 3-2 + Shortcut & $16 \times 16$ & $\left[3 \times 3, 128, 1\right] \times 2$ \\ \hline
Conv 4-1 + Shortcut & $16 \times 16$ & $\left[3 \times 3, 256, 1\right] \times 2$ \\ \hline
Local Pooling 2            & $8 \times 8$   & $\left[2 \times 2, 256, 2\right]$          \\ \hline
Conv 4-2 + Shortcut & $8 \times 8$   & $\left[3 \times 3, 256, 1\right] \times 2$ \\ \hline
Conv 5-1 + Shortcut & $8 \times 8$   & $\left[3 \times 3, 512, 1\right] \times 2$ \\ \hline
Local Pooling 3            & $4 \times 4$   & $\left[2 \times 2, 512, 2\right]$          \\ \hline
Conv 5-2 + Shortcut & $4 \times 4$   & $\left[3 \times 3, 512, 1\right] \times 2$ \\ \hline
Global Pooling             & $1 \times 1$   & $\left[4 \times 4, 512, 4\right]$          \\ \hline
Fully Connected            & $1 \times 1$   & $\left[1 \times 1, 10, 1\right]$           \\
\bottomrule
\end{tabular}
\label{tab2}
\end{table}

\begin{table*}
\centering
\caption{Classification performances of the evaluated CNNs on the CIFAR10 dataset (averaged over 10 experiments). Best result bold, second-best underlined.}
\setlength{\tabcolsep}{3pt}
\begin{tabular}{cccccc}
\toprule
\multirow{2}{*}{Pooling method}                          & \multicolumn{2}{c}{Pooling Function} & \multirow{2}{*}{Indicator} & \multicolumn{2}{c}{Backbone Network} \\
                                                         & Local pooling    & Global pooling    &                            & VGG               & ResNet           \\
\midrule
\multirow{3}{*}{Standard pooling methods}                & Max              & Max               & V1                         & 93.218            & 94.335           \\
                                                         & Average          & Average           & V2                         & 93.352            & 95.140           \\
                                                         & Stride           & Average           & V3                         & 92.607            & 94.639           \\ \hline
\multirow{3}{*}{Conventional pooling methods \cite{b15}} & Mixed            & Mixed             & V4                         & 93.388            & 94.127           \\
                                                         & Gated-channel    & Gated-channel     & V5                         & 93.197            & 94.548           \\
                                                         & Gated-layer      & Gated-layer       & V6                         & 93.375            & 94.974           \\ \hline
\multirow{5}{*}{Module architecture in proposed methods} & 1 FC             & 1 FC              & P1                         & 93.255            & 95.192           \\
                                                         & 1 FC             & 2 FC              & P2                         & 93.354            & \textbf{95.231}  \\
                                                         & 2 FC             & 2 FC              & P3                         & \underline{93.405}      & \underline{95.201}     \\
                                                         & 1 FC             & Conv              & P4                         & \textbf{93.418}   & 95.199           \\
                                                         & 2 FC             & Conv              & P5                         & 93.260            & 95.163           \\
\bottomrule
\end{tabular}
\label{tab3}
\end{table*}

\begin{figure*}[!t]
    \centering

    \subfloat[]{
        \includegraphics[width=0.45\textwidth]{figures/figure10-a.pdf}
        \label{fig10a}
    }

    \quad\quad\quad
    \subfloat[]{
        \includegraphics[width=0.45\textwidth]{figures/figure10-b.pdf}
        \label{fig10b}
    }
    \caption{Boxplot of experiments on the CIFAR10 dataset, performed in the (a)VGG architecture, and (b)the ResNet architecture. Red lines denote the experimental averages. The proposed method outperformed the other pooling methods.}
    \label{fig10}
\end{figure*}

\subsection{Places2}
\label{5C}

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure11.pdf}}
\caption{Examples from the Places2 dataset, which contains images of various places.}
\label{fig11}
\end{figure}

The Places2 dataset is a classification dataset containing 365 place classes. The dataset includes 1,803,460 training images and 50 validation images per class. Example images from the Places2 dataset are given in Figure 11. The network was trained by the SGD optimizer and the entire learning period was 120 epochs with a batch size of 256. The initial learning rate, momentum and weight decay of the SGD were set to 0.1, 0.9, and 0.0001, respectively. The learning rate was reduced by one tenth at 30-epoch intervals. The Places2 experiment was conducted in the pre-trained ResNet18 model, and the results are given in Table 2. Clearly, the classification accuracy was improved by the universal pooling.

\begin{table*}
\centering
\caption{Classification performance in the Places2 experiment.}
\setlength{\tabcolsep}{3pt}
\begin{tabular}{ccccc}
\toprule
\multirow{2}{*}{Pooling method} & \multicolumn{2}{c}{Pooling Function} & \multicolumn{2}{c}{ResNet18}      \\
                                & Local pooling    & Global pooling    & top1            & top5            \\
\midrule
Standard pooling methods        & Stride           & Average           & 53.963          & 83.778          \\ \hline
Proposed pooling methods        & 2 FC             & Conv              & \textbf{54.693} & \textbf{84.781} \\
\bottomrule
\end{tabular}
\label{tab4}
\end{table*}

The universal pooling function was qualitatively evaluated by varying the global pooling weights in the universal pooling method. The results of training with the different weights are visualized in Figures \ref{fig12}-\ref{fig15}. As implied in the three cases below, the universal pooling can train various pooling functions, and appropriately selects the pooling function to enhance the classification performance.

\subsubsection{Training by Average Pooling}
\label{5CA}
In almost one-third of the channels (129 out of 512 channels), the pooling weights are the same over the whole feature map, and the pooling reduces to average pooling. The first of these 129 channels is depicted in Figure \ref{fig12}. The trained pooling weights exhibited the same values over the input feature map.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure12.pdf}}
\caption{Pooling weights trained by average pooling (top) and the corresponding feature maps (bottom). The trained pooling weights have the same values over the input feature map.}
\label{fig12}
\end{figure}

\subsubsection{Training by Flexible Pooling}
\label{5CB}
In almost a fifth of the channels (95 out of 512 channels), the pooling weights depended on the input (see Figure \ref{fig13}). This pooling, called flexible pooling, includes max pooling as a special case. Figure \ref{fig13} shows how the pooling weights change with the feature maps. The trained universal pooling operates similarly, but not identically, to max pooling. In certain data, such as data \#1 and \#4 in the figure, the trained pooling exhibited similar characteristics to fixed pooling (discussed next).

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure13.pdf}}
\caption{Pooling weights trained by flexible pooling (top) and the corresponding feature maps (bottom). The pooling weights depend on the input feature map. The pooling weights resemble max pooling weights except those of data \#1 and data \#4, which are reminiscent of fixed pooling.}
\label{fig13}
\end{figure}

\subsubsection{Training by Fixed Pooling}
\label{5CC}
In more than one-half of the 512 channels (288 out of 512 channels), the trained pooling weights were almost fixed, and were nonresponsive to changes in the input feature map. This pooling, called fixed pooling, includes both average and stride pooling. However, average pooling is classified separately from other fixed poolings because the conditions of average and stride pooling differ in the universal pooling. The results of two channels with fixed training weights are visualized in Figures \ref{fig14} and \ref{fig15}. In both cases, the pooling weights depended only on the pooling position in the image, and were independent of the input feature map.

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure14.pdf}}
\caption{Pooling weights trained by fixed pooling in the proposed method (top) and the corresponding feature maps (bottom). In this channel, the pooling takes the features from the center of the image.}
\label{fig14}
\end{figure}

\begin{figure}[htbp]
\centerline{\includegraphics[width=0.9\linewidth]{figures/figure15.pdf}}
\caption{As in Figure \ref{fig14}, but with the pooling features taken from the surrounding area of the image.}
\label{fig15}
\end{figure}

\section{Conclusion}
\label{6}
This paper proposed a new pooling method named universal pooling, which dynamically learns the pooling weights from the given dataset, and which includes the existing popular pooling methods as special cases. Thus, replacing the existing pooling with universal pooling will likely improve the network performance. The proposed pooling was applied to two benchmark datasets and two pretrained models and was validated in comparison experiments.

{\small
\bibliographystyle{ieee_fullname}


\section{Introduction}

Quantum Field Theories (QFT) are celebrated for being the framework of the Standard Model and for making predictions which coincide with experimental data to extreme accuracy. Yet QFT is perceived as difficult to learn by many and this seems to be at least partially due to the lack of visualization options which would help to develop an intuition for the theory. In  quantum mechanics (or `first quantization'), for example, a localized particle is associated with a $\delta$-function and a uniformly moving particle with a plane wave, both of which are relatively easy to depict in a graph. Many textbooks show the energy eigenstates of hydrogen atoms, harmonic oscillators or quantum wells as plotted wave functions. A survey of visualization methods in quantum mechanics can be found in Thaller's books \cite{Thaller}$^,$\cite{Thaller2}. In QFT, on the other hand, visualizations of the fundamental concepts seem to be scarce. One prominent exception are the Feynman diagrams, which are an important tool in studying processes of interacting particles. Yet these diagrams are firmly connected to the specific computational method of perturbation theory and they are not meant as a graphical representation of the quantum field itself. Another approach is to compute and plot statistical quantities, like the particle density or the two point correlation function of a many body system. But a lot of information is lost when moving from the full description of its state to these aggregated quantities. Therefore they tend to be more useful to understand macroscopic properties of the system rather than the underlying details.

The difficulties in visualizing a quantum field obviously stem from the fact that such a system has too many degrees of freedom: Even if the field is simplified to a discrete lattice with $N$ atoms oscillating in only one space dimension each (like, e.g., in the `quantized mattress' model in Zee's book\cite{Zee}), the quantum state is given by a wave function  $\Psi: \mathbb{R}^N \rightarrow \mathbb{C}$, which cannot be plotted in a straight-forward way for any $N > 2$.

Such lattices of coupled quantum oscillators are also studied in several other areas of physics and especially in solid state physics they are applied to model phonons in a crystal.
Johnson and Gutierrez\cite{JohnsonGutierrez} visualize phonon states of a one-dimensional quantum lattice by projecting the probability density of the system to each of the one-dimensional spaces in which the atoms oscillate.

In the present article we propose a new way of visualizing the state of a bosonic quantum field and apply it to the example of the harmonic chain, which is essentially a discretized, one-dimensional model for a boson field. The target audience for this material are mainly graduate students who completed the courses on quantum mechanics and who are now about to advance towards quantum field theory. We assume that the reader is familiar especially with the treatment of the quantum harmonic oscillator and base changes in Hilbert space via Fourier transformation.

The article is organized as follows: In Section \ref{SecMonteCarlo} we present a Monte Carlo method for visualizing quantum wave functions in (first quantization) quantum mechanics. Then we introduce a quantum harmonic chain in Section \ref{SecLinHarChain} and we compute its dynamics by diagonalizing its Hamiltonian via a Fourier transformation. In Section \ref{ConnectingDots} we explain how the harmonic chain can be interpreted in the pictures of solid state physics (where each atom in the chain is a particle) and in quantum field theory (where a particle is an excitation of the chain). Setion \ref{SecVisualizationOfChain} contains our main result - a new visualization method for the state of a quantum harmonic chain and the application to several interesting special cases. We conclude with a discussion of our findings and an outlook to future research in Section \ref{SecDiscussion}.

\section{Monte Carlo plot of wave functions} \label{SecMonteCarlo}

De Aquino et. al. have proposed the following method to visualize the wave function $\Psi: \mathbb{R}^2 \rightarrow \mathbb{C}$ of a single quantum particle in two dimensions\cite{DeAquinoEtAl}: Associate $\mathbb{R}^2$ with the two orthogonal axes of a scatter plot and then draw many dots into this plane, with the location of each dot being randomly chosen according to the probability density function $|\Psi|^2$. The resulting charts show a high density of dots where $|\Psi|^2$ is relatively large and only very few dots where $\Psi$ is close to zero.

In our visualization method we follow a similar Monte Carlo approach, but we impose the following rules:
\begin{enumerate}
\item Choose the locations of the dots randomly according to a uniform probability distribution in a rectangular window within $\mathbb{R}^2$.
\item Represent the value of $\Psi$ at each dot via its color. For energy eigenstates, which can always be written as real-valued functions, such a visualization is possible on a black and white scale. For complex-valued functions a color spectrum can be used.
\end{enumerate}

To start with a simple example from ordinary quantum mechanics, consider a two-dimensional harmonic quantum oscillator with the Hamiltonian
\begin{equation}
\hat H = \frac{\hat{p_1}^2}{2m}+\frac{\hat{p_2}^2}{2m} + \frac{1}{2} \kappa (\hat{q_1}^2 + \hat{q_2}^2),
\end{equation}
where $q_l$ with $l \in \{1;2\}$ are the position coordinates of the particle with mass $m$, $\hat{p_l} = -i \partial/\partial q_l$ are the momentum operators, and $\kappa$ controls the strength of the attractive potential. Here and in the rest of the article we use units of measure such that $\hbar = 1$. It is well known\cite{Liboff} that the Hamiltonian is separable and the (non-normalized) energy eigenstates of the system can be written as
\begin{equation}
\Psi_{\nu_1,\nu_2} = \prod_{l=1,2}{H_{\nu_l}(\sqrt{m\omega}\,q_l)e^{-m\omega q_l^2/2}},
\end{equation}
where $\omega^2 = \kappa/m$,  $H_\nu$ is the $\nu$-th order Hermite polynomial and $\nu_{1,2} \in \mathbb{N}_0$ are the two quantum numbers which enumerate the energy eigenstates of the system.
As an example, Fig. \ref{2DHarmonic} shows a Monte Carlo plot of the energy state $\Psi_{2,1}$ created according to the two rules above. \BW{Black (white)}\color{Blue (red)} dots represent negative (positive) values of $\Psi_{2,1}$ and dots where $\Psi_{2,1} \approx 0$ blend in with the chart's \BW{grey}\color{white} background\BW{, i.e., the brighter a dot is, the larger is the value of $\Psi_{2,1}$ in the respective point of $\mathbb{R}^2$}. We have omitted a quantitative color scale since the normalization of the wave function is not important. It is easy to recognize the general shape of $\Psi_{2,1}$ and its nodal lines.

Obviously such a Monte Carlo plot is only of limited relevance in two dimensions, since we could have simply colored the whole chart area according to the value of $\Psi_{2,1}$ instead of just picking a few points, or we could have drawn the wave function in a 3D plot. We will see later that such a Monte Carlo visualization actually becomes quite useful when moving beyond two dimensions.

\begin{figure}[h!]
\centering
\BW{\includegraphics[width=180mm]{2DHarmonic-BW.eps}}
\color{\includegraphics[width=180mm]{2DHarmonic-COL.eps}}
\caption{Monte Carlo representation of the two-dimensional harmonic oscillator energy eigenstate $\Psi_{2,1}$. \BW{Bright (dark)}\color{Red (blue)} dots correspond to positive (negative) values of the wave function. Horizontal and vertical nodal lines of the function can be recognized as areas without any visible dots. As expected, the state $\Psi_{2,1}$ has two nodal lines intersecting with the $q_1$ axis and one intersecting with the $q_2$ axis.}
\label{2DHarmonic}
\end{figure}

\section{Quantized Linear Harmonic Chain} \label{SecLinHarChain}

Consider a linear chain of $N$ coupled harmonic oscillators whose positions are defined by the coordinates $q_n$ with $n = 1,...,N$, respectively. In classical mechanics this system's Hamiltonian is
\begin{equation}\label{EqClassicalH}
H = \sum_{n=1}^N \left(\frac{p_n^2}{2m} + \frac{1}{2} \kappa q_n^2 + \frac{1}{2} \gamma (q_n - q_{n+1})^2\right),
\end{equation}
where $p_n = m \dot q_n$. $\kappa > 0$ determines the strength of binding between the oscillating masses and their respective neutral positions, while $\gamma > 0$ controls the coupling between neighbours. For reasons of simplicity we impose periodic boundary conditions, i.e. $q_0$ is the same as $q_N$ and $q_{N+1}$ means just $q_1$. This system is the one-dimensional version of the `mattress' that Zee\cite{Zee} uses to introduce quantum field theory, or it can be seen as a simple model for phonons traveling through a crystal lattice.

We shall now follow the standard procedure: We will decompose the system's dynamics into normal modes which can be excited and evolve in time independently of one another. Mathematically speaking, this means to diagonalize the Hamiltonian $H$ via a Fourier transformation, revealing the system's equivalence to $N$ decoupled harmonic oscillators. One option is to expand $\vec q$ in real-valued sine and cosine modes, like in the article of  Johnson and Gutierrez\cite{JohnsonGutierrez}. Alternatively, a complex-valued decomposition into modes of $e^{ikn}$ can be applied like in the book of Greiner\cite{Greiner} on p.18. We shall also follow the latter approach and we will exhibit the technical steps in detail in order to make this article as self-contained as possible. We will deviate slightly from Greiner's path since we also need to find (among other things) the explicit coordinate transformation between the harmonic chain and the decoupled oscillators (eq. \ref{EqQtildeq} below).

We write the Hamiltonian (\ref{EqClassicalH}) as
\begin{equation}
H = \frac{|\vec{p}|^2}{2m} + \frac{1}{2} {\vec{q}}^T D \vec{q},
\end{equation}
where $D$ is an $N \times N$ matrix taking care of the terms with $\kappa$ and $\gamma$ in (\ref{EqClassicalH}):
\begin{equation}
D =
\begin{pmatrix}
\kappa+2\gamma & -\gamma & 0 & 0 & \dots & 0 & -\gamma\\
-\gamma & \kappa+2\gamma & -\gamma & 0 & \dots &0 & 0\\
0 & -\gamma & \kappa+2\gamma  & -\gamma & \dots &0 & 0\\
\dots  & \dots & \dots & \dots & \dots &\dots & \dots
\end{pmatrix}
\end{equation}
The matrix $D$ is a circulant, which means that each row vector is rotated one element to the right compared to the preceding row vector. In this sense $D$ is generated from the $N$-vector $(\kappa + 2\gamma, -\gamma, 0, 0, ..., 0 , -\gamma)$ which forms the top line of the matrix. Now remember that any circulant matrix can be diagonalized via a discrete Fourier transform\cite{Davis} which, in turn, means to expand the position coordinates $q_n$ with respect to the orthonormal basis $\{\vec{f^{(k)}}\}$ with
\begin{equation} \label{EqONB}
f^{(k)}_n = \frac{1}{\sqrt N} e^{-\frac{2\pi i}{N}k n} \quad \textrm{for} \quad k = -\frac{(N-1)}{2}, \dots, \frac{(N-1)}{2}.
\end{equation}
The vectors  $\vec{f^{(k)}}$ are also known as the `normal modes' of the system. We have made the assumption that $N$ is an odd number in order to save us a cumbersome treatment of special cases. Since $N$ is considered to be a large number in all practically relevant cases, this is no real restriction.
We call the position coordinates in the new basis $Q_k$ and they are related to the old ones by
\begin{equation}\label{EqFourierDecomp}
q_n =  \sum_k Q_k f^{(k)}_n = \frac{1}{\sqrt N} \sum_k Q_k e^{-\frac{2\pi i}{N}k n}.
\end{equation}
Here and in the rest of the article a sum over $k$ always means that $k$ runs from $-(N-1)/2$ to $+(N-1)/2$ unless explicitly stated otherwise. Note that the $Q_k$ are complex numbers. In order to ensure that the coordinates $q_n$ remain real-valued, we have to impose the constraint
\begin{equation} \label{EqConstraint}
Q_k = \overline{Q_{-k}}.
\end{equation}
An analogous transformation maps the momentum $\vec p$ to the momentum $\vec P$ with $P_k = m \dot Q_k$. In the new coordinates the Hamiltonian takes the form
\begin{equation} \label{EqHnew}
H = \sum_k \left(\frac{|P_k|^2}{2m} + \frac{1}{2} \omega_k |Q_k|^2 \right).
\end{equation}
Here we have used the fact that our transformation was essentially a discrete Fourier transformation and it diagonalized the circulant matrix $D$. This implies that the eigenvectors of $D$ are the base vectors $\vec{f^{(k)}}$. We have chosen to call the respective eigenvalues $\omega_k$. Since $D$ is real symmetric, the $\omega_k$ are real. We have
\begin{equation}
\omega_k \vec{f^{(k)}} = D \vec{f^{(k)}} = D \overline{\vec{f^{(-k)}}} = \overline{D \vec{f^{(-k)}}} = \overline{\omega_{-k} \vec{f^{(-k)}}} = \omega_{-k} \vec{f^{(k)}}
\end{equation}
and thus $\omega_k=\omega_{-k}$. In order to get rid of the complex coordinates before doing a canonical quantization, we perform a second coordinate transformation, splitting up each $Q_k$ into two real-valued coordinates:
\begin{equation} \label{EqQtilde}
Q_k = \frac 12 (1+i) \tilde Q_k + \frac 12 (1-i) \tilde Q_{-k} \quad\textrm{with } \tilde Q_k \in \mathbb{R}.
\end{equation}
The constraint (\ref{EqConstraint}) is automatically observed then. By putting (\ref{EqQtilde}) into (\ref{EqFourierDecomp}) we get the transformation rule between the $q_n$ and the $\tilde Q_k$:
\begin{equation} \label{EqQtildeq}
q_n = \sum_k \tilde Q_k \frac{1}{2} \left((1+i) f^{(k)}_n+(1-i)f^{(-k)}_n\right) =: \sum_k \tilde Q_k \tilde f^{(k)}_n
\end{equation}
We note that the $\{\tilde f^{(k)}_n\}$, in which we have developed $\vec q$, form another orthonormal basis in $\mathbb{R}^N$. This can be verified by direct computation, using the fact that the $\{f^{(k)}_n\}$ are also orthonormal.

The system's Hamiltonian (\ref{EqHnew}) looks very similar in the new coordinates:\\
\begin{equation} \label{EqHnew2}
H = \sum_k \left(\frac{\tilde P_k^2}{2m} + \frac{1}{2} \omega_k \tilde Q_k^2\right),
\end{equation}
where $\tilde P_k = m  \partial/\partial_t \tilde Q_k$. $H$ is now separated into $N$ one-dimensional harmonic oscillators with the eigenfrequencies $\omega_k$. Next, we perform a canonical quantization on $H$, promoting the location and momentum coordinates to Hilbert space operators. As a result we get a set of $N$ decoupled quantum oscillators. Solving the harmonic oscillator is one of the most common introductory problems in quantum mechanics and it is presented in many textbooks\cite{Liboff}. It can be achieved either by solving the oscillator's differential (Schroedinger) equation, or in an algebraic way based on commutation relations. Following the latter approach, we introduce the usual creation/annihilation operators $a^\dagger_k$ and $a_k$, so that the quantum version of (\ref{EqHnew2}) becomes
\begin{equation} \label{EqHnew3}
\hat H = \sum_k \omega_k \left(a^\dagger_k a_k + \frac 12 \right).
\end{equation}
The excitations of the harmonic chain created by the $a_k^\dagger$ are also called `particles' since this is how they often manifest themselves in experiments.

\section{Connecting the dots} \label{ConnectingDots}

By now we have established all the machinery needed to advance to the main ideas of this article, the Monte Carlo visualization of bosonic quantum fields. But let's first step back  and review how our work so far is related to first and second quantization. This section is meant to remind ourselves of the `big picture' and provide some orientation to novices in the field. It can be omitted by the advanced reader.

We had started with one quantum particle which was free to move in two dimensions $q_1$ and $q_2$ subject to an harmonic potential. As known from introductory quantum mechanics courses, the particle's wave function (or, more precisely, the square of its absolute value) tells us the probability density of finding the particle at a given point in space when performing a measurement. We have visualized such a wave function for the specific eigenstate $\Psi_{2,1}$ in Fig. \ref{2DHarmonic}. The coordinates $q_1$ and $q_2$ are linked to location operators $\hat q_1$ and $\hat q_2$ which in turn describe a quantum measurement of the particle's position  in the respective coordinate.

Then we have shown how to solve the dynamics of the harmonic oscillator chain, which can serve for at least two physical models:

On the one hand, we can think of each point in the chain as being an actual physical particle like the atoms in a crystal lattice. In this case $q_n$ is still the spacial coordinate of each mass point relative to its equilibrium position, just like in the two dimensional case before, and $n$ simply enumerates the atoms of the chain. We are still in the `first quantization world' and if we were to measure the exact positions of the $N$ particles at the same time, the wave function of the system would tell us the probability of finding the chain is a certain configuration. The $\hat q_n$ remain simple location operators which describe a position measurement of the $n$th particle. For example, if the system's normalized state is $\Psi$ and we measure the position of the third particle, then the expectation value of the result is $\langle \Psi, \hat q_3 \Psi \rangle$. The normal modes of oscillation in the chain, which we have identified in Section \ref{SecLinHarChain}, are called phonons - an important concept in solid state physics.

On the other hand, the chain can serve as a model for a bosonic quantum field, which belongs to the realm of `second quantization'. In this case $n$ becomes the discrete version of a spacial coordinate and $q_n$ is the field strength at the point $n$. The wave function of the system would give us the probability to encounter the field in a given classical state if it were possible to measure the field strength at all points in space at the same time. The operators $\hat q_n$ are now called `field operators' and they correspond to a measurement of the field strength in the spacial point $n$. In quantum field theory one would typically proceed to the continuum limit, replacing $n$ by $\vec x \in \mathbb{R}^3$ and $\hat q_n$ by some $\hat\phi(\vec x)$, for example. So all we have is a field, but where are particles - the bosons - which we wanted to describe? It turns out that what we call `particles' are simply excitation modes of the quantum field. Such excitations can be limited to a region in space and they can travel through the field very similarly to how particles in first quantization can be more or less localized and move in space. But, in contrast to first quantization, excitations can also become stronger or weaker - which means that particles can be created or destroyed.

In the following section we shall see some examples of such excitations, how they can be visualized, and how they correspond to physical particles.

\section{Visualization of the Quantum Harmonic Chain} \label{SecVisualizationOfChain}

\subsection{Visualization Method}

We have noted in Section \ref{SecLinHarChain} that the Hamiltonian of the quantum harmonic chain of length $N$ is equivalent to the $N$-dimensional quantum harmonic oscillator (\ref{EqHnew2}) up to a coordinate transformation. The latter has the convenient property of being separable, so that its energy eigenstates are simply products of the familiar one-dimensional harmonic oscillator states:
\begin{equation}
\Psi_{\{\nu_k\}} = \prod_{k}{H_{\nu_k}(\sqrt{m\omega_k}\,\tilde Q_k)e^{-m\omega_k\tilde Q_k^2/2}},
\end{equation}

where $H_\nu$ is the $\nu$-th order Hermite polynomial and $\nu_k \in \mathbb{N}_0$ for $k = -(N-1)/2,\dots,(N-1)/2$ are the quantum numbers.
Any wave function $\Psi$ which is given in terms of the coordinates $\tilde Q_k$ can be transformed back to the original coordinates $q_k$ via (\ref{EqQtildeq}). Once $\Psi(\vec q)$ is known, we can apply the Monte Carlo visualization which we have introduced in section \ref{SecMonteCarlo}. Compared to our simple example in Fig. \ref{2DHarmonic}, the main difference is that we now have $N$ coordinates instead of only two. Therefore we have to adapt our visualization prescription as follows:
\begin{enumerate}
\item Choose points $\vec q$ randomly according to a uniform probability distribution in a rectangular window within $\mathbb{R}^N$.
\item Each point $\vec q$ is visualized as a polyline in a parallel axes plot and the value of $\Psi(\vec q)$ is represented by the color of the polyline.
\end{enumerate}
The visualization in a parallel axes diagram has the advantage that each point $\vec q$ is represented by a polyline which can be intuitively associated with the corresponding state of a classical (i.e., non-quantum) oscillator chain.

In order to plot such visualizations we have developed a small program in the numerically oriented programming language Scilab\cite{Scilab}. The source code is available on request from the author of this article.

\subsection{Ground State}

As a first example of the results, Fig. \ref{GroundState} represents the ground state $\Psi_0 = \Psi_{\{\nu_k=0\,\forall k\}}$ of a quantum chain consisting of $N = 15$ oscillators. In the language of QFT, this state is called the `vacuum' since no excitations (=particles) are present. The wave function of the ground state is unique only up to multiplication with a complex number and in this case it has been chosen such that $\Psi_0$ is real and positive.  Each line in the plot represents one possible configuration of the 15 oscillators' positions. The \BW{brightness}\color{color} of each line corresponds to the value of the wave function $\Psi_0$ for the respective configuration. The color scale is chosen such that the value $\Psi_0 = 0$ blends in with the \BW{grey}\color{white} background and the \BW{brighter}\color{more colorful} the line, the higher the value of $\Psi_0$. Each polyline is plotted on top of the lines with lower $|\Psi_0|$, so that the most important lines (i.e. those with high $|\Psi_0|$) are more clearly visible in the chart. Those lines represent the most probable configurations of the chain in the sense of a quantum mechanical measurement and they tend to be relatively close to the $\vec q = 0$ line, which is the classical equilibrium state of the system.

\begin{figure}[h!]
\centering
\resizebox{0.99\textwidth}{!}{
  \BW{\includegraphics{GroundState-BW-v2.eps}}
  \color{\includegraphics{GroundState-COL-v2.eps}}
}
\caption{Monte Carlo representation of the ground state of a quantum harmonic chain. The wave function has been chosen real-valued and positive, thus all the lines are \BW{brighter than the background}\color{of the same color}. The \BW{brightest} lines \color{in deepest red }correspond to those $\vec q$ where the wave function assumes the highest values. Not surprisingly, in the ground state these are the lines close to $\vec q = 0$. Actually, a straight horizontal line corresponding exactly to $\vec q = 0$ would be drawn in the \BW{brightest}\color{deepest} color since it represents the global maximum of the ground state wave function. But due to the random nature of the Monte Carlo approach this line happens not to be drawn in the chart.}
\label{GroundState}
\end{figure}

\subsection{Particles at rest}

\begin{figure}[h!]
\centering
\resizebox{0.99\textwidth}{!}{
  \BW{\includegraphics{OnePartRest-BW-v2.eps}}
  \color{\includegraphics{OnePartRest-COL-v2.eps}}
}
\caption{Monte Carlo representation of the first excited state of the harmonic chain with wave number zero. In the most prominent configurations of the chain, i.e., those corresponding to \BW{very bright}\color{deeply red} or \BW{very dark}\color{deeply blue} lines, the $q_n$'s tend to be either jointly positive or all negative. Thus, if we could perform a quantum measurement of the exact position of the harmonic chain, we would most likely either find it entirely shifted towards the positive $q$-direction or entirely the other way. In the particle interpretation of a quantum field this state would be called `one particle at rest'.}
\label{OnePartRest}
\end{figure}

\begin{figure}[h!]
\centering
\resizebox{0.99\textwidth}{!}{
  \BW{\includegraphics{TwoPartRest-BW-v2.eps}}
  \color{\includegraphics{TwoPartRest-COL-v2.eps}}
}
\caption{The second excited state of the harmonic chain with wave number zero is similar to the first one shown in Fig. \ref{OnePartRest}. But where the first excited state has only one `nodal line' at $\vec q = 0$ (or, much more precisely, one nodal hypersurface which contains the point $\vec q = 0$), we can now recognize two such `nodal lines' where the graph turns from being mostly \BW{bright}\color{red} to mostly \BW{dark}\color{blue} and back again.  In the particle interpretation of a quantum field, this state would be called `two particles at rest'.}
\label{TwoPartRest}
\end{figure}

As further examples we consider excitations of the chain which can be interpreted as one or two particles at rest, respectively. The one-particle state $\Psi = a_0^\dagger \Psi_0$ is visualized in Fig. \ref{OnePartRest} and the two-particle state $\Psi = (a_0^\dagger)^2 \Psi_0$ in Fig. \ref{TwoPartRest}. As before, the color scale is chosen such that configurations with $\Psi = 0$ blend in invisibly with the \BW{grey}\color{white} background, and the \BW{brightest/darkest}lines \color{in deepest red/blue }correspond to the largest (positive) / smallest (negative) values of $\Psi$. Similar to the ground state, the configurations with the  highest absolute value of $\Psi$ tend to be those with little fluctuation along the $n$-axis. But in the one-particle case we note that the polylines with mostly positive (negative) $q_n$ tend to correspond to positive (negative) values of $\Psi$. In the two-particle case the polylines close to the $\vec q = 0$ line tend to correspond to negative values of $\Psi$ while those farther away from that line are rather associated with positive values of $\Psi$. This pattern obviously resembles the positive and negative parts of the energy eigenfunctions of a one-dimensional harmonic oscillator. The classical analog to Fig. \ref{OnePartRest} and Fig. \ref{TwoPartRest} is a chain whose points masses oscillate synchronously and a stronger excitation of the chain as a whole corresponds to a higher number of particles in the quantum field.

\subsection{Particle with Non-Zero Wave Number}

\begin{figure}[h!]
\centering
\resizebox{0.99\textwidth}{!}{
  \BW{\includegraphics{OnePartKOne-BW-v2.eps}}
  \color{\includegraphics{OnePartKOne-COL-v2.eps}}
}
\caption{The first excited state of the harmonic chain with wave number $k=1$ corresponds to a standing wave with two nodes in classical mechanics. It can be interpreted as one quantum particle in an energy eigenstate of the lowest non-zero kinetic energy.}
\label{OnePartKOne}
\end{figure}

Now we turn to excitations of the chain with non-zero wave number. For example, the state $\Psi = a_1^\dagger \Psi_0$ (which we can also write as $\Psi_{\{\nu_k\}}$ where $\nu_k$ is equal to $1$ for $k=1$ and zero else) is represented in Fig. \ref{OnePartKOne} for $N=15$. Its classical analog is a standing wave in the oscillator chain. The location of the nodes near $n = 6$ and $n = 13$ can be explained by transforming back to the $q_n$ coordinates: In our state $a_1^\dagger \Psi_0$ the oscillating system is excited along the $\tilde Q_1$ coordinate, since this is how we defined the creation operators in (\ref{EqHnew3}). Using (\ref{EqQtildeq}) we see that an oscillation along $\tilde Q_1$ transforms back to an oscillation along the vector $\vec q$ with
\begin{equation}
q_n = \frac{1}{2}\left((1+i) f_n^{(1)} + (1-i) f_n^{(-1)}\right)= \frac{1}{\sqrt{N}}\left(\cos\frac{2\pi n}{N} + \sin \frac{2\pi n}{N}\right),
\end{equation}
which is zero for $n = 3N/8 \approx 5.6$ or for $n = 7N/8 \approx 13.1$.

A graph of the state $a_{-1}^\dagger \Psi_0$ would look very similar to Fig. \ref{OnePartKOne} with the exception that the standing wave is shifted by an offset of $N/4$ along the $n$-axis compared to the state $a_1^\dagger \Psi_0$. In the particle language of excitations $a_1^\dagger \Psi_0$ and $a_{-1}^\dagger \Psi_0$ would both be one-particle states, which is clear from how it is constructed via one creation operator.

We can also construct the quantum analog to a progressive wave in the chain, for example the state $(a_1^\dagger + i a_{-1}^\dagger) \Psi_0$. But the resulting wave function would not be real-valued anymore. A complex-valued $\Psi$ can be visualized with our method by using different colors, but for the present article we limit ourselves to real-valued $\Psi$ functions which can be visualized \BW{in gray scales}\color{with two colors only}.

\subsection{Localized Particle}

\begin{figure}[h!]
\centering
\resizebox{0.99\textwidth}{!}{
  \BW{\includegraphics{OnePartLocalized-BW-v6.eps}}
  \color{\includegraphics{OnePartLocalized-COL-v6.eps}}
}
\caption{Monte Carlo plot of a local excitation corresponding to one localized particle at the position $n=5$. As opposed to the other states depicted so far, this state is not an energy eigenstate of the harmonic chain, which means that it is not essentially time-independent but would rather evolve into some spread-out oscillation pattern rapidly.}
\label{OnePartLocalized}
\end{figure}

Our next example is a state with one spatially localized particle. Recall that we had expressed $\vec q$ in the orthonormal basis $\{\tilde f^{(k)}_n\}$ in (\ref{EqQtildeq}) and note that the function $\delta_{n,n_0}$ (which is equal to one for $n=n_0$ and else zero) can be expanded in this basis as
\begin{equation} \label{EqDelta}
\delta_{n,n_0} = \sum_k \langle \tilde f^{(k)}_\cdot ,\delta_{\cdot,n_0} \rangle \tilde f^{(k)}_n =  \sum_k \tilde f^{(k)}_{n_0} \tilde f^{(k)}_n.
\end{equation}

Remember that our operators $a^\dagger_k$ create a particle `in the state' $\tilde f^{(k)}_n$. Thus, in analogy to (\ref{EqDelta}) we construct new creation operators
\begin{equation} \label{EqCreatorsLocalized}
b^\dagger_n := \sum_k \tilde f^{(k)}_n a^\dagger_{k}
\end{equation}
and we expect $b^\dagger_n$ to create a particle localized in space at the position $n$. Fig. \ref{OnePartLocalized} is a visualization of the state
\begin{equation} \label{EqPsiLoc}
\Psi = b^\dagger_{5} \Psi_0.
\end{equation}
It is apparent that the sign of $\Psi(\vec q)$ depends only on the sign of the single coordinate $q_5$. This is exactly what one should expect from a quantum oscillator which has $N = 11$ (in the case of Fig. \ref{OnePartLocalized}) degrees of freedom, but which is excited only along the fifth coordinate. Yet the effect of the coupling between the point masses is also visible: At least at the positions $n=4$ and $n=6$ there is a notable correlation of the polylines' \BW{brightness}\color{color} with those at $n=5$ in the sense that \BW{white (black)}\color{red (blue)} lines dominate where $q > 0$ $(q<0)$. It looks like the excitation of the quantum field is `blurred' around that point. Where does this blur come from, given that $\Psi$ in (\ref{EqPsiLoc}) is constructed in close analogy to the $\delta$ function (\ref{EqDelta}) with its sharp peak?
The deeper reason is that $a^\dagger_{k} \Psi_0$ is not only associated with a wave $\tilde f^{(k)}_n$, which stems from the Fourier transformation and which gives rise to the analogy between (\ref{EqDelta}) and (\ref{EqPsiLoc}), but also with an oscillation amplitude in the direction of the $\tilde Q_k$ coordinate. This amplitude depends on the strength of the harmonic potential in the $\tilde Q_k$ direction, which in turn depends (via the diagonalization of the system's Hamiltonian) on the coupling constant $\gamma$ between the mass points. Only in the trivial case $\gamma = 0$ the potential is radially symmetric in the space spanned by the $q_n$ or the $\tilde Q_k$ coordinates, which implies that the oscillation amplitudes of excited states don't depend on $n$ or $k$, so that the analogy between (\ref{EqDelta}) and (\ref{EqPsiLoc}) perfectly holds. In this case the `blurring' of the particle's position in Fig. \ref{OnePartLocalized} vanishes, which is also clear from the fact that for $\gamma = 0$ there is not even the concept of two `neighboring' point masses built into the physical system.

\subsection{Two Localized Particles}

\begin{figure}[h!]
\centering
\resizebox{0.99\textwidth}{!}{
  \BW{\includegraphics{TwoPartLocalized-large-dist-BW-v3.eps}}
  \color{\includegraphics{TwoPartLocalized-large-dist-COL-v3.eps}}
}
\caption{A state with two particles, one of them localized at $n=3$ and the other at $n=8$. While the peaks in excitation at those two locations are clearly visible, the detailed structure of the plot may seem chaotic at first glance. A closer inspection reveals a pattern in the \BW{white and black}\color{red and blue} lines which becomes much more evident in Fig. \ref{TwoPartLocalized-small-dist} below.}
\label{TwoPartLocalized-large-dist}
\end{figure}

\begin{figure}[h!]
\centering
\resizebox{0.99\textwidth}{!}{
  \BW{\includegraphics{TwoPartLocalized-small-dist-BW-v2.eps}}
  \color{\includegraphics{TwoPartLocalized-small-dist-COL-v3.eps}}
}
\caption{Two particles are localized very closely to each other at $n=5$ and $n=6$, respectively. The wave function assumes relatively large positive values for configurations of the harmonic chain where the atoms at $n=5$ and $n=6$ are strongly and synchronously displaced (i.e. where $q_5$ and $q_6$ are either both positive or both negative and of relatively large absolute value). It assumes negative values especially for configurations where the same two atoms are rather weakly and asynchronously displaced (i.e. $q_5$ is slightly positive and $q_6$ is slightly negative, or vice versa).}
\label{TwoPartLocalized-small-dist}
\end{figure}

\begin{figure}[h!]
\centering
\resizebox{0.99\textwidth}{!}{
  \BW{\includegraphics{TwoPartLocalized-no-dist-BW-v3.eps}}
  \color{\includegraphics{TwoPartLocalized-no-dist-COL-v3.eps}}
}
\caption{Two particles are in the same state, localized at $n=5$. The wave function $\Psi$ has two changes of sign along the $q_5$ axis: $\Psi$ is negative where $q_5 \approx 0$ and it is positive where the absolute value of $q_5$ is sufficiently large. This pattern resembles the second excited state of a one-dimensional quantum harmonic oscillator.}
\label{TwoPartLocalized-no-dist}
\end{figure}

In our final example we consider two localized particles, which are described by the state
\begin{equation} \label{EqPsi2Loc}
\Psi = b^\dagger_{n_1} b^\dagger_{n_2} \Psi_0.
\end{equation}
This case is particularly interesting because we can visualize how two bosons move gradually from two different states into one and the same. The Figures \ref{TwoPartLocalized-large-dist}, \ref{TwoPartLocalized-small-dist} and \ref{TwoPartLocalized-no-dist} contain Monte-Carlo visualizations for different values of the pair $(n_1,n_2)$, namely for $(3,8)$, $(5,6)$ and $(5,5)$, respectively. The two particles in Fig \ref{TwoPartLocalized-large-dist} are relatively far apart and each of the two localized excitations resemble the localized one-particle state from Fig. \ref{OnePartLocalized}, even though the phase structure of the $\Psi$ function (i.e. \BW{white}\color{red} lines vs. \BW{black}\color{blue} lines) is much more complicated. In Fig. \ref{TwoPartLocalized-no-dist}, on the other hands, two particles are located at the same point in space. Note the similarity of this plot with the two-particle state shown in Fig. \ref{TwoPartRest}: Both exhibit the pattern of an harmonic oscillator's second excited state (a function with two changes of sign), but in this case the excitation is limited to a small neighborhood of the point $n=5$. Fig. \ref{TwoPartLocalized-small-dist} shows how the transition between the two extreme cases (two particles far apart vs. localized in the same point) comes about.

\section{Discussion of Results} \label{SecDiscussion}

We have presented a new way of visualizing states of a quantum field and we conclude this article with a discussion of advantages and limitations of the new method, as well as an outlook for future investigation.

We hope that our visualizations can help students of QFT and solid state physics to develop a better intuitive understanding of some of their fundamental concepts. Compared to the method presented by Johnson and Gutierrez \cite{JohnsonGutierrez} we feel that our graphs are a slightly more direct representation of the quantum state, since they indicate the value of $\Psi$ for (a randomly chosen set of) points in the state space directly instead of showing a projection of $\Psi$ to the $q_n$ coordinates. For example, the spatial `blurriness' of a localized particle is clearly visible in Fig. \ref{OnePartLocalized} but not in the analogous Fig. 23 of Johnson's and Gutierrez' paper.

An inherent limitation of our method is that the graphs are not exactly reproducible. Due to the random nature of the Monte-Carlo approach two graphs generated with the exact same parameters will look similar, but not exactly the same. Also, the graphs are quite sensitive to parameters chosen in the Monte-Carlo algorithm. For example, if the number of random points chosen is too small, then none of these points might be in a region where $|\Psi|$ is large and the resulting graph will look quite noisy and without any clear structure. Whether that number is 'small' also depends on the number $N$ of point masses. The volume from which interesting points can be chosen (e.g. the set of all $\vec q \in \mathbb{R}^N$ with $|q_j| < 1$ for all $j = 1,\dots,N$) grows exponentially with $N$. Consequently, the necessary number of random points and the computation time of the algorithm also grow exponentially with $N$. Finally, for too `complex' states (e.g. those with more than just a few particles of different wave numbers) the graphs become too busy to be readable.

In principle, our method could be extended to two-dimensional fields, where the polylines from the plots in the present paper would be replaced by overlapping surfaces in a 3D plot. In order to avoid a very busy and confusing graph, the number of surfaces drawn would have to be limited strictly to just a few with the very highest values of $|\Psi|$.

An interesting subject for future work will be the extension of the method presented here to more complex phenomena in quantum fields, for example by considering interactions between fields or by making the shift from bosons to fermions.

\section{Acknowledgments}

The author is grateful to the anonymous referees for thorough reviews and several helpful suggestions which have made this article better.