\section{Introduction}

Commonsense reasoning is an important yet challenging task in artificial intelligence and natural language processing. Take commonsense question answering as an example, given a question and multiple choices, some commonsense knowledge is usually required to make the correct answer from the provided choices. Table~\ref{examples} show some typical commonsense question answering examples extracted from the dataset of commonsenseQA~\cite{talmor2018commonsenseqa}.

\begin{table}
\caption{Examples from CommonsenseQA dataset.}
\label{examples}
\begin{center}
\begin{small}
\begin{tabular}{l p{5cm}}
\toprule
Question: & Where is a good idea but not required to have a fire extinguisher?\\
Choices: & (A)~school bus (B)~boat \textbf{(C)~house} (D)~hospital (E)~school\\
\midrule
Question: & Where can you put a picture frame when it's not hung vertically?\\
Choices: & (A)~art show (B)~wall (C)~newspaper (D)~car \textbf{(E)~table}\\
\bottomrule
\end{tabular}
\end{small}
\end{center}
\vskip -0.25in
\end{table}

Existing commonsense reasoning methods mainly utilize raw texts to conduct the data representation and answer prediction process~\cite{talmor2018commonsenseqa, rajani2019explain}. However, the background knowledge required in the commonsense reasoning task, such as spatial relations, causes and effects, scientific facts and social conventions, are usually not explicitly provided by the text. Therefore, it is difficult to capture such knowledge solely from the raw texts. Some other works propose to leverage knowledge bases to extract related commonsense knowledge~\cite{lin2019kagnet, lv2019graph, kipf2016semi, ye2019align, li2019teaching, ma2019towards}. However, the construction of a knowledge base is expensive, and the contained knowledge is too limited to fulfill the requirement. Furthermore, most commonsense question answering datasets, such as CommonsenseQA, are constructed from an existing knowledge base, e.g., ConceptNet \cite{speer2017conceptnet}. So it is unfair to use the knowledge base in these tasks. To sum up, how to automatically learn commonsense remains a challenging problem in NLP.

Motivated by the fact that images usually contain richer scene information, which can be viewed as an important supplementary resource to perceive for commonsense knowledge, this paper proposes to learn commonsense from images and incorporate such knowledge into the commonsense reasoning process. Take the question `\emph{Where is a good idea but not required to have a fire extinguisher?}' shown in Table~\ref{examples} as an example. Solving this problem requires a strong background knowledge that fire extinguishers are usually equipped in public places, such as hospitals, schools, and school buses. We can see that such background knowledge is not explicitly provided by the raw texts, and meanwhile, too abstract and complex to be extracted by the current language model techniques. In this case, images will help. For example, we could find many images where fire extinguishers appear in these scenes of public places. Therefore, this commonsense knowledge could be learned by perceiving the scene information of these images, and the corresponding question will be well answered. These analyses are in accordance with Minsky's statement in \citet{minsky2000commonsense}, `perhaps a good architecture theory based on multiple representations and multi-modal reasoning would help us to design better systems that allow us to study and understand commonsense reasoning.'

Our approach, named Loire (\textbf{L}earning C\textbf{o}mmonsense from \textbf{I}mages for \textbf{Re}asoning), consists of two stages, i.e.~visual commonsense learning and knowledge-augmented reasoning. In the first stage, a scene layout generation task is conducted on a bi-modal data such as the representative benchmark COCO~\cite{lin2014microsoft}. Firstly, a text encoder Visual BERT (ViBERT for short) is employed to obtain the representation of a caption. ViBERT is then incorporated into the recurrent encoder-decoder structure for the labeled bounding box generation. This module is trained separately by a supervised learning approach, based on the ground-truth bounding boxes of images. In this way, the required visual commonsense knowledge will be encoded in ViBERT. In the following commonsense reasoning stage, the concerned text representations (such as question and answer in commonsenseQA) will be obtained by concatenating ViBERT and a traditional pre-trained language model, e.g. ~BERT. Then the language model is fine-tuned on the commonsense reasoning data, with ViBERT fixed as some prior knowledge. Experimental results on two commonsense reasoning tasks, i.e.~CommonsenseQA and WinoGrande \cite{sakaguchi2019winogrande}, demonstrate that the learnt commonsense from images brings improvements to traditional models, such as BERT fine-tune \cite{devlin2018bert} and RoBERTa fine-tune \cite{liu2019roberta}. We also give some case studies to show how the learned visual commonsense knowledge helps the reasoning process.

To the best of our knowledge, we are the first to propose learning commonsense knowledge from images to facilitate the commonsense reasoning in NLP. The proposed model of using scene layout generation as the supervision demonstrates a preliminary exploration in this direction. Other methods like learning commonsense from retrieved relevant images could also be investigated. We believe this novel approach may provide a new perspective for commonsense reasoning in NLP.

\section{Related Work} \label{related-work}

\subsection{Commonsense reasoning Methods}
There are mainly two kinds of commonsense reasoning methods, knowledge base approach and raw text approach.

Knowledge base approach makes use of the existing knowledge bases \cite{speer2017conceptnet, sap2019atomic} to conduct the commonsense reasoning process. Some methods regard knowledge base as a supplement and integrate extracted knowledge with information from the processed text.
For example, \citet{mihaylov2018knowledgeable} encodes external commonsense knowledge as a key-value memory. \citet{lv2019graph} and \citet{lin2019kagnet} extract knowledge from ConceptNet and Wikipedia to construct graphs, then use Graph Convolutional Network \cite{kipf2016semi} for modeling and inference. Other methods \cite{zhong2018improving, ma2019towards, ye2019align, li2019teaching} use knowledge bases as another corpus for pre-training, and then refining the models on task-specific contents.

Besides extracting knowledge from knowledge bases, some other methods directly learn commonsense knowledge from raw texts. A common way is to use pre-trained language models. Recently, \citet{talmor2018commonsenseqa, da2019cracking, sakaguchi2019winogrande, zhou2019evaluating} have made comprehensive empirical studies and shown that pre-trained language models significantly outperform traditional methods on the task of commonsense reasoning. In addition, \citet{da2019cracking} proves that pre-trained language models have the ability to encode some commonsense knowledge in the embedding space through the attribute classification evaluation. However, they also show that the encoded commonsense knowledge is limited, which could be improved by introducing some supplementary data, like ConceptNet.

Moreover, some methods propose to leverage additional text information/data for better commonsense reasoning. \citet{tandon2018reasoning} uses commonsense knowledge as constraints and large scale web corpus to steer the model away from unlikely predictions.
\citet{rajani2019explain} incorporates the generated explanations into the training of language models for enhancement. \citet{xia2019incorporating} leverages two auxiliary relation-aware tasks to better model the interactions between question and candidate answers. \citet{chalier2020joint} proposes a multi-faceted model of commonsense knowledge statements to capture more expressive meta-properties.

Different from the above approaches, we propose to learn commonsense from images and incorporate this visual knowledge into the following commonsense reasoning process.

\subsection{Bi-modal Language Models}
Recently, some transformer-based bi-modal language models~\cite{Su2019VL, li2019unicoder-vl:, alberti2019fusion, li2019visualbert, tan2019lxmert, lu2019vilbert} have been proposed to tackle with bi-modal reasoning problems in computer vision, such as visual question answering, visual commonsense reasoning, and image retrieval. They first encode image representation and text representation into a shared embedding space, then apply the joint embeddings for downstream reasoning. At first glance, these models are quite similar to ours. However, we should make it clear that they are totally different. The purpose of a bi-modal language model is to capture a cross-modal alignment between image and text to benefit the bi-modal task, which is only available when both image and text data are provided as input simultaneously. That is why they are usually popular in bi-modal scenarios like VQA. If we want to apply these models to commonsense reasoning in NLP, how to find corresponding images to the question, and how to employ the joint embeddings to the downstream NLP reasoning tasks is still unclear. Our model also adopts image data as a supplementary, but the modeling approach is different from bi-modal language models. We first encode the visual commonsense knowledge into ViBERT by the upstream layout generation process on a bi-modal data, then apply ViBERT as fixed prior knowledge to fine-tune the pre-trained language models for the downstream NLP reasoning tasks.

\section{Visual Commonsense Knowledge} \label{commonsense-in-image}

\begin{figure}
\centering
\includegraphics[width=0.31\columnwidth]{figure/eat_img_163540.pdf}
\includegraphics[width=0.31\columnwidth]{figure/eat_img_203416.pdf}
\includegraphics[width=0.31\columnwidth]{figure/eat_img_248478.pdf}
\caption{Images and the associated bounding boxes from COCO with captions similar to `a woman eats in the restaurant'.}
\label{restaurant_example}
\vskip -0.25in
\end{figure}

Images are made up of individual pixels, which are detailed but sometimes noisy. Therefore, how to extract useful commonsense knowledge from images remains a challenging problem. Inspired by the knowledge base in NLP, where knowledge is usually represented as a triple to demonstrate the relation between two entities, we focus on the attributes and relations of the objects in images. Clearly, such information can be well captured by the scene layout. Take the sentence `a woman eats in the restaurant' as an example. Images related to this sentence are shown in the Figure~\ref{restaurant_example}. We can see that the scene layouts of these images, including bounding boxes and labels, contains a lot of useful information for commonsense reasoning:

(1) Size attributes and relations can be easily obtained by the bounding boxes in images. For instance, the bounding boxes of tableware, e.g. ~fork, cup, spoon are always smaller than the bounding boxes of the dining table.

(2) Position can be accurately captured by the coordinate of each bounding box, to help understand some abstract commonsense. For instance, through the relative positions between the bounding boxes of person and table, one can figure out what "next to" means. Besides, since the bounding boxes of person and table are always close in the eating scene, one can learn that if a person is eating, he will be next to the table instead of standing far away, which provides some detailed information for the abstract word `eating'.

(3) Co-occurrence relations of objects are expressed by the labels of bounding boxes. For instance, images of `a woman eats in the restaurant' often contain labels of table, chair, person, food and tableware. So from the co-occurrence of these objects, one can infer that it is in a dinner or restaurant scenario, which offers rich context information for the abstract word `eating'.

From the above analysis, images usually contain rich scene information, such as size, position and co-occurrence relations, which are useful for understanding the commonsense knowledge hidden in language. So we propose to learn such visual commonsense knowledge and incorporate them into the commonsense reasoning models in NLP.

\section{Our Approach: Loire}
Now we introduce Loire, which includes two stages, i.e.~visual commonsense learning and knowledge-augmented reasoning.

\begin{figure}
\begin{center}
  \centerline{\includegraphics[width=0.9\columnwidth]{figure/recurrent-layout-generation.pdf}}
  \caption{The recurrent structure of the visual commonsense learning stage.}
\label{overall-for-layout}
\end{center}
\vskip -0.4in
\end{figure}

\subsection{Visual Commonsense Learning} \label{scene-layout-generation}
The visual commonsense learning stage is conducted on bi-modal data, like the typical image caption data COCO. For a given image, the required scene layout is generated by a sequence-to-sequence approach, shown in Figure~\ref{overall-for-layout} and ~\ref{onestep-for-layout}. This module consists of a text encoder, namely ViBERT, to map the input
sentence to a latent representation, a layout encoder to encode the current generated scene layout, and a bounding box decoder to generate the next bounding box and its label.

Specifically, we make the following notations. Let the input image caption be $S = \{w_1, w_2, \dots, w_L\}$, where $w_i$ stands for the $i$-th word in the sentence, and $L$ is the sentence length. The output is a set of labeled bounding boxes $B_{1:T} = \{B_1, ..., B_T \}$, with each labeled bounding box $B_t$ contains the position, size and category label of a corresponding object at the $t$-th step. So we denote $B_t = (b_t,l_t)$, where $b_t = [b_t^x, b_t^y, b_t^w, b_t^h] \in \mathbb{R}^4$ stands for 2-dimensional coordinates, width and height, respectively. $l_t \in \{0,1\}^{C+1}$ is a one-hot vector to indicate the category label for an object, and the additional $C+1$ class is defined as a special indicator for the end of generation.

\subsubsection{ViBERT: Text Encoder}
The text encoder ViBERT is fine-tuned from BERT, which is a popular pre-trained language model introduced in \citet{devlin2018bert}. The network structure is a typical transformer-based architecture containing multiple transformer blocks of multi-headed scaled dot product attention and fully connected layers~\cite{vaswani2017attention}. It has been proven to be effective in many natural language processing tasks.

To adapt to the setting of BERT, the image caption is preprocessed as follows. The special tokens `[CLS]' and `[SEP]' are inserted into the beginning and the end of the sentence, to obtain $S = \{w_0, w_1, ..., w_{L+1}\}$, where $w_0, w_{L+1}$ stands for $\textrm{[CLS]}$ and $\textrm{[SEP]}$, respectively. After that, each word $w_i$ is mapped to its word embedding vector $e_i^S$ by ViBERT, so that $e(S) = \{e_0^S,e_1^S,...,e_{L+1}^S\}$.
With BERT, the output of `[CLS]' from the last transformer layer is fed into a pooler layer to obtain the representation of the whole sentence $e^S$,
\begin{equation}
    e^S = \textrm{tanh}(f(e_0^S)),
\end{equation}
where $f$ is a single-layer perceptron.

\subsubsection{Layout Encoder}
At each time step $t$, a layout encoder is utilized to encode the state of the current generated layout $B_{0:t-1}$. Specifically, we construct a layout matrix $I_{t-1} \in \{0,1\}^{C \times W \times H}$, where $W, H$ are width and height of this layout respectively. The value of $i_{lwh}$ in $I_{t-1}$ indicates whether the bounding box of object $l$ covers the pixel at coordinate $[w, h]$. A blank layout without any object is used to initialize $B_0$. A layout encoder takes layout matrix and previous layout representation as inputs, and uses a convolutional GRU architecture to output the representation of the current layout $e^I_t$ as follows:
\begin{equation}
  e^I_t = \textrm{ConvGRU}(I_{t-1}, e^I_{t-1}).
\end{equation}

\begin{figure}
\begin{center}
  \centerline{\includegraphics[width=0.9\columnwidth]{figure/one-step-generation.pdf}}
  \caption{An illustration of the $t$-th layout generation.}
\label{onestep-for-layout}
\end{center}
\vskip -0.4in
\end{figure}

\subsubsection{Bounding Box Decoder}
At each time step $t$, a bounding box decoder is used to predict the labeled bounding box of next object, based on the caption representation $e^S$ from ViBERT and the current layout representation $e_t^I$ from the layout encoder. Specifically, we decompose the conditional joint bounding box probability as $p(b_t,l_t|s, B_{0:t-1}) = p(l_t|S, B_{0:t-1}) p(b_t|S, B_{0:t-1},l_t)$. The decoder firstly samples a class label $l_t$ according to $p(l_t|S, B_{0:t-1})$:
\[
  p(l_t|s,B_{0:t-1}) = \textrm{Softmax}(g(u_t, c_t)),
\]
\[
  u_t^l = \phi^l(e^I_t, e^S),\;\;
  c_t^l = \varphi^l([u_t^l; l_{1:t-1}], e(S)),
\]
where $g$ is a two-layer perceptron, $\phi^l$ is a Convolution network~\cite{xu2015show} with spatial attention on $e^I_t$, and $\varphi^l$ is a text-based attention module~\cite{luong2015effective}, which is used to focus on different parts of the caption.

After that, the decoder tries to find out $b_t$ for object $l_t$ based on $p(b_t|S,B_{0:t-1},l_t)$, which is obtained by a regression network $\theta$ with $\hat{b}_t = (\hat{x}_t, \hat{y}_t, \hat{w}_t, \hat{h}_t) = \theta(c_t^b, u_t^b)$. The parameters $c_t^b$ and $u_t^b$ are represented similarly to $u_t$ and $c_t$. That is,
\[
   u_t^b = \phi^b (e^I_t, c_t^b), \;\;c_t^b = \varphi^b([u_t^l; l_t], e(S)),
\]
where $\phi^b$ is an image-based attention module to find an appropriate position, and $\varphi^b$ is another text-based attention module but focuses more on the contents related to the current object.

\subsubsection{Training}
To reduce the expensive training ViBERT from scratch, we initialize ViBERT with the parameter weights of BERT$_{BASE}$ released by Google~\footnote{\url{https://github.com/google-research/bert}}. Then the scene layout generator can be trained by minimizing the negative log-likelihood of the ground-truth object labels and the mean-square error of the ground-truth bounding box coordinates as follows:
\begin{equation}
  \mathcal{L}_{layout} = \sum_t\big(||\hat{b_t}-b_t^*||_2 - \log p(l_t^*) \big),
\end{equation}
where $b_t^*$ and $l_t^*$ stands for the ground-truth bounding box and label, respectively. As for the generation order, we have observed that the model is difficult to converge with unfixed order, which may be caused by some dependencies among different bounding boxes. So we follow the existing image generation methods and simply fix the order from bottom to top, left to right.

It should be noted that although we use BERT as a text encoder on image captions, we do not optimize the objective of the language model, i.e. ~ the masked language model (MLM) objective. This is to avoid the possibility that the improvement of downstream reasoning task is due to the use of more text data, instead of visual commonsense knowledge from images. In our experiments, we have conducted some ablation studies to validate this point.

\subsection{Knowledge-Augmented Reasoning} \label{knowledge-argumented-reasoning}
After using scene layout generation to encode visual commonsense knowledge into ViBERT, we can apply ViBERT as a fixed prior to enhance the downstream commonsense reasoning tasks.

Here we use commonsenseQA as an example to demonstrate our method. For a given question $q_i \in Q$, where $Q$ is the question set, and its candidate answers $A_i = \{a_i^1, \dots, a_i^n\}$, where $n$ denotes the number of choices, a common existing method is to first concatenate question and each candidate answer to a raw representation $[q_i;a_i^j]$. Then a pre-trained language model is applied to obtain a semantic representation, denoted as $E_{i,j}^{(1)} = \textrm{LM}([q_i;a_i^j])$. In our method, ViBERT is applied on the raw representation $[q_i;a_i^j]$ to obtain a image scene-aware text representation, denoted as $E_{i,j}^{(2)} = \textrm{ViBERT}([q_i;a_i^j])$. Since the two representations are not always in the same space, we use a projection matrix $M$ to project $E_{i,j}^{(2)}$ to the space of $E_{i,j}^{(1)}$. After that, they are simply concatenated and fed into a linear layer to compute the probability $p(a_i^j|p_i)$ as follows.
\[
  \textrm{score}(a_i^j) = h(E_{i,j}^{(1)}; M^T E_{i,j}^{(2)}]),
\]
\[
  p(a_i^j|p_i) = \textrm{Softmax}(\{\textrm{score}(a_i^j)\}_j),
\]
where $h$ is a simple linear layer for classification, and the parameters of both language model and the linear layer will be fine-tuned on the downstream commonsense reasoning task. Specifically, in the training process, the objective is to minimize the negative log-likelihood of ground-truth answers $a_i^*$ as follows. After that, the choice with the highest score will be selected as the answer.
\begin{equation}
  \mathcal{L}_{qa} = -\sum_i \log p(a_i^*|q_i).
\end{equation}

\section{Experiments} \label{experiment}
This section demonstrates our experiments on two commonsense reasoning tasks, i.e.~comonsense question answering and pronoun resolution.

\subsection{Datasets}

\textbf{CommonsenseQA}\footnote{\url{https://www.tau-nlp.org/commonsenseqa}}~\cite{talmor2018commonsenseqa} is a typical commonsense question answering dataset, which consists of 12,102 natural language questions generated from ConceptNet. It covers various types of commonsense knowledge, including spatial, causal, social, and activity, etc. Each question has five candidate answers. Table~\ref{csqa-examples} shows 3 question-answering examples. In our experiments on this dataset, we use the official random-split setting for fair comparisons with the reported results on CommonsenseQA’s leaderboard.

\textbf{WinoGrande}\footnote{\url{https://leaderboard.allenai.org/winogrande/submissions/get-started}}~\cite{sakaguchi2019winogrande} is a challenging pronoun resolution dataset extended from the original Winograd Schema Challenge \cite{levesque2012winograd}. The task is about resolving a pronoun (represented as a blank line) to one of its two probable co-referents in the sentence. For this task, each sentence is treated as a fill-in-the-blank question with binary choices. The line in the sentence is replaced by each option, and the model is required to provide the likelihood for the two resulting sentences for determination. In the training set of WinoGrande, there are five different sizes, i.e.~XS(160), S (640), M (2,558), L (10,234) and XL (40,398). We experiment on all the five sizes and report their results for analysis.

\subsection{Experimental Settings}
For the upstream scene layout generation module, we train our \textbf{ViBERT} on 2 Nvidia K80 GPUs with a batch size of 32 for 15 epochs. The learning rate is $5e^{-5}$, and the optimizer is Adam with StepLR schedule, where the step size is 3 and $\gamma$ is 0.8. In the training process, the bi-modal data COCO~\cite{lin2014microsoft} is used to train our layout generation model. COCO consists of 123,287 images over 80 object categories, and each image is associated with instance-wise annotations and 5 image captions. For better training, we ignore small objects and filter images with more than 20 objects. This leaves us 119,146 images. We use the official train and validation splits, and set a max sequence length as 128.

For the downstream commonsense reasoning module, we choose BERT and RoBERTa as our baseline models, which are the fundamental and competitive models for NLP tasks.

\textbf{BERT}~\cite{talmor2018commonsenseqa} is a powerful contextualized word representation model and has been proven helpful in many NLP tasks. We apply uncased BERT$_{\textrm{BASE}}$ to downstream commonsense reasoning tasks by encoding each question and its candidate answers as a series of delimiter-separated sequences, i.e. ~`[CLS] question [SEP] choice [SEP]' for CommonsenseQA and `[CLS] segment1 [SEP] option segment2 [SEP]' for WinoGrande. Then the representation of `[CLS]' is then fed into a BERT-Pooler and converted to predictions by a linear classification layer.

\textbf{RoBERTa}~\cite{liu2019roberta} is similar to BERT, but is usually pre-trained with a larger amount of training data and different techniques such as dynamic masking. Besides RoBERTa$_{\textrm{BASE}}$, we also compare with a fine-tuned RoBERTa$_{\textrm{LARGE}}$ following the implementation released in fairseq~\footnote{\url{https://github.com/pytorch/fairseq/tree/master/examples/roberta/commonsense_qa}}. And according to fairseq, we prepend a prefix of Q: to the question and A: to the answer for CommonsenseQA, which was found to be helpful.

\textbf{Loire} By using BERT and RoBERTa as a language model for text, we concatenate the representations from ViBERT and the pre-trained language model, and obtain two versions of our model, denote as \textbf{Loire-BERT} and \textbf{Loire-RoBERTa}, respectively. Since ViBERT is a static feature extractor and doesn't need to be fine-tuned in the downstream reasoning tasks, our running time is similar to the baselines without extra time cost.

We train all models on 2 Nvidia K80 GPUs using AdamW~\cite{loshchilov2018fixing} with WarmupLinearSchedule approach~\cite{he2016deep} for optimization, where the warmup percentage is set to 0.1 and 0.05 for BERT and RoBERTa, respectively.
We use grid-search for hyper-parameters tuning. The learning rate, number of epochs and batch-size are chosen from $\{1, 2\} \times e^{-5}$, $\{3, 5, 8\}$, and $\{8, 16, 32\}$.
The best development set accuracy from 5 random restarts of fine-tuning is reported, with the standard deviation. The best models on the development dataset are then submitted to the official private test dataset to return the test results. All our code and data are publicly available at \url{https://github.com/VickiCui/Loire}.

\subsection{Experimental Results}
On the dev set, the accuracy of the layout generation for label prediction is 63.4\%, and the mean square error for bbox prediction is 0.015 (the coordinates of bbox have been standardized between 0 and 1). This shows that the layout generator has a good performance and can generate good quality scene layouts, and the model does learn the corresponding knowledge.

\begin{table}
\caption{Results on CommonsenseQA (\%), where `*' indicates the reported result from the leaderboard.}
\label{csqa_result}
\begin{center}
\begin{small}
\begin{adjustbox}{max width=1.\linewidth}
\begin{tabular}{lccc}
    \toprule
    Model &  Dev Acc. & Dev Avg.& Test Acc.\\
    \midrule
    \citet{ott2019fairseq}  & - & - & \textbf{72.1$^*$} \\
    \midrule
    RoBERTa$_{\textrm{LARGE}}$  & 77.47 & 76.65$\pm{0.58}$ & 71.58 \\
    Loire-RoBERTa$_{\textrm{LARGE}}$  & \textbf{77.94} & \textbf{77.56$\pm{0.28}$} & \textbf{71.93} \\

    RoBERTa$_{\textrm{BASE}}$  & 65.47 & 64.96$\pm{0.62}$ & 59.82 \\
    Loire-RoBERTa$_{\textrm{BASE}}$   & 66.67 & 66.12$\pm{0.47}$ &60.61 \\

    BERT$_{\textrm{BASE}}$  & 59.71 & 58.95$\pm{0.65}$ & 53.00$^*$ \\
    Loire-BERT$_{\textrm{BASE}}$  & 61.19 & 60.07$\pm{0.58}$ & 54.91  \\
    \midrule
    Human & - & - & 88.00 \\
    \bottomrule
\end{tabular}
\end{adjustbox}
\end{small}
\end{center}
\vskip -0.15in
\end{table}

\begin{table
\caption{Results on WinoGrande with 5 training sizes, where `*' indicates the reported result from the leaderboard.}
\label{wsc_result}
\vskip -0.25in
\begin{center}
\begin{small}
\begin{adjustbox}{max width=1.\linewidth}
\begin{tabular}{lccccc}
    \toprule
    Model & XS & S & M & L & XL\\

    \midrule
     & \multicolumn{5}{c}{Dev Acc. (\%)} \\
    BERT$_{\textrm{BASE}}$ & 50.76 & 51.61 & 52.81 & 55.26 & 60.19   \\
    Loire-BERT$_{\textrm{BASE}}$ & \textbf{51.61} & \textbf{52.34} & \textbf{53.9} & \textbf{56.74} & \textbf{61.50}   \\
    RoBERTa$_{\textrm{BASE}}$ & 51.72 & 54.71 & 57.91 & 62.52 & 67.94   \\
    Loire-RoBERTa$_{\textrm{BASE}}$ & \textbf{53.26} & \textbf{55.18} &\textbf{ 58.93} & \textbf{64.09} & \textbf{69.21 }  \\
    RoBERTa$_{\textrm{LARGE}}$ &  52.40 & 61.95 & 68.67 & 75.14 & 79.08   \\
    Loire-RoBERTa$_{\textrm{LARGE}}$ &  \textbf{52.64} & \textbf{63.06} & \textbf{70.40} & \textbf{76.56} & \textbf{81.06}   \\

    \midrule
     & \multicolumn{5}{c}{Test Acc. (\%)} \\
    BERT$_{\textrm{BASE}}$ & 49.75 & \textbf{49.75} & 49.01 & 51.50 & 54.73   \\
    Loire-BERT$_{\textrm{BASE}}$ & \textbf{49.86} & 49.29 & \textbf{52.07} & \textbf{53.88} & \textbf{59.54}   \\
    RoBERTa$_{\textrm{BASE}}$ & 50.93 & 52.01 & \textbf{57.67} & 61.35 & 65.42 \\
    Loire-RoBERTa$_{\textrm{BASE}}$ & \textbf{53.42} & \textbf{53.42} & 56.82 & \textbf{62.31} & \textbf{67.12}\\
    \citet{levesque2012winograd} &  50.37 & 58.63 & 67.57 & 74.70 & \textbf{79.12}   \\
    \citet{yang2020g} &  \textbf{55.04} & 62.37 & 66.72 & 74.19 & 78.21   \\
    Loire-RoBERTa$_{\textrm{LARGE}}$ &  53.14 & \textbf{63.27} & \textbf{70.51} & \textbf{76.12} & 77.99   \\
    \bottomrule
\end{tabular}
\end{adjustbox}
\end{small}
\end{center}
\vskip -0.25in
\end{table}

Table~\ref{csqa_result} shows the experimental results on CommonsenseQA. From the results, we can see that our approach leads to a 1.91\%, 0.79\% and 0.35\% improvement in terms of accuracy on the test set, as compared with BERT$_{\textrm{BASE}}$, RoBERTa$_{\textrm{BASE}}$ and RoBERTa$_{\textrm{LARGE}}$ respectively. Similar results have been observed on the development set. Besides, the standard deviation of several random results on the development set becomes smaller when using Loire, which demonstrates better stability. Someone may argue that the improvement is marginal as compared with RoBERTa$_{\textrm{LARGE}}$, and our result is worse than the best result of RoBERTa$_{\textrm{LARGE}}$ on the leaderboard \cite{ott2019fairseq}. It should be noted that the best result of RoBERTa$_{\textrm{LARGE}}$ on the leaderboard is based on validation performance after 100 trials. However, we only conducted five trials in our experiments due to our limited computing resources. The purpose of this paper is to propose a new perspective of learning commonsense from the image, rather than achieving a SOTA result. We can clearly see some improvement from the comparison with the baseline models. It is acceptable that when using more complicated language models, the effect of visual knowledge will be weakened. However, there are indeed some methods to improve the current results, which will be investigated in our future work. For example, we have filtered out small objects to make the training easier, which may result in insufficient details. Besides, the adopted bi-modal data COCO is very limited, with only 80 categories of objects. On the one hand, the coverage of the commonsense may be restricted. On the other hand, the layouts generated by our model may not be very accurate for some objects. For instance, the generated layout of `laundry' is `a suitcase' since COCO does not contain clothes in our case study. We plan to employ larger data such as Visual Genome~\cite{krishna2017visual} to tackle this problem.

Table~\ref{wsc_result} shows the experimental results on WinoGrande. Specifically, five models are trained on five different training data sizes separately, and the development set and test set are identical for all models. As for the accuracy of the development set, We can see that Loire achieves consistent performance improvements across different sizes of training data, as compared with both BERT$_{\textrm{BASE}}$, RoBERTa$_{\textrm{BASE}}$ and RoBERTa$_{\textrm{LARGE}}$. While for the test accuracy (\citet{levesque2012winograd} and \citet{yang2020g} are two test results of RoBERTa$_{\textrm{LARGE}}$ from the leaderboard), except for a few ones, Loire consistently outperforms the corresponding baselines on across different sizes of training data. These results show the effectiveness of incorporating visual scene knowledge for commonsense reasoning.

\subsection{Ablation Study}

\begin{table}
\caption{Accuracy (\%) of different models on CommonsenseQA development set.}
\label{ablation_study}
\begin{center}
\begin{small}
\begin{tabular}{lccc}
    \toprule
      Model. & Dev Acc. & Dev Avg\\
    \midrule
    BERT$_{\textrm{BASE}}$  & 59.71 & 58.95$\pm{1.03}$\\
    \hspace{0.9em}+BERT$_{\textrm{BASE}}^*$  & 59.89 & 59.12$\pm{0.65}$ \\
    \hspace{0.9em}+BERT$_{\textrm{CAPTION}}$  & 60.29 & 59.47$\pm{0.60}$\\
    \hspace{0.9em}+ViBERT (ours)  & \textbf{61.19} & \textbf{60.07$\pm{0.58}$} \\
    \bottomrule
\end{tabular}
\end{small}
\end{center}
\vskip -0.25in
\end{table}

In order to validate that the performance improvement owes to the introduction of learned visual commonsense knowledge, rather than using more parameters or data, we conduct the following ablation studies on CommonsenseQA. The results are shown in Table~\ref{ablation_study}, where `+ViBERT' denotes Loire.

Firstly, we study whether the improvement owes to the use of additional parameters. To this end, we compare with the BERT$_{\textrm{BASE}}$ concatenated with freeze BERT$_{\textrm{BASE}}^*$ features, in which the parameters are set to be the same as BERT$_{\textrm{BASE}}$+ViBERT. From the results, we can see that, under the same setting, the accuracy of BERT$_{\textrm{BASE}}$ concatenated with freeze BERT$_{\textrm{BASE}}^*$ features is $59.89\%$ on dev set, which is worse than ours.

Then we study whether the improvement owes to the use of additional text data, i.e.~captions in COCO. We first fine-tune a BERT$_{\textrm{BASE}}$ model on COCO captions with MLM objective, denoted as BERT$_{\textrm{CAPTION}}$. Then we concatenate it with BERT$_{\textrm{BASE}}$, to perform a similar downstream fine-tuning as in Loire-BERT$_{\textrm{BASE}}$. We also randomly initialized the model 5 times. The best dev result is $60.29\%$, which is worse than Loire.

In summary, these ablation studies prove that the commonsense knowledge learned form images, rather than the introduction of more parameters or text data, is responsible for the improvements.

\subsection{Case Study}

\begin{table}[t]
\caption{Case study examples from the dev set of CommonsenseQA.}
\label{csqa-examples}
\begin{center}
\begin{small}
\begin{tabular}{l p{5cm}}
\toprule
Question1: & The man got a pail to catch the draining motor oil, where was he likely doing this at home?\\
Choices1: & \textbf{(A)~garage} (B)~hardware store (C)~utility room (D)~wishing well (E)~laundry\\
\midrule
Question2: & Where would a person be doing when having to wait their turn?\\
Choices2: & (A)~have patience (B)~get in line (C)~sing \textbf{(D)~stand in line} (E)~turn left\\
\midrule
Question3: & Where would you find magazines along side many other printed works?\\
Choices3: & (A)~doctor \textbf{(B)~bookstore} (C)~market (D)~train station (E)~mortuary\\
\bottomrule
\end{tabular}
\end{small}
\end{center}
\vskip -0.25in
\end{table}

To understand what type of commonsense knowledge is learned by Loire, we analyze the relations between question concept and answer concept in CommonsenseQA according to ConceptNet. For the part of the questions that are done right by our model but wrong by the text-only model, which can be seen benefit from images, the top three relation types are AtLocation (36.4\%), Causes (12.7\%) and RelatedTo (8.5\%). These relationships can indeed be expressed through the scenes shown in the images. So this is accordant with our motivation, and the introduction of images can indeed play a complementary role. For complete statistics of relation types, please see Appendix A.

Table~\ref{csqa-examples} gives three examples in the development set of CommonsenseQA that benefit from visual commonsense knowledge. To better understand how visual commonsense helps, we generate the layout for each pair of question and choice by the trained upstream layout generator. Figure~\ref{case-study-layout} shows the layouts of Question1 and its choices, and others can be found in Appendix B due to space limitations.

Take the first question as an example, language models mainly rely on word co-occurrence or semantics for modeling, so they are easy to wrongly choose 'utility room' as the answer. That is because it is difficult to capture the commonsense of `got a pail to catch the draining motor oil in garage' from the pure language. From Figure~\ref{case-study-layout}, we can see that the layout of question, the correct answer `garage' and the wrong answer `utility room' are 'a person' with `a truck', `cars', and 'chairs' and `old televisions', respectively. That is to say, we can learn from images that `got a pail to catch the draining motor oil' usually happen with the scene that a person is with a truck. By encoding this knowledge into ViBERT, it is easy for the language model to connect the similarity between `truck' and `cars', so Loire is able to choose the correct answer 'garage', instead of 'utility room'.

\begin{figure}
\begin{center}
  \centerline{\includegraphics[width=0.9\columnwidth]{figure/layout_example.pdf}}
  \caption{Scene layout of the first example in Table~\ref{csqa-examples}.}
\label{case-study-layout}
\end{center}
\vskip -0.4in
\end{figure}

\section{Conclusion} \label{conclusion}
In this paper, we propose a novel two-stage pipeline approach Loire to learn commonsense from images. In the first stage, a text representation model ViBERT is trained in the bi-modal sequence-to-sequence approach for scene layout generation on COCO. Therefore, visual commonsense knowledge like spatial relations will be encoded in ViBERT by the supervision of caption and image layout. After that, ViBERT is concatenated with a pre-trained language model to perform a knowledge-augmented reasoning process. Experimental results show that Loire outperforms the current state-of-the-art language models BERT and RoBERTa on two NLP commonsense reasoning tasks, i.e.~commonsense question answering data CommonsenseQA and pronoun resolution data WinoGrande. The ablation and case study further show that the improvements are truly owing to the learned visual commonsense knowledge, and how this knowledge helps the NLP reasoning process.

The current approach is a preliminary study on the proposed direction of using images to automatically learn commonsense knowledge to facilitate the NLP reasoning tasks, which could be modified from the following aspects to further improve the empirical performances. Firstly, larger bi-modal data could be employed to learn more commonsense required in the reasoning task. Secondly, other bi-modal methods instead of training ViBERT by the supervision of scene layout generation may be investigated. Thirdly, how to design intrinsic evaluation to help to understand what is learned by Lorie is still challenging and will be considered in the future.

\section*{Acknowledgement}
This work was supported by the National Key R\&D Program of China (2020AAA0105200), the National Natural Science Foundation of China (NSFC) under Grants No. 61722211, 61773362, 61872338, and 61906180, the Lenovo-CAS Joint Lab Youth Scientist Project, and the Foundation and Frontier Research Key Program of Chongqing Science and Technology Commission (No. cstc2017jcyjBX0059), the Tencent AI Lab Rhino-Bird Focused Research Program (No. JR202033).


\section{Introduction}
The ESA/Rosetta mission was the first mission  to approach and follow a comet for a large part of its orbit, as well as to land on its nucleus. The Rosetta orbiter and the Philae lander were equipped with a series of instruments to study in-situ the gas and dust coma surrounding comet 67P/Churyumov-Gerasimenko (hereafter 67P) and the properties of its nucleus. The mission provided unprecedented insights into the structure, composition, and activity evolution of 67P, and comets in general \citep{Taylor2017}.

In order to complement the mission, an ambitious ground-based campaign was put together, with dozens of telescopes across the globe and in space performing observations of the comet at the same time as the Rosetta mission \citep{Snodgrass2017}. Those observation probed the large scale coma of 67P and provided context to the measurements performed by the Rosetta mission. More importantly, they provide a link with ground-based observations of a large number of comets that have been performed over the last century. As part of this campaign, observations were performed with the Multi Unit Spectroscopic Explorer (MUSE) instrument at the Very Large Telescopes (VLT). The observations were performed in March 2016, after the comet's perihelion passage, when the comet was moving away from the Sun. In this paper, we present the results of the MUSE observations of comet 67P.

\section{Observations and data reduction}

MUSE is an integral field unit spectrograph mounted on the UT4 telescope of the VLT, in Chile \citep{Bacon2010}, nominally covering the 480-930 nm range. In Wide Field mode, which was used for our observations, MUSE has a field of view (FoV) of 1\arcmin$\times$1\arcmin$ $ covered without gaps. The instrument has a platescale of 0.2\arcsec /pix and a spectral resolving power between 1770 (at 480 nm) and 3590 (at 930 nm).
Observations of comet 67P were performed between 2016 March 3 and March 7. In total, 38 datacubes were obtained over 5 nights. For 4 of those cubes, the comet was not visible, or very badly centred, so that we discarded the data. The sky during those 5 nights was either clear or photometric, but the seeing was variable, between 0.6 and 2.3\arcsec. For all observations, we used an exposure time of 600s, and the position angle of the instrument was set to $0\degree$ (we did not apply any rotation between the exposures). At the time of the observations, even though the comet was still active, and thus extended, it did not fill the entire MUSE FoV. Because of that, no dedicated sky observations were performed. The observing circumstances are presented in Table \ref{TableObs}.

\begin{table*}
\caption{\label{TableObs} Observing circumstances of the 67P MUSE campaign}
\centering
\begin{tabular}{lcccccc}
\hline\hline
Date (UT) & N & Exposure time (s) & r (au) & $\Delta$ (au) & Airmass & Phase Angle (\degree)\\
\hline
2016-03-03 03h27-06h08 & 12 & 600 & 2.49 & 1.52 & 1.2-1.6 & 6 \\
2016-03-04 05h34 & 1 & 600 & 2.50 & 1.53 & 1.2 & 6 \\
2016-03-05 04h50-06h29 & 8 & 600 & 2.51 & 1.53 & 1.2-1.3 & 5\\
2016-03-06 03h25-04h00 & 4 & 600 & 2.52 & 1.54 & 1.4-1.5 & 5 \\
2016-03-07 04h39-07h38 & 13 & 600 & 2.52 & 1.54 & 1.3-1.6 & 4\\
\hline
\end{tabular}
\end{table*}

The data reduction was performed using the ESO pipeline \citep{Weilbacher2016}, with the sky estimated from regions near the edge of the cubes that are free from the comet contamination. In addition to the cube reconstruction, bias subtraction, flatfield correction, wavelength calibration, and sky subtraction, the ESO MUSE pipeline also corrects for the telluric absorption and flux-calibrates the science data, using a standard star observed the same night as the science observations.

Even though the sky was estimated directly on the science cubes and subtracted by the pipeline, while examining the reduced cubes we noticed relatively strong sky residuals. In order to reduce those, we then used the ZAP (Zurich Atmosphere Purge) software \citep{Soto2016}. ZAP is a Principal Component Analysis-based software designed to perform sky subtraction on IFU data. The use of ZAP allows us to remove most of the sky residuals left after the pipeline reduction. As will be discussed later, we also performed a full data reduction without applying any sky subtraction, in order to search for forbidden oxygen emission lines in the coma of 67P.

As mentioned before, the ESO pipeline performs a correction of the telluric absorption using a standard star observed on the same night as the science data. However, the standard stars are primarily  used for flux calibration and are not optimal for telluric correction. Also, they are not always observed just after the science observations nor at the same airmass. Because of that, in the reduced cubes, there are residuals of the strong O$_2$ telluric band around 760 nm if the telluric correction is performed with the pipeline. For all the analysis focused on the 2D structure of the dust coma, and the detection of forbidden oxygen lines, this does not impact the quality of our measurements, as those residuals are relatively constant over the field and we extracted the cubes over wavelength ranges that mostly avoid that region of the spectrum. However, in Section \ref{Ref_Spectrum}, we focus on the spectrum of the dust in the coma of 67P and residual from the telluric correction could impact the quality of the spectrum presented. For the spectra presented in that section, we thus used the Molecfit software \citep{Smette2015,Kausch2015} to perform a better telluric correction. Molecfit has proven to be very useful to provide an accurate telluric correction even when dedicated observations of a standard star were not available. We applied Molecfit directly on the extracted 1D spectra. The fit of the telluric features was performed for each extracted spectrum individually.

\section{Analysis}

\subsection{A reference dust spectrum}
\label{Ref_Spectrum}

   \begin{figure*}[h!]
   \centering
   \includegraphics{Dust_Spectrum.png}
   \caption{Top: Spectrum of 67P dust coma (black), together with the solar spectrum used to compute the reflectance spectrum (red). The blue dotted line is the spectrum extracted over a 10\arcsec$ $ aperture.  The solar spectrum has been shifted arbitrarily over the y-axis for better visibility. Bottom: Relative reflectance spectrum of 67P (black), compared to the one measured with X-Shooter in November 2014 \citep[red;][]{Snodgrass2016}.}
    \label{Dust_Spect}
    \end{figure*}

We first study the spectrum of comet 67P over the optical range. To do so, we extract all the data cubes over a 5 pixels radius (1\arcsec) aperture around the comet optocenter. We chose such a small aperture to focus on the part of the coma where the signal is the strongest, to avoid having strong sky residuals in the presented spectrum. After the extraction, the spectra are corrected for the telluric features using the Molecfit software, as explained above. All 34 spectra are then median-combined to produce a high quality spectrum of the dust in the coma of 67P. This spectrum is shown in the upper part of Fig. \ref{Dust_Spect}. We do not see any emission band in the spectrum, because of the large heliocentric distance and low activity level of the comet at the time of the observations. In the same part of the Figure, we show a reference solar spectrum obtained using the SOLar SPECtrometer (SOLSPEC) instrument of the SOLAR payload on board the International Space Station (ISS) \citep{Meftah2018}, re-sampled to match the sampling of the 67P spectrum. Given that the aperture we chose is very small, the spectrum shown in Fig. \ref{Dust_Spect} probably contains a non-negligible nucleus contribution. For comparison purposes, we also extract the dust spectrum over a 10\arcsec$ $ aperture, in which the dust coma dominates (blue dashed curve in the top part of Fig. \ref{Dust_Spect}). The two spectra match well, except for a small difference around 880 nm. This region is severely affected by telluric absorption and has strong sky emissions, that might explain this difference. In general, we can thus say that the spectrum extracted over the 1\arcsec$ $ aperture is representative of the coma dust and we will discuss this spectrum only in the following text.

In the bottom part of Fig. \ref{Dust_Spect}, we divide the comet spectrum by the Solar reference spectrum to compute the relative reflectance of the dust in the coma of 67P. The reflectance spectrum is normalised at 600 nm. Overlaid to the reflectance spectrum measured with MUSE, is the reflectance measured using the X-Shooter spectrograph in November 2014 over the same wavelength range \citep{Snodgrass2016}. The reflectance spectra measured in November 2014 and March 2016 are consistent with each other, indicating that the dust reflectance as measured from the ground is similar at large heliocentric distance pre- and post-perihelion. We do not see any sign of absorption bands in the optical spectrum of 67P. We measure a spectral slope for the reflectance, or dust reddening, of 10\%$/100$ nm in the 480-900 nm interval. We also notice that the slope becomes shallower at longer wavelengths. We measure a slope of 13 \%$/100$ nm in the 500-700 nm interval but only of 5\%$/100$ nm in the 700-900 nm interval. This is fully consistent with what is reported from X-Shooter observations performed in 2014, from which values between 10\%$/100$ nm and 20\%$/100$ nm were reported in the 550 to 1000 nm interval, the shallower values corresponding to the red end of the wavelength range.

Our measurements are also consistent with in-situ measurements of the dust from the ESA/Rosetta mission. \cite{Bertini2017} report a reddening in the interval 376-744 nm ranging between 11 and 14\%$/100$ nm from measurements with the OSIRIS cameras. Similarly, \cite{LaForgia2019} report slopes measured in the inner coma ranging from 12 to 16\%$/100$ nm between 480 and 649 nm. Those slopes are similar to those measured for the nucleus of 67P, both in-situ and from ground-based observations. \cite{Fornassier2015} report average slopes of 11$-$16\%$/100$ nm over 250$-$1000 nm with shallower slopes toward longer wavelengths, from observation of 67P's nucleus with the Rosetta/OSIRIS camera. From ground-based observations, \cite{Tubiana2011} measure a slope for the nucleus of $12\pm1$\%$/100$ nm over 430$-$850 nm, slightly shallower in the 500$-$850 nm range. Finally, the dust reddening measured in the coma of 67P is consistent with what is measured usually in the coma of active comets, typically between 0 and 20\%$/100$ nm, and with a shallower slope towards the near-IR \citep{Solontoi2012,Jewitt1986}.

\subsection{Dust coma morphology and activity}
\label{Morphology}

For each cube in our dataset, we extract maps over the bandpasses of the V, R, and I Johnson-Cousins filters. All maps in the same bandpass are re-centred and then co-added. The resulting co-added maps are displayed in the top part of Fig. \ref{DustMaps}. The dust coma morphology is the same for all three bandpasses. It is asymmetrical, probably due to the presence of dust jets. To investigate further the presence of jets, we divide the dust maps by an azimuthal median profile \citep{Samarasinha2013}. Enhanced maps are displayed on the bottom part of Fig. \ref{DustMaps}. On those we can clearly see two jets, which are most likely the cause of the apparent asymmetry of the dust coma. The first jet is located close to the anti-sunward direction, towards the South-West. The second jet is located about 90\degree$ $ away, towards the South-East. Finally, we see a faint feature towards the sunward direction. Those jets are consistent with what was observed on previous passages and what is reported by \cite{Knight2017} and \cite{Snodgrass2017} from observations at the same epoch. They are also consistent with the modelling of the pole orientation and active region location done by \cite{Vincent2013} prior to the Rosetta mission. In addition to the two jets mentioned above, we see an enhancement towards the North-West. This corresponds to the dust trail, that was reported to be at least two degrees long at that epoch \citep{Snodgrass2017,Boehnhardt2016,Knight2017}. We do not see changes of the coma morphology over the 5 nights during which we have observations, nor over a single night.

\begin{figure*}[h!]
\centering
\includegraphics[width=5.2cm,height=5.2cm]{MapV_Tag.png}
\includegraphics[width=5.2cm,height=5.2cm]{MapR_Tag.png}
\includegraphics[width=5.2cm,height=5.2cm]{MapI_Tag.png}
\includegraphics[width=5.2cm,height=5.2cm]{MapsV_En.png}
\includegraphics[width=5.2cm,height=5.2cm]{MapR_En.png}
\includegraphics[width=5.2cm,height=5.2cm]{MapI_En.png}
\caption{Top row: combined maps of 67P in V, R, and I bands. The maps are centred on the comet and the FoV is 1\arcmin$\times$1\arcmin. Bottom row: V, R, and I maps enhanced by dividing by an azimuthal median profile \citep{Samarasinha2013}.}
\label{DustMaps}
\end{figure*}

To constrain the comet activity at that time, we compute the Af$\rho$ value, which is a proxy for dust production, as defined by \cite{A'Hearn1984}. We compute the $\mathrm{Af\rho}$ for the V, R, and I bands over a 2500 km physical aperture. We use a 2500 km aperture instead of the more commonly used 10000 km aperture because the comet signal drops significantly at 10000 km. We obtain values of 65$\pm$4 cm, 75$\pm$4 cm, and 82$\pm$4 cm in the V, R, and I bands respectively. Those values have not been corrected for the phase angle effect and are consistent with those reported by \cite{Boehnhardt2016} at the same epoch, and comparable (even though slightly lower) to those reported by \cite{Knight2017} once the phase angle effect is taken into account. Using the Af$\rho$ values in the different bands, we can also compute the dust reflectivity gradient defined as \citep{Jewitt1986}:

\[ \mathrm{S' = \frac{Af\rho_2- Af\rho_1}{Af\rho_2+ Af\rho_1} \frac{20000}{\lambda_2-\lambda_1}}\]

where $\mathrm{Af\rho_x}$ and the corresponding wavelength $\lambda_x$ are expressed in cm and nm, respectively. For the I-V combination, this gives a reflectivity gradient of 7$\pm$1 \%$/100$ nm, while it is of 11$\pm$2 \%$/100$ nm for the I-R and 7$\pm$1 \%$/100$ nm for the R-V combination. This confirms the trend outlined in section \ref{Ref_Spectrum} of higher reflectivity gradient at lower wavelengths. To check for trends with the aperture size, we computed the Af$\rho$ in apertures of 5000 km, 7500 km, and 10000 km. Within the error bars, we do not see significant changes in the $\mathrm{Af\rho}$ or reflectivity gradient values with the aperture size in our data, as expected for a steady-state coma.

\subsection{Gas detection}
\label{Gas}
Among the features usually observed in the optical spectrum of comets, the oxygen forbidden line at 630 nm is one of the brightest. There are three forbidden oxygen lines that can be detected in the coma of comets at optical wavelengths: the green line at 557.7 nm, and the red doublet at 630 and 636.4 nm. Those lines are emitted by the decay of oxygen atoms in a metastable ($^1$D) or ($^1$S) state. Exited atomic oxygen is mainly produced through the photo-dissociation of $\mathrm{H_2O}$, $\mathrm{CO}$, $\mathrm{CO_2}$, or even $\mathrm{O_2}$ (see, e.g. \cite{Cessateur2016}). Observing forbidden oxygen lines then represents an opportunity to constrain the production of those potential parent species from optical observations. However, the cometary forbidden oxygen lines are often blended with the equivalent atmospheric lines, unless observed at high spectral resolution.

The spectral resolution of MUSE does not allow us to resolve the telluric and cometary lines for the geocentric velocity of comet 67P at the time of our observations. In theory, the sky subtraction should subtract the atmospheric contribution, allowing us to recover the cometary signal. However, the noise introduced by the subtraction of the very strong atmospheric features prevents us from detecting any cometary forbidden oxygen emission lines. At the time of our observations the comet was at 2.5 au post-perihelion and was only weakly active, as confirmed by the fact that we do not detect any emission lines in the spectrum presented in Fig. \ref{Dust_Spect}. If the oxygen lines are present, they are thus very faint and masked by the noise introduced by the sky subtraction.

We therefore attempt to detect the forbidden oxygen lines using another method. We reduce the full dataset without performing any sky subtraction. For each spaxel of each cube, we then subtract the continuum sky contribution and the dust contribution underlying the oxygen forbidden lines by defining continuum region on both sides of the red doublet and the green oxygen lines, fitting a line through those regions and subtracting it. Finally, we extract the datacubes over a very narrow wavelength range centred on the wavelength of the three oxygen lines. In low/medium resolution spectra, forbidden oxygen lines can be contaminated by $\mathrm{NH_2}$ or $\mathrm{C_2}$ lines, but this is highly unlikely in the case of those observations since the comet is weakly active and no emission lines are detected in Fig. \ref{Dust_Spect}. The result is maps of the flux contained in the forbidden lines at 557.7, 630, and 636.4 nm over the whole MUSE FoV. Since the atmospheric and cometary lines cannot be resolved by MUSE, the maps contain the sum of the atmospheric and potential cometary contribution. All maps are re-centred so that the optocenter of the comet is placed at the same position. We average all 34 maps for each line, performing a 3-$\sigma$ clipping. The result in shown in Fig. \ref{OMAPS}

In those maps, we expect the atmospheric contribution to be relatively uniform over the 1\arcmin$\times$1\arcmin$ $ FoV of MUSE. The comet contribution, on the other side, would be concentrated around the optocenter of the comet, with a limited spatial extension since the states responsible for the emission lines are metastable states. In Fig \ref{OMAPS}, we can see that for the 557.7 and 636.4 nm lines, the maps are mostly uniform over the whole FoV. Small variations are observed but they most likely come from non-optimal corrections of the detector-to-detector effects. We thus conclude that the signal in those maps comes from the atmospheric oxygen lines. In the 630 nm map, however, in addition to inhomogeneities similar to those seen in the other two maps, we have a clear over-density located around the position of the comet optocenter. This indicates that, in addition to the atmospheric line, we detect signal from the cometary forbidden oxygen line at 630 nm. It is not surprising that we only detect the 630 nm line, since it is the strongest of the three oxygen lines. The signal we measure in the comet aperture for the 630 nm map is barely 3-$\sigma$ above the sky background variation (measured from four adjacent sky apertures, see later). Since the ratio between the 630 nm and 636.4 nm maps is expected to be 3 (as it is the case for the background in the two maps, see Fig. \ref{OMAPS}), it is then consistent with the fact that we do not detect significant cometary signal in the 636.4 nm map. We note that the [OI] cometary signal is enhanced towards the South, similarly to what is observed for the dust. However, we do not believe the signal is due to residuals from the dust continuum. Indeed, the technique used for the dust subtraction is the same of all 3 lines and only the 630 nm map (with the strongest oxygen feature) shows an over-density around the optocenter. Control maps built by using a wavelength range adjacent to the 630 nm line do not show a similar signal.

\begin{figure*}[h!]

\includegraphics[width=18cm]{Plot_OMAPS.png}
\caption{Combined maps of the sky+comet flux at the wavelengths of the 3 forbidden oxygen lines at 557.7 nm (left), 630 nm (middle), and 636.4 nm (right). The FoV is 52\arcsec$\times$52\arcsec and is oriented North up and East left and the black arrow in the centre points to the position of the comet optocenter. The unit of the color bars is $10^{-20}$ erg/s/cm$^2$.}
\label{OMAPS}
\end{figure*}

The cometary signal we detect in the 630 nm map is faint, but we use this signal to estimate the water production rate of the comet at that time. The 630 nm line comes from the decay of atomic oxygen in ($^1$D) state. As mentioned before, oxygen in that state is mainly produced by the photo-dissociation of $\mathrm{H_2O}$, $\mathrm{CO}$, $\mathrm{CO_2}$, or $\mathrm{O_2}$. For comets at 1 au from the Sun, $\mathrm{H_2O}$ photodissociation is the dominant source for producing metastable oxygen, but at distance above 2.5 au other molecules such as CO and $\mathrm{CO_2}$ start to contribute \citep{Decock2013,McKay2015}. Long-term measurements of the production rate of all 4 species of interest derived from  the ROSINA instrument onboard Rosetta are presented by \cite{Laeuter2020} and \cite{Combi2020}. At the time of our observations, the CO and $\mathrm{O_2}$ production rates are more than a factor 10 and 100 lower than that of $\mathrm{H_2O}$, respectively. We can thus assume that those species will contribute little to the production of metastable oxygen. The $\mathrm{CO_2}$ production rate is only about a factor 2 to 4 lower than the $\mathrm{H_2O}$ production rate. The emission rate for O($^1$D) production from $\mathrm{CO_2}$ is 1.5 times higher than that of $\mathrm{H_2O}$ \citep{Bhardwaj2012}. Taking this into account, a significant part (up to half) of O($^1$D) atoms could be produced by the photo-dissociation of $\mathrm{CO_2}$. For the measurements presented below, we will assume that water is the main source for the production of metastable oxygen, so that the water production rate we derive might be overestimated.

In order to derive water production rates, we follow the procedure described in \cite{Schultz1992,Morgenthaler2001}. We consider a photo-chemical model including the following three reactions,
\begin{equation}
\mathrm{H_2O + h\nu \longrightarrow H_2 + O(^1D)}
\label{equ1}
\end{equation}
\begin{equation}
\mathrm{H_2O + h\nu \longrightarrow H + OH}
\label{equ2}
\end{equation}
\begin{equation}
\mathrm{OH + h\nu \longrightarrow H + O(^1D)}
\label{equ3}
\end{equation}
the water production rate is given by:

\[ \mathrm{Q(H_2O) = \frac{Q([OI])}{BR1+(BR2)(BR3)} }\]
where BR1 and BR2 are the branching ratios for reactions \ref{equ1} and \ref{equ2} \citep{Huebner1992} (equal to 0.05 and 0.855 for the quiet Sun, respectively) and BR3 is the branching for reaction \ref{equ3} from \cite{Morgenthaler2001} (equal to 0.094). The production rate of oxygen atoms in ($^1$D) state, Q([OI]), can be obtained from:
\[ \mathrm{Q([OI]) = \frac{4}{3}4\pi\Delta^2I_{630}AC }\]
where $\mathrm{\Delta}$ is the comet geocentric distance (in cm), $\mathrm{I_{630}}$ is the intensity of the [OI] emission (in $\mathrm{photons/s/cm^2}$), and AC is the aperture correction factor correcting for the [OI] emission not encompassed in the aperture.

We measure the comet flux in the 630 nm map using a 25 pixels (5\arcsec$ $) radius circular aperture centred on the comet optocentre. In order to measure and subtract the sky flux, we measure the flux in four 5\arcsec$ $ apertures located at four different positions with respect to the comet optocenter ((+50 pix,+50 pix), (+50 pix,-50 pix), (-50 pix,+50 pix), (-50 pix,-50 pix)), computed the average, and subtracted it from the flux measured in the comet aperture. Given the faintness of the cometary emission, the main uncertainty in the determination of the water production rate comes from the measurement of the sky. To estimate this uncertainty, we used the upper and lower values measured in the individual sky apertures. The resulting intensity of the comet 630 nm line is $1.8\pm1.0 \times 10^{-3}  \mathrm{photons/cm^2/s}$. The aperture we used encompasses all the visible emission from the comet, so that we set the aperture correction factor to 1. This results in a water production rate of $1.5\pm0.7 \times 10^{26} \mathrm{molec./s}$. This represents the actual water production rate only if all the oxygen atoms in ($^1$D) state are produced by the photo-dissociation of water.

No other measurement of the comet water production rate was reported from ground-based observations at a similar epoch because the comet was too faint and weakly active. In-situ measurements made with the ROSINA DFMS instruments determined a water production rate around $2.3\times 10^{26}$ molec./s at the same heliocentric distance \citep{Hansen2016}. Similar values are reported by \cite{Combi2020,Laeuter2020} using the same instrument. \cite{Biver2019} report a water production rate of $8.5\pm2.5 \times 10^{25} \mathrm{molec./s}$ from measurements with the MIRO instrument. Given that $\mathrm{CO_2}$ could potentially contribute up to half of the 630 nm [OI] line brightness, the actual water production rate we measure could be as low as $0.7 \times 10^{26} \mathrm{molec./s}$, which is in very good agreement with the \cite{Biver2019} measurement. Given the difference in technique, scale of the observations and models used, our measurement is close to those reported using Rosetta instruments in particular the MIRO instrument.

\section{Discussion and Conclusions}

\begin{figure*}[h!]
\centering
\includegraphics[width=17.2cm]{Rosetta_Plot.png}
\caption{Images of the nucleus of 67P taken by the OSIRIS instrument on-board Rosetta during the 19/02/2016 outburst and at the same epoch as the MUSE observations presented here. Credit: ESA/Rosetta/MPS for OSIRIS Team MPS/UPD/LAM/IAA/SSO/INTA/UPM/DASP/IDA.}
\label{Plot_Rosetta}
\end{figure*}

Comet 67P was observed with the MUSE IFU over 5 nights, between March 3 and 7, 2016 when the comet was at 2.5 au from the Sun and 1.5 au from the Earth. The spectrum of the dust coma of 67P presented here is of high quality. It has a good signal to noise ratio given the faintness and the distance of the comet at the time of our observations. Due to the use of the Molecfit software it is little affected by telluric features at near-IR wavelengths, which is not the case of most comet spectra in the 800-900 nm range. It matches very well other ground-based observations of 67P, as well as in-situ measurements from the Rosetta spacecraft. It is also representative of what is usually observed for active comets in general. Finally, it does not contain gas emission features above the noise level (some [OI] signal is present, merged with the sky line; see below). For all those reasons, this spectrum could be used in the future as a 'template' dust spectrum to help perform the subtraction of the dust-reflected continuum for observations of other comets with the MUSE IFU. This spectrum could replace the observation of solar analogs since so far very few good solar analogs have been observed with MUSE. The dust reddening can vary from one comet to another, so that this spectrum would need to be corrected from any slope difference between the target and 67P. Nonetheless, this dust-only spectrum of 67P represents a good tool to help for the analysis of future comet observations with MUSE.

The morphology of the coma, the dust activity, and the dust reflectivity gradient as measured with MUSE are consistent with other measurements performed with ground-based telescopes at the same epoch as well as with in-situ measurements from the Rosetta mission.

About two weeks before our observations, on February 2 2016, an outburst was detected by several instruments on-board the Rosetta orbiter \citep{2016Grun}. A similar but less intense event was also detected in NAVCAM images on 1 March 2016. Since the outburst(s) happened shortly before our observations, we could have expected to detect some indications of such event, either in terms of activity level or coma morphology, in our observations. However, we do not detect any sign of the outburst in the MUSE data presented here. In Fig. \ref{Plot_Rosetta}, we show images obtained with the OSIRIS WAC on-board Rosetta. In the left part, gas jets are easily detected during the February 19 outburst. The right side shows the comet at the same epoch as our observations. No gas jets are visible. This is consistent with the comet being back to a quieter state and the fact that we do not detect the outburst in the MUSE observations.

Through careful reprocessing and separating comet and sky signals spatially (rather than by resolving them spectroscopically), we detect the 630 nm forbidden oxygen line and derive a water production rate of $1.5\pm0.7 \times 10^{26} \mathrm{molec./s}$ in the coma of 67P, if all [OI] atoms in ($^1$D) state are produced by water. This value is consistent with the Rosetta measurements \citep{Hansen2016,Combi2020,Laeuter2020,Biver2019}. It is the only measurement of the water production rate of 67P at such large heliocentric distance from remote observations. Ground-based detection of water (or water products) in the coma of 67P only happened close to perihelion \citep{Snodgrass2017}, so that this point is precious to compare measurements and models from the Rosetta spacecraft to ground-based measurements over a larger portion of the comet orbit. In general, detecting water at such low production rate for comets at 2.5 au from the Sun is extremely difficult (see e.g. discussion on water detection in \cite{Snodgrass2017b}).
MUSE had been shown in the past to have a huge potential to study species parentage in the coma of brighter comets \citep{Opitom2019}.
Our present work demonstrate the efficiency of MUSE to detect low level of water production around distant solar system bodies. This opens up future opportunities, at a time where an intensive search for evidence of water ice in the main asteroid belt and elsewhere in the Solar System is being performed.

\begin{acknowledgements}
 Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme 096.C-0160(A). Datasets of the ESA/Rosetta OSIRIS instrument have been downloaded from the ESA Planetary Science Archive. The authors acknowledge the OSIRIS Principal Investigator H. Sierks (MPS, Goettingen, Germany) and the OSIRIS Team for providing images and related datasets and the ESA Rosetta Project for enabling the science of the mission.

-
\end{acknowledgements}

\bibliographystyle{aa}


\section{The conjecture}
We use $\lambda_i(G)$ to denote the $i$th largest eigenvalue of the adjacency matrix of a simple graph $G$ on $n$ vertices,
for $i=1,2,\cdots, n$. By the Perron-Frobenius Theorem, $\lambda_1$ is always positive (unless $G$ has no edges) and
$|\lambda_i|\le \lambda_1$ for all $i\ge 2$. Let $\lambda = \max_{2\le i\le n} |\lambda_i| =\max\{|\lambda_2|, |\lambda_n|\}$,
that is, $\lambda$ is the second largest absolute eigenvalue. For a $d$-regular graph, it is well known that $\lambda_1=d$.

Let $c(G)$ denote the number of components of a graph $G$.
The {\bf toughness} $t(G)$ of a connected graph $G$ is defined as $t(G)=\min\{\frac{|S|}{c(G-S)}\}$,
where the minimum is taken over all proper subsets $S\subset V(G)$ such that $c(G-S)>1$.
A graph $G$ is {\bf $k$-tough} if $t(G)\ge k$.
Graph toughness was introduced by Chv\'atal \cite{Chva73} in 1973 and is closely related to
many graph properties, including connectivity, Hamiltonicity, pancyclicity, factors, spanning trees, etc.

The study of toughness from eigenvalues was initiated by Alon~\cite{Alon95} who showed that
for any connected $d$-regular graph $G$, $t(G)>\frac{1}{3}(\frac{d^2}{d\lambda+\lambda^2}-1)$, through which,
Alon was able to show that for every $t$ and $g$ there are $t$-tough graphs of girth strictly greater than $g$.
This strengthened a result of Bauer, Van den Heuvel and Schmeichel \cite{BaVS95} who showed the same for $g=3$, and
thus disproved in a strong sense a conjecture of Chv\'atal \cite{Chva73} that there exists a constant $t_0$ such that
every $t_0$-tough graph is pancyclic.

Brouwer~\cite{Brou95} independently discovered a better bound and showed that $t(G)>\frac{d}{\lambda}-2$ for a connected
$d$-regular graph $G$. He mentioned in~\cite{Brou95} that the bound might be able to be improved a little to $\frac{d}{\lambda}-1$,
and then posed the exact conjecture as an open problem in \cite{Brou96}.
Some partial results have been provided in \cite{CiGu16}. However, no substantial progress has been made for
more than two decades. Most recently, the author improved the bound to $t(G) > \frac{d}{\lambda} - \sqrt{2}$ in \cite{Gu21}.

\begin{conjecture}[Brouwer~\cite{Brou95, Brou96}]
\label{tconj}
For any connected $d$-regular graph $G$, $t(G)\ge \frac{d}{\lambda}-1$.
\end{conjecture}

It is mentioned by Brouwer~\cite{Brou95} that there are infinitely many graphs $G$ such that $t(G)\le d/\lambda$
with equality in many cases (for example, strongly regular graphs constructed in  \cite{Brou95} and \cite{CiWo14},
and many Kneser graphs \cite{POHBWWC}). Brouwer~\cite{Brou95} also pointed out that Conjecture~\ref{tconj}, if true, is tight.
To see sharpness, one may notice that the toughness can be arbitrarily close to $0$, while $d/\lambda\ge 1$ holds for all
$d$-regular graphs. Cioab\u{a} and Wong~\cite{CiWo14} briefly described an explicit construction of such extremal graphs.

In this paper, we confirm Conjecture~\ref{tconj}.
\begin{theorem}
\label{tofg}
Let $G$ be a connected $d$-regular graph. Then $\displaystyle t(G)\ge \frac{d}{\lambda}-1.$
\end{theorem}

\section{The proof}
Our main tool is the expander mixing lemma. A $d$-regular graph on $n$ vertices with the second largest absolute eigenvalue
at most $\lambda$ is called an $(n, d, \lambda)$-graph.
It is well known that an $(n, d, \lambda)$-graph  for which $\lambda = \Theta(\sqrt{d})$ is a very good pseudo-random graph
behaving, in many aspects, like a truly random graph $G(n, p)$. The quantitative definition of pseudo-random graphs was introduced
by Thomason~\cite{Thom85, Thom87} who defined jumbled graphs.

The celebrated expander mixing lemma is usually attributed to Alon and Chung~\cite{AlCh88}.
This idea appeared earlier with a different form in the PhD thesis~\cite{Haem79} of Haemers.
We present this lemma  in Theorem~\ref{expmix} (see \cite{KrSu06} for a complete proof).
A variation can be found in \cite[Chapter 9]{AlSp16} by Alon and Spencer.
We refer readers to the informative survey \cite{KrSu06} by Krivelevich and Sudakov for more about the expander mixing lemma
and pseudo-random graphs. As mentioned in \cite{KrSu06}, the expander mixing lemma is a truly remarkable result, connecting
edge distribution and graph spectrum, and providing a very good quantitative handle for the uniformity of edge distribution based
on graph eigenvalues.

For every two subsets $A$ and $B$ of $V(G)$, let $e(A, B)$ denote the number of edges with
one end in $A$ and the other one in $B$ (edges with both ends in $A\cap B$ are counted twice).
We use $e(A)$ to denote the number of edges with both ends in $A$, and thus $e(A, A) = 2 e(A)$.
\begin{theorem}[Expander Mixing Lemma]
\label{expmix}
Let $G$ be an $(n, d, \lambda)$-graph. Then for every two subsets $A$ and $B$ of $V(G)$,
\begin{equation}
\label{genmix}
\left| e(A, B) - \frac{d}{n} |A| |B| \right| \le \lambda \sqrt{|A| |B| \left(1-\frac{|A|}{n}\right)\left(1-\frac{|B|}{n}\right)}.
\end{equation}
In particular,
\begin{equation}
\label{oneset}
\left| e(A) - \frac{d}{2n} |A|^2 \right| \le \frac{\lambda}{2} |A| \left(1-\frac{|A|}{n}\right).
\end{equation}
\end{theorem}

The following lemma has been proved and used by Brouwer and Haemers~\cite{BrHa05}. For the sake of completeness,
we include a short proof here.
\begin{lemma}[\hspace{-1sp}\cite{BrHa05}]
\label{indexlem}
Let $x_1,\cdots, x_c$ be positive integers such that $\sum_{i=1}^c x_i \le 2c-1$. Then for every integer $\ell$ with
$0\le \ell \le \sum_{i=1}^c x_i$, there exists an $I\subset \{1,\cdots, c\}$ such that $\sum_{i\in I} x_i =\ell$.
\end{lemma}
\begin{proof}
The proof provided by Brouwer and Haemers~\cite{BrHa05} goes by induction on $c$. The case $c=1$ is trivial.
Let $c\ge 2$ and assume that $x_1\le \cdots \le x_c$. Suppose it is true for $c-1$ integers.
Let $\ell' = \ell$ if $\ell \le c-1$ and $\ell' = \ell - x_c$ otherwise.
Since $c-1\le \sum_{i=1}^{c-1} x_i = \sum_{i=1}^c x_i - x_c\le 2(c-1)-1$, by inductive hypothesis,
there exists an $I'\subset \{1,\cdots, c-1\}$ such that $\sum_{i\in I'} x_i =\ell'$.
Let $I = I'$ or $I = I'\cup \{c\}$, as desired.
\end{proof}

Now we are ready to present the proof of the main theorem.
\begin{proof}[\bf Proof of Theorem~\ref{tofg}]
Suppose to the contrary that
\begin{equation}
\label{contr}
t(G) < \frac{d}{\lambda}-1.
\end{equation}
By definition, suppose $S$ is a subset of $V(G)$ such that $\frac{|S|}{c(G-S)}=t(G)$. Let $B = V(G-S)$.
Denote $|V(G)|=n$, $c(G-S)=c$ and $t(G)=t$. Then $|S| = tc$ and so $|B| = n - tc$.

First we show that
\begin{equation}
\label{cupb}
c\le \frac{\lambda n}{d + \lambda}.
\end{equation}
In fact, (\ref{cupb}) can be obtained from the well-known Hoffman ratio bound. Here we give a direct proof.
To see this, let $U$ be a set of vertices that consists of exactly one vertex from each component of $G-S$.
Then $|U|=c$ and $e(U)=0$. By (\ref{oneset}), $\left| e(U) - \frac{d}{2n} |U|^2 \right| \le \frac{\lambda}{2} |U| \left(1-\frac{|U|}{n}\right)$,
and thus $\frac{d}{n} |U| \le \lambda \left(1-\frac{|U|}{n}\right)$, which implies that $c=|U|\le  \frac{\lambda n}{d + \lambda}$.\\

If $|B| \le \frac{2\lambda n}{d+ \lambda}$, then $tc = n - |B| \ge \frac{(d - \lambda)n}{d+\lambda}$. Together with (\ref{cupb}), we have
$$t= (n - |B|)/c \ge \frac{(d - \lambda)n}{d+\lambda} \cdot \frac{d + \lambda}{\lambda n}
= \frac{d - \lambda}{\lambda} = \frac{d}{\lambda} -1,$$
contrary to (\ref{contr}). Thus, we may assume that
\begin{equation}
\label{B-lowb}
|B| > \frac{2\lambda n}{d+ \lambda},
\end{equation}
By (\ref{cupb}) and (\ref{B-lowb}), we have $|B| > 2c$, that is
\begin{equation}
\label{bvsc}
|B|\ge 2c+1.
\end{equation}

Let $V_1,\cdots, V_c$ denote the vertex sets of the $c$ components of $G-S$.
Without loss of generality, we may assume that $|V_1|\le \cdots\le |V_c|$.

\par\medskip
\noindent
{\bf Claim 1:} $\sum_{i=1}^{c-1} |V_i|\ge c$.
\\Proof of Claim 1: Otherwise, $\sum_{i=1}^{c-1} |V_i| = c-1$, that is, each $V_i$ is a single vertex for $i=1,\cdots, c-1$.
Let $A= \cup_{i=1}^{c-1} V_i$ and thus $|A|=c-1$ and $e(A, B) = 0$. By (\ref{genmix}),
$$\frac{d}{n} |A| |B| \le \lambda \sqrt{|A| |B| \left(1-\frac{|A|}{n}\right)\left(1-\frac{|B|}{n}\right)}
< \lambda \sqrt{|A| |B|\left(1-\frac{|B|}{n}\right)},$$
and so
$$\left(\frac{d}{n} |A| |B| \right)^2 < \lambda^2 |A| |B|\left(1-\frac{|B|}{n}\right),$$
which implies that
$$\frac{d^2 |A| |B|}{\lambda^2 n} < n-|B| = tc.$$

Thus
$$t > \frac{|A|}{c}\cdot \frac{d^2}{\lambda^2}\cdot \frac{|B|}{n} = \frac{c-1}{c}\cdot \frac{d^2}{\lambda^2}\cdot \frac{|B|}{n}
= \left(1-\frac{1}{c} \right) \cdot \frac{d^2}{\lambda^2}\cdot \frac{|B|}{n}\ge \frac{d^2}{2\lambda^2}\cdot \frac{|B|}{n}.$$

By (\ref{B-lowb}), $\frac{|B|}{n} > \frac{2\lambda}{d+ \lambda}$ and we have
$$t> \frac{d^2}{2\lambda^2}\cdot \frac{2\lambda}{d+ \lambda} = \frac{d^2}{\lambda(d+ \lambda)}
=\frac{d^2 - \lambda^2 + \lambda^2}{\lambda(d+ \lambda)}
= \frac{d^2 -\lambda^2}{\lambda(d+ \lambda)} + \frac{\lambda^2}{\lambda(d+ \lambda)}
> \frac{d}{\lambda} -1,$$ contrary to (\ref{contr}). This completes the proof of Claim 1.

\par\medskip
\noindent
{\bf Claim 2:} $V_1,\cdots, V_c$ can be partitioned into two disjoint sets $X$ and $Y$ such that $e(X, Y)=0$, $|X|\ge c$ and  $|Y|\ge c$.
\\Proof of Claim 2:
If $|V_c|\ge c$, then  let $X = \cup_{i=1}^{c-1} V_i$ and $Y = B-X = V_c$. Then $|Y|\ge c$ and by Claim 1, $|X|\ge c$.
Thus we may assume that $|V_c| \le c-1$. We also notice that $|V_c|\ge 3$, for otherwise, $|B|=\sum_{i=1}^c |V_i| \le 2c$, contradicting
(\ref{bvsc}). Thus $3\le |V_c| \le c-1$. Let $\ell = c - |V_c|$. Then $1\le \ell \le c-3$.

If $\sum_{i=1}^{c-1} |V_i|\le 2(c-1) -1 = 2c-3$, then by Lemma~\ref{indexlem}, there exists an $I\subset \{1,\cdots, c-1\}$ such that
$\sum_{i\in I} |V_i| =\ell$. Let $X = \cup_{i\in I} V_i \cup V_c$ and $Y = B-X$. Clearly $|X| = \ell + |V_c| =c$ and $|Y| = |B| -|X|\ge c+1$.

If $\sum_{i=1}^{c-1} |V_i| > 2(c-1) -1 = 2c-3$, then let $V'_i$ be a nonempty subset of $V_i$ for each $i=1,2,\cdots c-1$ such that
$\sum_{i=1}^{c-1} |V'_i| = 2c-3$. We can easily do it by removing vertices from some $V_i$'s one by one until the total number of
remaining vertices is exactly $2c-3$ (keep at least one vertex in each subset). By Lemma~\ref{indexlem}, there exists an
$I\subset \{1,\cdots, c-1\}$ such that $\sum_{i\in I} |V'_i| =\ell$. Then $\sum_{i\not\in I, i< c} |V'_i| = (2c-3)-\ell \ge (2c-3) - (c-3) = c$.
Let $X = \cup_{i\in I} V_i \cup V_c$ and $Y = B-X = \cup_{i\not\in I, i< c} V_i$.
Clearly $|X| \ge \ell + |V_c| =c$ and $|Y| = \sum_{i\not\in I, i< c} |V_i|\ge \sum_{i\not\in I, i< c} |V'_i| \ge c$.
This completes the proof of Claim 2.

\par\medskip
By Claim 2, $V_1,\cdots, V_c$ can be partitioned into two disjoint sets $X$ and $Y$ such that $e(X, Y)=0$, $|X|\ge c$ and  $|Y|\ge c$.

Since $e(X, Y)=0$, by (\ref{genmix}),
$\frac{d}{n} |X| |Y| \le \lambda \sqrt{|X| |Y| \left(1-\frac{|X|}{n}\right)\left(1-\frac{|Y|}{n}\right)},$
and so $$ d^2 |X|^2 |Y|^2 \le \lambda^2 |X| |Y| \left(n- |X|\right)\left(n-|Y|\right),$$
which implies that
\begin{equation}
\label{xynoeg}
|X| |Y| \le \frac{\lambda^2}{d^2} (n - |X|)(n - |Y|).
\end{equation}

Without loss of generality, we may assume that $|X|\le |Y|$. By (\ref{xynoeg}), we have
$$ |X|^2 \le |X|\cdot |Y| \le \frac{\lambda^2}{d^2}(n-|X|)(n-|Y|)\le \frac{\lambda^2}{d^2}(n-|X|)^2,$$
that is $$|X| \le \frac{\lambda}{d} (n - |X|),$$
and hence
\begin{equation}
\label{xxbound}
|X| \le \frac{\lambda n}{d+\lambda}.
\end{equation}

Also, since $|Y| = n - |S| - |X|$, by (\ref{xynoeg}) again, we have
$$ |X| (n - |S| - |X|) = |X|\cdot |Y| \le \frac{\lambda^2}{d^2}(n - |X|)(n - |Y|) = \frac{\lambda^2}{d^2}(n - |X|)(|S| + |X|),$$
which implies that
\begin{equation}
\label{xnxn}
|X|n \le \left(\frac{\lambda^2}{d^2}(n - |X|) + |X| \right) \left(|S| + |X| \right)
= \left(\frac{\lambda^2}{d^2}n + \frac{d^2 - \lambda^2}{d^2}|X| \right) \left(|S| + |X| \right).
\end{equation}
By (\ref{xxbound}), we have
$$\frac{d^2 - \lambda^2}{d^2}|X| \le \frac{d^2 - \lambda^2}{d^2}\cdot \frac{\lambda n}{d+\lambda} = \frac{d\lambda - \lambda^2}{d^2}n,$$
plugging in (\ref{xnxn}),
$$|X|n \le \left(\frac{\lambda^2}{d^2}n + \frac{d\lambda - \lambda^2}{d^2}n \right) \left(|S| + |X| \right) = \frac{\lambda n}{d} \left(|S| + |X| \right),$$
and we have
$$ |X| \le \frac{\lambda}{d} \left(|S| + |X| \right).$$
Hence
\begin{equation*}
tc = |S| \ge \left(\frac{d}{\lambda} -1 \right) |X|,
\end{equation*}
implying that
$$t \ge \left(\frac{d}{\lambda} -1 \right)\cdot \frac{|X|}{c} \ge \frac{d}{\lambda} -1,$$
completing the proof of the theorem.
\end{proof}