\section{Introduction}

\begin{figure}[t]
	\begin{center}
		\includegraphics{pull_figure}
	\end{center}
	\caption{More robust models tolerate larger noise levels. Thus, To assess the robustness of VQA models, we apply noise at a controllable level to the main question. We sort a dataset of basic questions based on their similarity to the main question and append three basic questions at a time as noise. The robustness is then measured by }\def\rscore/{$R_{\text{score}}$/ as the deterioration in accuracy over a given testing dataset.}
	\label{fig:pull_figure}
\end{figure}

Visual Question Answering (VQA) is one of the most challenging computer vision tasks in which an algorithm is given a natural language question about an image and tasked with producing an answer for that image-question pair. Recently, various VQA models \cite{4,9,31,37,41,42,57,58,59} have been proposed to tackle this problem, and their main performance measure is accuracy. In general, any model must have in some sense few aspects of quality relevant to the problem it solves. Accuracy, robustness, stubbornness, and myopia are some examples of such qualities for VQA models that were studied by the community \cite{29}. Obviously, some qualities are more important than others and an overall evaluation of a model should rationalize the interactions of these qualities if possible. In this work, we are interested only in the robustness of VQA models to small \emph{noise} or \emph{perturbation} to the input question, dubbed main question (MQ).

\subsection{Motivation}
The idea of analyzing model robustness as well as training robust models is already a rapidly growing research topic for deep learning models applied to images, }\def\eg/{{\em e.g.}/ adversarial attacks in image classification \cite{61,63}. However, most VQA research study other quality aspects that test certain signs of intelligence, such as stubbornness, compositionality, and myopia. Robustness is a concept that despite its importance is heavily misinterpreted and commonly neglected. To the best of our knowledge, we are the first to propose a measure for the robustness of VQA models.

\subsection{Assumptions}
Robustness to some noises (}\def\eg/{{\em e.g.}/, motion blur in images) is more important than others (}\def\eg/{{\em e.g.}/, salt and pepper) because they are more common. Additionally, some noises have controllable strength, which allows for more sophisticated robustness measures (}\def\eg/{{\em e.g.}/, area under the performance curve versus noise level). In this work, we are establishing a framework for analyzing robustness against a type of \emph{controllable additive noise} that is relevant to VQA but is not common. Namely, plain-text concatenation of the MQ with semantically similar questions, dubbed basic questions (BQs); refer to Figure \ref{fig:pull_figure} for a visual illustration. The approach, of considering those BQs as noise, and similar concepts have merit to them and they are studied in-depth in psychology under ``deductive reasoning in human thinking'' despite being a controversial topic \cite{74}. For instance, a person quizzed on the ``color of the bus'' in an image of a car will write down the color of the car while acknowledging the discrepancy. Comparably, if they were asked multiple similar questions, they might be unhinged but answer correctly nevertheless.

\subsection{Robustness Evaluation Framework}
Figure \ref{fig:pull_figure}, depicts the noise generation part of our robustness evaluation. Given a plain-text question (}\def\ie/{{\em i.e.}/, the MQ) and a plain-text basic questions dataset (BQD), we start by ranking the BQs in BQD by their similarity to the MQ using some text similarity ranking method; we formulate this ranking as }\def\lasso/{\textit{LASSO}/ optimization problem. Then, depending on the required level of noise, we take the top, }\def\eg/{{\em e.g.}/ $n=3$, ranked BQs and we append them one after the other. The concatenation of these BQs with MQ is the generated noisy question. To measure the robustness of any VQA model, the accuracy over a testing dataset with and without adding this generated noise is compared using our proposed }\def\rscore/{$R_{\text{score}}$/.

This framework is posed such that extrapolating to other types of noise is straightforward. For example, \emph{question-words-shuffling noise} can be controlled by the hamming distance (}\def\eg/{{\em e.g.}/, the larger the hamming distance the higher level of noise) whereas \emph{question-rephrasing} is a noise that is not trivially controllable. We apply the noise of choice to the questions in a testing dataset, then evaluate the deterioration in accuracy using }\def\rscore/{$R_{\text{score}}$/. A more comprehensive robustness evaluation of a VQA model should employ multiple types of noise, possibly ones that could be jointly applied on the image and the MQ at the same time.

\subsection{Contributions}
\begin{enumerate}[i]
	\item We propose a novel framework to measure the robustness of VQA models and test it on six different models.
	\item We propose a new text-based similarity ranking method, }\def\ie/{{\em i.e.}/ }\def\lasso/{\textit{LASSO}/, and compare it against BLEU-1, BLEU-2, BLEU-3, BLEU-4 \cite{49}, ROUGE \cite{68}, CIDEr \cite{67} and METEOR \cite{69}.
	\item We provide the similarity ranking of two large-scale BQDs to test the robustness of VQA models: General Basic Question Dataset (GBQD) and Yes/No Basic Question Dataset (YNBQD).
\end{enumerate}

\section{Related Work}
There is a big effort from the community to tackle the problem of VQA \cite{26,33,54,65,66}. It is a multidisciplinary task that involves natural language progressing (NLP), computer vision and machine learning.

\subsection{Sentence Evaluation Metrics}
Sentence evaluation metrics have been widely used in different areas such as \emph{text summarization} and \emph{machine translation}. In our work, we exploit these metrics to measure the similarity between the MQs and the BQs. BLEU \cite{49} is one of the most popular metrics in machine translation that is based on precision. However, its effectiveness was questioned by \cite{70,71}. METEOR \cite{69}, however, is based on precision and recall. In addition, ROUGE \cite{68} is a popular recall-based metric in the text summarization. It tends to reward longer sentences with higher recall. Moreover, a consensus-based metric, CIDEr \cite{67}, rewards a sentence for being similar to the majority of descriptions written by a human. In our experiments, we take all of the aforementioned metrics to rank BQDs, and our experimental results show that our proposed }\def\lasso/{\textit{LASSO}/ ranking method achieves better ranking performance.

\subsection{Sentence Embedding}
There exists many method that analyze the relationship between words, phrases and sentences by mapping text into some latent vector spaces \cite{11,43,47}. It was shown that if two phrases share the same context in the corpus, their embedded vectors will be close to each other in the latent space. Skip-thoughts \cite{43} can map text to embedding space using an encoder-decoder architecture of recurrent neural networks (RNNs). The encoder is an RNN with gated recurrent unit (GRU) activation \cite{10}. While the decoder is an RNN with a conditional GRU. We use this model because it performs well on embedding long sentences.

\subsection{Attention Mechanism in VQA}
Attention-based VQA models can attend to local image regions related to the query question \cite{13,15,35,42}. In the pooling step of \cite{42}, the authors use an image-attention mechanism to help determine the relevance of updated question representations to the original. As far as we know, no work has tried to apply the mechanism of language attention to VQA models before \cite{41}. They propose a mechanism of co-attention that jointly performs language and image attention.

\subsection{Multiple Modality Fusion Approaches in VQA}
The VQA task considers features from both the question and image and can be viewed as a multimodal feature fusion task. The authors of \cite{26,57,58,59} have tried to focus on modeling the interactions between two different embedding spaces. In \cite{26}, the authors demonstrate the success of the bilinear interaction between two embedding spaces in deep learning for fine-grained classification and multimodal language modeling. Multimodal Compact Bilinear (MCB) pooling \cite{58} exploits an outer product between visual and textual embedding. Moreover, Multimodal Low-rank Bilinear (MLB) pooling \cite{59} uses a tensor to parametrize the full bilinear interactions between question and image spaces. MCB and MLB are efficiently generalized in \cite{57}.

\subsection{Robustness of Neural Network}
In \cite{61}, the authors analyze the robustness of deep models by adding some perturbations into images and observe how the prediction result is affected. In \cite{72}, the authors mention that VQA models can produce different answers by using slight variations of a query question implying that VQA models do not actually comprehend the asked question.

\section{Methodology}
Our goal is to evaluate the robustness of VQA models over the MQs in the testing set of the popular VQA dataset \cite{4}. For each (image, MQ) pair in the testing set, we concatenate to MQ the most three similar questions (}\def\ie/{{\em i.e.}/, BQs) to MQ, obtained by a certain similarity ranking method, over a large-scale basic questions dataset (BQD) as in Figure \ref{fig:pull_figure}. Finally, we compute the accuracy before and after adding this noise and compare them using }\def\rscore/{$R_{\text{score}}$/. In this section, we will go through the details of this procedure.

\subsection{Datasets Preparation}
Our BQD is the combination of only unique questions from the training and validation sets of ``real images'' in \emph{the VQA dataset} \cite{4}, which is a total of $186027$ questions. However, because there will be questions in the testing set of $244302$ questions that is also in this BQD, we will exclude them only when executing the similarity ranking.

We will need feature representations of our questions in both the testing set and the BQD. We will use the Skip-thought vector \cite{43} which exploits an RNN encoder with GRU \cite{10} activation to map any English sentence to a feature vector $\mathbf{v} \in \mathbf{R}^{4800}$.

\subsection{Similarity Ranking}
We can project the problem of finding similar questions to a given MQ using this BQD and a similarity ranking method. We will consider two categories of question similarity ranking methods; \emph{direct similarity} and \emph{sparse combination}.

\subsubsection{Direct Similarity} We will need a text similarity (or dissimilarity) scoring method, that takes two sentences and return a similarity (or dissimilarity) score of how ``close'' (or ``far'') are those two semantically from each other. By computing the score of a given MQ to all BQs in the BQD we can directly rank the questions from the most to the least similar. Here, we will use seven scoring methods; BLEU-1, BLEU-2, BLEU-3, BLEU-4, ROUGE, CIDEr and METEOR.

\subsubsection{Sparse Combination} We will need the feature representations of all the questions and a distance measure to determine a sparse combination of BQs needed to represent a given MQ. This can be modelled as follows:
\begin{equation}
	\min_{\mathbf{x}}~\text{distance}\left(\mathbf{A}\mathbf{x}, \mathbf{b} \right)+\lambda \left \| \mathbf{x} \right \|_{1},
	\label{eq:lasso}
\end{equation}
where $\mathbf{b}$ is the feature vector of the MQ, $\mathbf{A}$ is the matrix of feature vectors of all questions in BQD as its columns, and $\lambda$ is a trade-off parameter that controls the sparsity of the solution. Note that all these vectors are normalized to a unit $\ell_2$-norm. To cast this as a }\def\lasso/{\textit{LASSO}/ optimization problem
$$\text{distance}(\mathbf{q}, \mathbf{p}) = \frac{1}{2}\|\mathbf{q} - \mathbf{p}\|_2^2$$
Different distance metric could give different results especially if they make sense for the latent embedding space. The $i$\textsuperscript{th} element of the solution $\mathbf{x}$ is the similarity score of the BQ encoded in the $i$\textsuperscript{th} column of $\mathbf{A}$. We reiterate the importance of making sure that $\mathbf{b}$ is not one of the columns of $\mathbf{A}$. Otherwise, and because we are encouraging sparsity, the ranking will give all other BQs a zero similarity score.

\subsection{Robustness Evaluation}
In the VQA dataset \cite{4}, a predicted answer can be considered partially correct if it matched the answers of less than three human annotators. The overall accuracy is:
\begin{equation}
	\text{Accuracy}_\text{VQA}=\frac{1}{N}\sum_{i=1}^{N}\min\left \{ \frac{\sum_{t\in T_{i}}\mathbb{I}[a_{i}=t]}{3},1 \right \}
\end{equation}
where $\mathbb{I}[\cdot]$ is the indicator function, $N$ is the total number of examples, $a_{i}$ is the predicted answer, and $T_{i}$ is human annotators answer set of the $i^{th}$ image-question pair.

\subsubsection{The }\def\rscore/{$R_{\text{score}}$/ of VQA Models}
First, we measure the accuracy of the model on the clean testing set and denote it as $\text{Acc}_\text{VQA}$. Then, we append the top ranked $k$ BQs to each of the MQs and recompute the accuracy of the model on this noisy dataset and we call it $\text{Acc}_\text{BQD}$. Finally, we compute the absolute difference $\text{Acc}_\text{diff} = \left|\text{Acc}_\text{VQA}-\text{Acc}_\text{BQD}\right|$ and report the robustness
\begin{equation}
	R_\text{score} = \text{clamp}_{0}^{1}\left(\frac{\sqrt{m}-\sqrt{{\text{Acc}_\text{diff}}}}{\sqrt{m}-\sqrt{t}}\right)
\end{equation}
Here, we apply this min-max clipping of the score
$$\text{clamp}_{a}^{b}(x) = \max\left(a,{\min\left(b,x\right)}\right)$$
where $0 \leq t < m \leq 1$. The parameters $t$ and $m$ are the tolerance and maximum robustness limit, respectively. In fact, }\def\rscore/{$R_{\text{score}}$/ is designed to decrease smoothly from $1$ to $0$ as $\text{Acc}_\text{diff}$ moves from $t$ to $m$ and remains constant outside this range. The rate of change of this transition is exponentially decreasing from exponential to sub-linear in the range $[t, m]$. The reasoning behind this is that we want the score to be sensitive if $\text{Acc}_\text{diff}$ is small, but not before $t$, and less sensitive if it is large, but not after $m$. Note that }\def\rscore/{$R_{\text{score}}$/ will be penalized even if $\text{Acc}_\text{VQA}$ was less than $\text{Acc}_\text{BQD}$, }\def\ie/{{\em i.e.}/  when we observe an improvement in the accuracy after adding the noise.

\subsection{Basic Questions Datasets}
The size of the BQD has a great impact on the noise generation method. Basically, the more questions we have, the more chance it has to contain similar questions to any given MQ. On the other hand, solving Eq \ref{eq:lasso} for each question in the VQA testing set of $244302$ questions for the $186027$ unique questions in both training and validation set can become computationally expensive. Luckily, we need to do it at most once. We set $\lambda = {10}^{-6}$ and keep the top-$k$ BQs for each MQ produced by solving Eq \ref{eq:lasso}, where $k = 21$. The resultant BQD will have $244302$ instance of this format:
$$\{\text{Image},~MQ,~21~(BQ + \text{corresponding~ similarity~score})\}$$
We keep the same output 10 answers, annotated by different workers through AMT (Amazon Mechanical Turk), for each instance as specified by the VQA dataset \cite{4} in the open-ended and multiple-choice (18 choice) tasks. Approximately $90\%$ of the answers only have a single word and $98\%$ of the answers have no more than three words. In addition, these datasets contain $81434$ testing images from the MS COCO dataset \cite{51}.

Because most of the VQA models have the highest accuracy performance in answering yes/no questions, we have collected two BQDs; Yes/No Basic Question Dataset (YNBQD) and General Basic Question Dataset (GBQD). Now, in practice, our proposed BQDs can be used directly to test the robustness of VQA models without running the }\def\lasso/{\textit{LASSO}/ ranking method again. In our experiments, we will compare using }\def\lasso/{\textit{LASSO}/ (}\def\ie/{{\em i.e.}/, \emph{sparse combination}) to using other seven \emph{direct similarity} metrics in building BQDs.

\section{Experiments and Analysis}
To analyze our proposed framework, we will perform our experiments on six VQA models; LQI denoting LSTM Q+I \cite{4}, HAV denoting HieCoAtt (Alt, VGG19) and HAR denoting HieCoAtt (Alt, Resnet200) \cite{41}, MU denoting MUTAN without Attention and MUA denoting MUTAN with Attention \cite{57}, and MLB denoting MLB with Attention \cite{4}. On top of that, we will limit ourselves to the open-ended task on the test-dev partition from the 2017\textsuperscript{th} VQA Challenge \cite{4}, denoted here as dev, unless otherwise specified like using the test-std partition, denoted std.

\subsection{Comparing Similarity Ranking Methods}
In addition to building GBQD and YNBQD using }\def\lasso/{\textit{LASSO}/, we also generate them using the other similarity metrics (}\def\ie/{{\em i.e.}/, BLEU-1, BLEU-2, BLEU-3, BLEU-4, ROUGE, CIDEr and METEOR) as described in the Methodology section. Then, for each MQ, we split the top-21 ranked BQs into seven partition each of which contains three consecutive BQs to control the noise level. More concretely, partition 1, which is (BQ1, BQ2, BQ3), has smaller noise than partition 7, which is (BQ19, BQ20, BQ21), as illustrated in Figure \ref{fig:pull_figure}, with partition 0 being the empty partition. Finally, some VQA models limit the number of words in the MQ, which might lead to trimming the appended noise (}\def\ie/{{\em i.e.}/ the BQs in the partition).

\begin{figure}
	\begin{subfigure}[t]{0.49\linewidth}
		\includegraphics{LASSO_GBQD}
		\caption{LASSO\_GBQD}
		\label{fig:lasso_gbqd}
	\end{subfigure}
	\hfill
	\begin{subfigure}[t]{0.49\linewidth}
		\includegraphics{LASSO_YNBQD}
		\caption{LASSO\_YNBQD}
		\label{fig:lasso_ynbqd}
	\end{subfigure}
	\caption{Compares the accuracy of six VQA models with increasing noise levels from both GBQD and YNBQD. The x-axis is the partition index with 0 meaning MQ without noise. We can see a monotonous trend as the noise increases.}
	\label{fig:lasso}
\end{figure}

Figure \ref{fig:lasso} shows the accuracy of the six VQA models with the noise coming from GBQD and YNBQD generated by our proposed }\def\lasso/{\textit{LASSO}/ similarity ranking method. We can see that the accuracy is decreasing as the partition index increases, }\def\ie/{{\em i.e.}/ confirming our assumption of treating these BQs as noise and that the noise level is increasing as the similarity score of the BQs decreases. However, we repeated the same experiments on GBQD generated by the other similarity ranking methods in Figure \ref{fig:others_gbqd} and we couldn't see this trend anymore. The plots became less monotonous and acting randomly as we move from partition 1 to partition 7. Also, we observe a big drop in accuracy starting from partition 1 rendering these similarity measures ineffective in this context. In spite of that, various work \cite{1,19,27,28,48} still use them for sentence similarity evaluation because of their simple implementation. This fact signifies the need to develop better similarity metrics, which in turn can be directly evaluated using this flexible framework.

\subsection{Evaluating }\def\rscore/{$R_{\text{score}}$/ with }\def\lasso/{\textit{LASSO}/ Ranking}
}\def\lasso/{\textit{LASSO}/ Ranking prevails over the other seven similarity metrics. In Table \ref{tbl:lasso_examples}, We present two qualitative }\def\lasso/{\textit{LASSO}/ ranking examples to showcase the limitations. In Tables \ref{tbl:lasso_gbqd} and \ref{tbl:lasso_ynbqd}, we compare the six VQA models on GBQD and YNBQD and report their }\def\rscore/{$R_{\text{score}}$/ on only partition 1 with $t=5\times{10}^{-4}$ and $m=2\times{10}^{-1}$. As we can see, generally speaking, attention-based models (}\def\ie/{{\em i.e.}/ HAV, HAR, MUA and MLB) are more robust than non-attentive ones (}\def\ie/{{\em i.e.}/ LQI and MU) with HieCoAtt being the most robust. However, MU (non-attentive) is more robust on YNBQD than MUA (the attentive version), while it is not the case on GNBQD. This shows how }\def\rscore/{$R_{\text{score}}$/ is strongly tied to the domain of the BQD and that we can only compare models under the same BQD. Note that, we studied only one property of the BQD domain, which is the type of the answer (}\def\ie/{{\em i.e.}/, yes-no vs. general).

\begin{figure*}
	\begin{subfigure}[t]{0.49\linewidth}
		\includegraphics{ROUGE_and_BLEU_GBQD}
		\caption{ROUGE, BLEU-4, BLEU-3, BLEU-2 and BLEU-1}
	\end{subfigure}
	\hfill
	\begin{subfigure}[t]{0.49\linewidth}
		\includegraphics{METEOR_and_CIDEr_GBQD}
		\caption{CIDEr~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(c) METEOR}
	\end{subfigure}
	\caption{Compares the accuracy of six VQA models with increasing noise levels from GBQD generated by BLEU-1, BLEU-2, BLEU-3, BLEU-4, ROUGE, CIDEr and METEOR. We could not observe a trend as the noise level increase like in Figure \ref{fig:lasso}.}
	\label{fig:others_gbqd}
\end{figure*}

\begin{table*}
	\begin{subtable}{0.49\linewidth}
		\begin{tabular}{l|l|l}
			Rank & MQ: How old is the car?                 & Score  \\
			\hline
			BQ01 & How old is the truck?                   & 0.2952 \\
			BQ02 & How old is this car?                    & 0.2401 \\
			BQ03 & How old is the vehicle?                 & 0.1416 \\
			BQ04 & What number is the car?                 & 0.1195 \\
			BQ05 & What color is the car?                  & 0.0933 \\
			BQ06 & How old is the bedroom?                 & 0.0630 \\
			BQ07 & What year is the car?                   & 0.0630 \\
			BQ08 & Where is the old car?                   & 0.0368 \\
			BQ09 & How old is the seat?                    & 0.0330 \\
			BQ10 & How old is the cart?                    & 0.0320 \\
			BQ11 & What make is the blue car?              & 0.0281 \\
			BQ12 & How old is the golden retriever?        & 0.0280 \\
			BQ13 & What is beneath the car?                & 0.0239 \\
			BQ14 & Is the car behind him a police car?~~~~ & 0.0223 \\
			BQ15 & How old is the pilot?                   & 0.0198 \\
			BQ16 & How old are you?                        & 0.0172 \\
			BQ17 & How old is the laptop?                  & 0.0159 \\
			BQ18 & How old is the television?              & 0.0157 \\
			BQ19 & What make is the main car?              & 0.0149 \\
			BQ20 & What type and model is the car?         & 0.0148 \\
			BQ21 & What is lifting the car?                & 0.0145 \\
		\end{tabular}
	\end{subtable}
	\begin{subtable}{0.49\linewidth}
		\begin{tabular}{l|l|l}
			Rank & MQ: What is the cat sitting on?         & Score  \\
			\hline
			BQ01 & Where is the cat sitting on?            & 0.2812 \\
			BQ02 & What is this cat sitting on?            & 0.1079 \\
			BQ03 & What is cat sitting on?                 & 0.0548 \\
			BQ04 & What is the cat on the left sitting on? & 0.0526 \\
			BQ05 & What is the giraffe sitting on?         & 0.0499 \\
			BQ06 & What is the cat sitting in the car?     & 0.0469 \\
			BQ07 & That is the black cat sitting on?       & 0.0462 \\
			BQ08 & What is the front cat sitting on?       & 0.0423 \\
			BQ09 & What is the cat perched on?             & 0.0414 \\
			BQ10 & What's the cat sitting on?              & 0.0408 \\
			BQ11 & What is the cat leaning on?             & 0.0370 \\
			BQ12 & What object is the cat sitting on?      & 0.0351 \\
			BQ13 & What is the doll sitting on?            & 0.0289 \\
			BQ14 & How is the cat standing?                & 0.0226 \\
			BQ15 & What is the cat setting on?             & 0.0220 \\
			BQ16 & What is the cat walking on?             & 0.0217 \\
			BQ17 & What is the iPhone sitting on?          & 0.0210 \\
			BQ18 & What is the cat napping on?             & 0.0210 \\
			BQ19 & What is the dog sitting at?             & 0.0201 \\
			BQ20 & What is the birds sitting on?           & 0.0183 \\
			BQ21 & What is the sitting on?                 & 0.0182 \\
		\end{tabular}
	\end{subtable}
	\caption{Demonstrates the limitations of }\def\lasso/{\textit{LASSO}/ ranking with two examples and their corresponding similarity scores despite performing well in our experiments. From the table, we readily see how the misspellings and sentence re-phrasings can affect the quality of the generated BQD. This can be attributed to the limitation of the encoding and with better semantic encoders and appropriate distance metrics this framework should improve out-of-the-box without any significant changes. In addition, one can utilize the feature vectors to filter out the unique questions to form the BQD instead of relying on raw string comparisons.}
	\label{tbl:lasso_examples}
\end{table*}

\begin{table*}
	\begin{subtable}{0.49\linewidth}
		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = 0.30  & Other & Num   & Y/N   & All            & Diff           \\ [0.5ex]
			\hline
			dev: partition 1 & 37.78 & 34.93 & 68.20 & \textbf{49.96} & \textbf{10.20} \\
			dev: partition 2 & 37.29 & 35.03 & 65.62 & \textbf{48.67} & \textbf{11.49} \\
			dev: partition 3 & 34.81 & 34.39 & 62.85 & \textbf{46.27} & \textbf{13.89} \\
			dev: partition 4 & 34.25 & 34.29 & 63.60 & \textbf{46.30} & \textbf{13.86} \\
			dev: partition 5 & 33.89 & 34.66 & 64.19 & \textbf{46.41} & \textbf{13.75} \\
			dev: partition 6 & 33.15 & 34.68 & 64.59 & \textbf{46.22} & \textbf{13.94} \\
			dev: partition 7 & 32.80 & 33.99 & 63.59 & \textbf{45.57} & \textbf{14.59} \\
			\hline
			std: partition 1 & 38.24 & 34.54 & 67.55 & \textbf{49.93} & \textbf{10.52} \\
			\hline
			dev: partition 0 & 47.16 & 37.32 & 81.45 & \textbf{60.16} & -              \\
			std: partition 0 & 47.57 & 36.75 & 81.56 & \textbf{60.45} & -              \\
			\hline
		\end{tabular}
		\caption{MUTAN without Attention model evaluation results.}

		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = \textbf{0.48} & Other & Num   & Y/N   & All            & Diff           \\ [0.5ex]
			\hline
			dev: partition 1         & 44.44 & 37.53 & 71.11 & \textbf{54.63} & \textbf{5.85}  \\
			dev: partition 2         & 42.62 & 36.68 & 68.67 & \textbf{52.67} & \textbf{7.81}  \\
			dev: partition 3         & 41.60 & 35.59 & 66.28 & \textbf{51.08} & \textbf{9.4}   \\
			dev: partition 4         & 41.09 & 35.71 & 67.49 & \textbf{51.34} & \textbf{9.14}  \\
			dev: partition 5         & 39.83 & 35.55 & 65.72 & \textbf{49.99} & \textbf{10.49} \\
			dev: partition 6         & 39.60 & 35.99 & 66.56 & \textbf{50.27} & \textbf{10.21} \\
			dev: partition 7         & 38.33 & 35.47 & 64.89 & \textbf{48.92} & \textbf{11.56} \\
			\hline
			std: partition 1         & 44.77 & 36.08 & 70.67 & \textbf{54.54} & \textbf{5.78}  \\
			\hline
			dev: partition 0         & 49.14 & 38.35 & 79.63 & \textbf{60.48} & -              \\
			std: partition 0         & 49.15 & 36.52 & 79.45 & \textbf{60.32} & -              \\
			\hline
		\end{tabular}
		\caption{HieCoAtt (Alt, VGG19) model evaluation results.}

		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = 0.36  & Other & Num   & Y/N   & All            & Diff           \\ [0.5ex]
			\hline
			dev: partition 1 & 49.31 & 34.62 & 72.21 & \textbf{57.12} & \textbf{8.67}  \\
			dev: partition 2 & 48.53 & 34.84 & 70.30 & \textbf{55.98} & \textbf{9.81}  \\
			dev: partition 3 & 48.01 & 33.95 & 69.15 & \textbf{55.16} & \textbf{10.63} \\
			dev: partition 4 & 47.20 & 34.02 & 69.31 & \textbf{54.84} & \textbf{10.95} \\
			dev: partition 5 & 45.85 & 34.07 & 68.95 & \textbf{54.05} & \textbf{11.74} \\
			dev: partition 6 & 44.61 & 34.30 & 68.59 & \textbf{53.34} & \textbf{12.45} \\
			dev: partition 7 & 44.71 & 33.84 & 67.76 & \textbf{52.99} & \textbf{12.80} \\
			\hline
			std: partition 1 & 49.07 & 34.13 & 71.96 & \textbf{56.95} & \textbf{8.73}  \\
			\hline
			dev: partition 0 & 57.01 & 37.51 & 83.54 & \textbf{65.79} & -              \\
			std: partition 0 & 56.60 & 36.63 & 83.68 & \textbf{65.68} & -              \\
			\hline
		\end{tabular}
		\caption{MLB with Attention model evaluation results.}
	\end{subtable}
	\begin{subtable}{0.49\linewidth}
		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = 0.34  & Other & Num   & Y/N   & All            & Diff           \\ [0.5ex]
			\hline
			dev: partition 1 & 51.51 & 35.62 & 68.72 & \textbf{56.85} & \textbf{9.13}  \\
			dev: partition 2 & 49.86 & 34.43 & 66.18 & \textbf{54.88} & \textbf{11.10} \\
			dev: partition 3 & 49.15 & 34.50 & 64.85 & \textbf{54.00} & \textbf{11.98} \\
			dev: partition 4 & 47.96 & 34.26 & 64.72 & \textbf{53.35} & \textbf{12.63} \\
			dev: partition 5 & 47.20 & 33.93 & 64.53 & \textbf{52.88} & \textbf{13.10} \\
			dev: partition 6 & 46.48 & 33.90 & 64.37 & \textbf{52.46} & \textbf{13.52} \\
			dev: partition 7 & 46.88 & 33.13 & 64.10 & \textbf{52.46} & \textbf{13.52} \\
			\hline
			std: partition 1 & 51.34 & 35.22 & 68.32 & \textbf{56.66} & \textbf{9.11}  \\
			\hline
			dev: partition 0 & 56.73 & 38.35 & 84.11 & \textbf{65.98} & -              \\
			std: partition 0 & 56.29 & 37.47 & 84.04 & \textbf{65.77} & -              \\
			\hline
		\end{tabular}
		\caption{MUTAN with Attention model evaluation results.}

		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = 0.45  & Other & Num   & Y/N   & All            & Diff           \\ [0.5ex]
			\hline
			dev: partition 1 & 46.51 & 36.33 & 70.41 & \textbf{55.22} & \textbf{6.59}  \\
			dev: partition 2 & 45.19 & 36.78 & 67.27 & \textbf{53.34} & \textbf{8.47}  \\
			dev: partition 3 & 43.87 & 36.28 & 65.29 & \textbf{51.84} & \textbf{9.97}  \\
			dev: partition 4 & 43.41 & 36.25 & 65.94 & \textbf{51.88} & \textbf{9.93}  \\
			dev: partition 5 & 42.02 & 35.89 & 66.09 & \textbf{51.23} & \textbf{10.58} \\
			dev: partition 6 & 42.03 & 36.40 & 65.66 & \textbf{51.12} & \textbf{10.69} \\
			dev: partition 7 & 40.68 & 36.08 & 63.49 & \textbf{49.54} & \textbf{12.27} \\
			\hline
			std: partition 1 & 46.77 & 35.22 & 69.66 & \textbf{55.00} & \textbf{7.06}  \\
			\hline
			dev: partition 0 & 51.77 & 38.65 & 79.70 & \textbf{61.81} & -              \\
			std: partition 0 & 51.95 & 38.22 & 79.95 & \textbf{62.06} & -              \\
			\hline
		\end{tabular}
		\caption{HieCoAtt (Alt, Resnet200) model evaluation results.}

		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = 0.19  & Other & Num   & Y/N   & All            & Diff           \\ [0.5ex]
			\hline
			dev: partition 1 & 29.24 & 33.77 & 65.14 & \textbf{44.47} & \textbf{13.55} \\
			dev: partition 2 & 28.02 & 32.73 & 62.68 & \textbf{42.75} & \textbf{15.27} \\
			dev: partition 3 & 26.32 & 33.10 & 60.22 & \textbf{40.97} & \textbf{17.05} \\
			dev: partition 4 & 25.27 & 31.70 & 61.56 & \textbf{40.86} & \textbf{17.16} \\
			dev: partition 5 & 24.73 & 32.63 & 61.55 & \textbf{40.70} & \textbf{17.32} \\
			dev: partition 6 & 23.90 & 32.14 & 61.42 & \textbf{40.19} & \textbf{17.83} \\
			dev: partition 7 & 22.74 & 31.36 & 60.60 & \textbf{39.21} & \textbf{18.81} \\
			\hline
			std: partition 1 & 29.68 & 33.76 & 65.09 & \textbf{44.70} & \textbf{13.48} \\
			\hline
			dev: partition 0 & 43.40 & 36.46 & 80.87 & \textbf{58.02} & -              \\
			std: partition 0 & 43.90 & 36.67 & 80.38 & \textbf{58.18} & -              \\
			\hline
		\end{tabular}
		\caption{LSTM Q+I model evaluation results.}
	\end{subtable}
	\caption{Compares the accuracy and }\def\rscore/{$R_{\text{score}}$/ of six VQA models with increasing noise levels from GBQD generated by }\def\lasso/{\textit{LASSO}/ evaluated on dev and std. The results are split by the question type; Numerical (Num), Yes/No (Y/N), or Other. There is a multitude of things to consider here in order to give a well-informed interpretation of these scores. One of which is the way $\text{Acc}_\text{diff}$ is defined in Eq \ref{eq:lasso}. The absolute value can be replaced with a ReLU (}\def\ie/{{\em i.e.}/, $max(x, 0)$) but this is a design decision which is widely accepted in the literature of adversarial attacks. As we can see, the accuracy of all tested models decrease as they are evaluated on noisier partitions ranked by }\def\lasso/{\textit{LASSO}/ (i.e., partition 7 has a higher noise level than partition 1).}
	\label{tbl:lasso_gbqd}
\end{table*}

\begin{table*}
	\begin{subtable}{0.49\linewidth}
		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = 0.30  & Other & Num   & Y/N   & All            & Diff           \\ [0.5ex]
			\hline
			dev: partition 1 & 33.98 & 33.50 & 73.22 & \textbf{49.96} & \textbf{10.13} \\
			dev: partition 2 & 32.44 & 34.47 & 72.22 & \textbf{48.67} & \textbf{11.18} \\
			dev: partition 3 & 32.65 & 33.60 & 71.76 & \textbf{46.27} & \textbf{11.36} \\
			dev: partition 4 & 32.77 & 33.79 & 71.14 & \textbf{46.30} & \textbf{11.53} \\
			dev: partition 5 & 32.46 & 33.51 & 70.90 & \textbf{46.41} & \textbf{11.81} \\
			dev: partition 6 & 33.02 & 33.18 & 69.88 & \textbf{46.22} & \textbf{12.00} \\
			dev: partition 7 & 32.73 & 33.28 & 69.74 & \textbf{45.57} & \textbf{12.18} \\
			\hline
			std: partition 1 & 34.06 & 33.24 & 72.99 & 49.93          & 10.43          \\
			\hline
			dev: partition 0 & 47.16 & 37.32 & 81.45 & \textbf{60.16} & -              \\
			std: partition 0 & 47.57 & 36.75 & 81.56 & \textbf{60.45} & -              \\
			\hline
		\end{tabular}
		\caption{MUTAN without Attention model evaluation results.}

		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = 0.48  & Other & Num   & Y/N   & All            & Diff          \\ [0.5ex]
			\hline
			dev: partition 1 & 40.80 & 30.34 & 76.92 & \textbf{54.49} & \textbf{5.99} \\
			dev: partition 2 & 39.63 & 30.67 & 76.49 & \textbf{53.78} & \textbf{6.70} \\
			dev: partition 3 & 39.33 & 31.12 & 75.48 & \textbf{53.28} & \textbf{7.20} \\
			dev: partition 4 & 39.31 & 29.78 & 75.12 & \textbf{52.97} & \textbf{7.51} \\
			dev: partition 5 & 39.38 & 29.87 & 74.96 & \textbf{52.95} & \textbf{7.53} \\
			dev: partition 6 & 39.13 & 30.74 & 73.95 & \textbf{52.51} & \textbf{7.97} \\
			dev: partition 7 & 38.90 & 31.14 & 73.80 & \textbf{52.39} & \textbf{8.09} \\
			\hline
			std: partition 1 & 40.88 & 28.82 & 76.67 & 54.37          & 5.95          \\
			\hline
			dev: partition 0 & 49.14 & 38.35 & 79.63 & \textbf{60.48} & -             \\
			std: partition 0 & 49.15 & 36.52 & 79.45 & \textbf{60.32} & -             \\
			\hline
		\end{tabular}
		\caption{HieCoAtt (Alt, VGG19) model evaluation results.}

		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = 0.37  & Other & Num   & Y/N   & All            & Diff           \\ [0.5ex]
			\hline
			dev: partition 1 & 46.57 & 32.09 & 76.60 & \textbf{57.33} & \textbf{8.46}  \\
			dev: partition 2 & 45.83 & 32.43 & 75.29 & \textbf{56.47} & \textbf{9.32}  \\
			dev: partition 3 & 45.17 & 32.52 & 74.87 & \textbf{55.99} & \textbf{9.80}  \\
			dev: partition 4 & 45.11 & 32.31 & 73.73 & \textbf{55.47} & \textbf{10.32} \\
			dev: partition 5 & 44.35 & 31.95 & 72.93 & \textbf{54.74} & \textbf{11.05} \\
			dev: partition 6 & 43.75 & 31.21 & 72.03 & \textbf{54.00} & \textbf{11.79} \\
			dev: partition 7 & 43.88 & 32.59 & 71.99 & \textbf{54.19} & \textbf{11.60} \\
			\hline
			std: partition 1 & 46.11 & 31.46 & 76.84 & 57.25          & 8.43           \\
			\hline
			dev: partition 0 & 57.01 & 37.51 & 83.54 & \textbf{65.79} & -              \\
			std: partition 0 & 56.60 & 36.63 & 83.68 & \textbf{65.68} & -              \\
			\hline
		\end{tabular}
		\caption{MLB with Attention model evaluation results.}
	\end{subtable}
	\begin{subtable}{0.49\linewidth}
		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = 0.23  & Other & Num   & Y/N   & All            & Diff           \\ [0.5ex]
			\hline
			dev: partition 1 & 43.96 & 28.90 & 71.89 & \textbf{53.79} & \textbf{12.19} \\
			dev: partition 2 & 42.66 & 28.08 & 70.05 & \textbf{52.32} & \textbf{13.66} \\
			dev: partition 3 & 41.62 & 29.12 & 69.58 & \textbf{51.74} & \textbf{14.24} \\
			dev: partition 4 & 41.53 & 29.30 & 67.96 & \textbf{51.06} & \textbf{14.92} \\
			dev: partition 5 & 40.46 & 27.66 & 68.03 & \textbf{50.39} & \textbf{15.59} \\
			dev: partition 6 & 40.03 & 28.44 & 66.98 & \textbf{49.84} & \textbf{16.14} \\
			dev: partition 7 & 39.11 & 28.41 & 67.44 & \textbf{49.58} & \textbf{16.40} \\
			\hline
			std: partition 1 & 43.55 & 28.70 & 71.76 & 53.63          & 12.14          \\
			\hline
			dev: partition 0 & 56.73 & 38.35 & 84.11 & \textbf{65.98} & -              \\
			std: partition 0 & 56.29 & 37.47 & 84.04 & \textbf{65.77} & -              \\
			\hline
		\end{tabular}
		\caption{MUTAN with Attention model evaluation results.}

		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = \textbf{0.53} & Other & Num   & Y/N   & All            & Diff          \\ [0.5ex]
			\hline
			dev: partition 1         & 44.42 & 36.39 & 76.94 & \textbf{56.90} & \textbf{4.91} \\
			dev: partition 2         & 43.37 & 34.99 & 76.10 & \textbf{55.90} & \textbf{5.91} \\
			dev: partition 3         & 42.22 & 33.97 & 75.80 & \textbf{55.11} & \textbf{6.70} \\
			dev: partition 4         & 42.52 & 34.21 & 75.33 & \textbf{55.09} & \textbf{6.72} \\
			dev: partition 5         & 42.81 & 34.69 & 75.21 & \textbf{55.23} & \textbf{6.58} \\
			dev: partition 6         & 42.27 & 35.16 & 74.50 & \textbf{54.73} & \textbf{7.08} \\
			dev: partition 7         & 41.95 & 35.14 & 73.64 & \textbf{54.22} & \textbf{7.59} \\
			\hline
			std: partition 1         & 44.93 & 35.59 & 76.82 & 57.10          & 4.96          \\
			\hline
			dev: partition 0         & 51.77 & 38.65 & 79.70 & \textbf{61.81} & -             \\
			std: partition 0         & 51.95 & 38.22 & 79.95 & \textbf{62.06} & -             \\
			\hline
		\end{tabular}
		\caption{HieCoAtt (Alt, Resnet200) model evaluation results.}

		\begin{tabular}{c | c c c c | c}
			}\def\rscore/{$R_{\text{score}}$/ = 0.08  & Other & Num   & Y/N   & All            & Diff           \\ [0.5ex]
			\hline
			dev: partition 1 & 20.49 & 25.98 & 68.79 & \textbf{40.91} & \textbf{17.11} \\
			dev: partition 2 & 19.81 & 25.40 & 68.51 & \textbf{40.40} & \textbf{17.62} \\
			dev: partition 3 & 18.58 & 24.95 & 68.53 & \textbf{39.77} & \textbf{18.25} \\
			dev: partition 4 & 18.50 & 24.82 & 67.83 & \textbf{39.43} & \textbf{18.59} \\
			dev: partition 5 & 17.68 & 24.68 & 67.99 & \textbf{39.09} & \textbf{18.93} \\
			dev: partition 6 & 17.29 & 24.03 & 67.76 & \textbf{38.73} & \textbf{19.29} \\
			dev: partition 7 & 16.93 & 24.63 & 67.45 & \textbf{38.50} & \textbf{19.52} \\
			\hline
			std: partition 1 & 20.84 & 26.14 & 68.88 & 41.19          & 16.99          \\
			\hline
			dev: partition 0 & 43.40 & 36.46 & 80.87 & \textbf{58.02} & -              \\
			std: partition 0 & 43.90 & 36.67 & 80.38 & \textbf{58.18} & -              \\
			\hline
		\end{tabular}
		\caption{LSTM Q+I model evaluation results.}
	\end{subtable}
	\caption{Compares the accuracy and }\def\rscore/{$R_{\text{score}}$/ of six VQA models with increasing noise levels from YNBQD generated by }\def\lasso/{\textit{LASSO}/ evaluated on dev and std. The results are split by the question type; Numerical (Num), Yes/No (Y/N), or Other. We notice that the }\def\rscore/{$R_{\text{score}}$/ of some models under YNBQD is better than GBQD (}\def\eg/{{\em e.g.}/, HAR) and vice versa (}\def\eg/{{\em e.g.}/, MUA). So, the }\def\rscore/{$R_{\text{score}}$/ can only be compared in different models if they have the same BQD because it can show certain biases in the models under investigation that are related to the type of the BQD. For example, we can clearly see how HAR is more robust towards Yes/No questions than general questions. Whereas, HAV is apparently agnostic to this property (}\def\ie/{{\em i.e.}/, the answer type being Yes/No versus general).}
	\label{tbl:lasso_ynbqd}
\end{table*}

\section{Conclusion}
In this work, we propose a novel framework, a semantic similarity ranking method, two large-scale basic question datasets and robustness measure (}\def\rscore/{$R_{\text{score}}$/) as a benchmark to help the community build accurate \emph{and} robust VQA models.

\section{Acknowledgments}
This work was supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research and used the resources of the Supercomputing Laboratory at KAUST in Thuwal, Saudi Arabia.

{\small


\section{Introduction}
\label{sect:intro}

The chromosphere that lies between the photosphere and the hot corona is  often reveals jet like features, which are small in size and has very short lifetimes.  Though there have been a lot of
confusions with the relationship between the quiet sun limb
spicules, mottles observed on the quiet sundisk and the
dynamic fibrils (DFs)\citep{1992A&A...264..236G},  the three phenomena are related \citep{1994A&A...290..285T}.  The physical process of chromosphere is not well understood, especially the mass and energy transportation in dynamic fibrils. \citep[see][] {2012RSPTA.370.3129R}.The flows (up and down) from photosphere to chromosphere and the higher layer corona flaunt a variety of  oscillation modes like three- minute, five-minute and the seven-minute oscillations.  \citet{2014ApJ...786..137T}  presented strong evidence to the shock behavior of sunspot oscillations for the first time in sunspots in the transition region (TR) and in the chromosphere based on the IRIS observations.  The 5-minute oscillations \citep{1962ApJ...135..474L} happening over the sunspots have been studied extensively over a period of time and have been found that there is a decrease in amplitude with height and not easily detected in the higher layers above photosphere. The 3-minute oscillations could be due to the mode conversion of  magnetoacoustic waves \citep{2012A&A...539L...4S}.  The fibrilar structure of chromosphere reveals that magnetic fields play a significant role \citep{1908ApJ....28..315H}. The jet like features which are also known as spicules, mottles, and the DFs are observed at limb, quiet sun disk and in AR respectively. \\
Our main concentration is on the observations of DFs. The dynamics of DFs and its imprint on TR, and its oscillations are extensively studied \citep{1999SoPh..190..419D,2003ApJ...595L..63D,2003ApJ...590..502D,2004Natur.430..536D}.  A remarkable power of oscillations in the Dfs were determined by  \citet{2003ApJ...595L..63D}. This was discovered using high resolution data from SST. \citet{2004Natur.430..536D} and \citet{2006RSPTA.364..383D} build up a model which explains the reason for DFs to be the inclination of magnetic field lines. This inclination reduce acoustic cut off frequency permitting the longer period so called,  p-mode waves, to leak into the atmosphere where it forms shock and these shocks are the reason for driving the DFs. This was further investigated and proved by \citet{2006ApJ...647L..73H} using the similar data from SST that the DFs are caused by magnetoacoustic shock waves.
In this paper we performed an analysis on the motion of DFs using New Solar Telescope \citep[NST,][]{2010AN....331..636C} (now known as GST, Goode Solar Telecope, after July 2017 ) operating at Big Bear Solar Observatory (BBSO) by measuring its properties and described its temporal evolution. Our results supports the fact that the DFs are formed due to the shock waves in chromosphere that are driven by the flows that were convective and photospheric oscillations as reported by De Pontieu et al. (2004). Through analysis of dataset taken at NST, we find that the DFs and its properties could be a result of shock waves that are generated from the p-mode propagating upward.

   \begin{figure}
           \centering
           \includegraphics[width=15.0cm, angle=0]{ms0197fig1-eps-converted-to.pdf}
           \caption{\normalsize H$\alpha$ line core image showing the paths of DFs overplotted. The DF marked as '1' is used for further analysis. Tick marks are in Pixels. }
           \label{Fig1}
        \end{figure}

\section{Observations and data processing}
\label{sect:Obs}

The sunspot of NOAA AR 12132, located at S09E08 on 2014 August 5 from BBSO  \citep[NST,][]{2010AN....331..636C}, was chosen for studying the DFs (see Figure 1). The pointer  was (270$\arcsec$, −395$\arcsec$), targeting the sunspot of AR 12132. The observations were performed using NST during 18:20 UT$-$19:20 UT. The Chromospheric images were acquired every 23 s by scanning of the H$\alpha$ spectral line from the blue wing -1 \AA{} to red wing +1 \AA{} with a step of size 0.2 \AA{}. The FOV is about 70$\arcsec$, and has a pixel scale of 0.029$\arcsec$ pixel$^{-1}$. The data is used to investigate the umbral oscillations in chromospheric sunspot at  different solar altitudes observed on August 5. A  combination of a  5 \AA{} interference filter and a Fabry–Pérot etalon is used in the  VIS to get a  bandpass of  0.07 \AA{} at the H$\alpha$ line. We chose the first H$\alpha$ -1.0 \AA{} image as a  reference image to align all other images in this passband. The relative  shifts were recorded, and used  to register the images in the other passbands of H$\alpha$.

     \section{Data analysis and results}\label{sec:analysis}
            \begin{figure*}
            \centering
            \includegraphics[width=14.0cm, angle=0]{ms0197fig2-eps-converted-to.pdf}
            \caption{\normalsize Space time diagram of the H$\alpha$ line center image for DF marked as '1' in Figure 1. The three dotted lines indicates the umbral (a), penumbral (b) and superpenumbral (b) regions for which power spectrum is shown in Figure 7. }
            \label{Fig2}
            \end{figure*}
     \label{sect:des}

     A close look into the movie of H$\alpha$ line core images gives an insight into  fibrils located in the close proximity and in the vicinity of sunspot near the image center, a few of which were colligated with the superpenumbra. We investigate upon 40 such DFs for this study, wherein most AR DFs from the line core of H$\alpha$ images ensues almost a perfect parabolic paths. We opt for chosing the direction of a DF, or a bunch of DFs, manually using CRisp SPectral EXplorer, \citet[CRISPEX,][]{2012ApJ...750...22V} which helps in analyzing the multidimensional data cubes and an additional software called TANAT (Timeslice ANAlysis tools) for anlayzing and also for various measurements of data. A comprehensive use of this software not only made the seeking of sequence of images easier but also helped in tracking the events that were very clear while ascending and descending along their individual paths in 584 jet like features in H$\alpha$. The CRISPEX is a widget based versatile IDL tool for visual inspection and analysis of high resolution data. CRISPEX also provides  space-time diagram for all the linear as well as curved paths. This information could be stored and could be retrieved later for further analysis. We have produced the  space-time diagram for all the 40 DFs detected. The  space-time diagram for one of the DF is shown in Figure 2 and its evolution is shown in Figure 3. The figure shows the receding of the already risen DF along the same path. We have detected almost 584 such trajectories and did a parabolic fitting on it. Figure 4 gives an example of the fitting. From the fit, we calculated the parameters like deceleration, maximum velocity, maximum height, duration and kinetic energy as well.  We assume the electron density to be same throughout the chromosphere.
     This electron density  (${n_e}$) multiplied with the mass of the electron  (${m_e}$) gives
     density ($\rho$ = ${n_e}$ ${m_e}$ ) which means that the kinetic energy per unit volume would be
     equal to the square of the maximum velocity  (KE = ρ${v_{max}^2}$/2). The kinetic
     energy per unit volume of these jet like features may transform further into
     heat and might play a significant role in the heating of corona.

                       \begin{figure}
                       \centering
                       \includegraphics[width=13.0cm, angle=0]{ms0197fig3-eps-converted-to.pdf}
                       \caption{\normalsize The snapshots of H$\alpha$ line core showing the temporal evolution of DFs (marked as '1' in Figure 1) as it recedes. [This event is available as an animation in the electronic edition of the Journal]}
                       \label{Fig3}
                       \end{figure}

                                \begin{figure*}
                                  \centering
                                  	\hskip 0.2cm
                                  \includegraphics[width=5.35cm,height=5.3cm, angle=0]{ms0197fig4a-eps-converted-to.pdf} \includegraphics[width=5.35cm,height=5.3cm, angle=0]{ms0197fig4b-eps-converted-to.pdf} \\
                                 \includegraphics[width=5.2cm,height=4.8cm, angle=0]{ms0197fig4c-eps-converted-to.pdf}   \includegraphics[width=5.2cm,height=4.8cm, angle=0]{ms0197fig4d-eps-converted-to.pdf}

                                 \caption{\normalsize The extracted portion of x-t plot for sample DFs are shown in bottom panels. The DFs follows  almost perfect parabolic path. The yellow line indicates the bestfit used to derive the deceleration, maximum velocity, duration and the maximum height. The Parabolic fitting for corresponding trajectories are shown in top panels }
                                 \label{Fig4}
                                 \end{figure*}

                    \begin{figure*}
                    \centering
                   \includegraphics[width=9.0cm, angle=0]{ms0197fig5-eps-converted-to.pdf}
                   \caption{\normalsize Scatter plots of deceleration vs. maximum velocity (a),  maximum velocity vs. duration (b), and decceleration vs. duration (C). A clear linear correlation exists between the deceleration and the maximum velocity.}
                  \label{Fig5}
                   \end{figure*}
                   \clearpage

                 \begin{figure}
                 \centering
                 \includegraphics[width=15.0cm, angle=0]{ms0197fig6-eps-converted-to.pdf}
                 \caption{\normalsize Distributions of the maximum velocity (panel a), duration (panel b), decceleration (panel c), Kinetic energy (panel d) respectively. The histograms are constructed from the observations of DFs.}
                 \label{Fig6}
                 \end{figure}
                 \clearpage

        \begin{figure*}
        \centering
       \includegraphics[width=14.0cm, angle=0]{ms0197fig7-eps-converted-to.pdf}
       \caption{\normalsize The fourier power spectra at the positions marked in white dotted lines in Figure 2.}
       \label{Fig7}
       \end{figure*}

       We obtained the maximum velocity is between 10 km s$^{-1}$ and 30 km s$^{-1}$ and deceleration between 50 m s$^{-2}$ and 200 m s$^{-2}$. We noticed that  these correlations are quite similar to the measurements on dynamic fibrils done by \citet{2006ApJ...647L..73H}.
        But we found that the correlation between the duration and the deceleration/maximum velocity is comparatively weaker, showing a high scatter and nonlinearity in almost all the DFs.
                \begin{figure}
               \includegraphics[width=17.0cm, angle=0]{ms0197fig8-eps-converted-to.pdf}

              \caption{\normalsize Shock waves along the dynamic fibrils. The wavelength-time maps of the H$\alpha$ spectral line at different points along the DF as it propagates upward is shown in (a), (b) and (c)}
              \label{Fig8}
              \end{figure}
      It could be seen in \citep{2007ApJ...655..624D} that the correlation between the duration and  the deceleration/ maximum velocity is not so clear or in other words weak. Similar results have been obtained by \citet{2008ApJ...673.1194L} in DFs with the reason being the effect of projection. It could be probable  that there is no effect of projection seen in duration whereas the effect is more prominent in velocity and deceleration. It is very clear from scatter plot in Figure 5 that the correlation between duration and deceleration/ maximum velocity is ambiguous with huge scatter and difference in slopes.
       By tracing all the jet like features we obtain the histograms of max velocity, duration, deceleration, and the kinetic energy shown in Figure 6. From the histograms, it is detected that the maximum velocity, duration and decelerations follows a Gaussian distribution. The maximum velocities are within 10 km s$^{-1}$ to 30 km s$^{-1}$ and the average is about 15 km s$^{-1}$, the duration ranges from 2 min to 20 min and the average is about 11 min, the deceleration ranges from 10 m s$^{-2}$ to 200 m s$^{-2}$ and the average is around 100 m s$^{-2}$. We obtain  the average maximum height to be around 0.025 Mm. We find the kinetic energy of all the 584 jet like features and estimated its average to be about 220 gm cm $^{-1}$s$^{-2}$ and found it to follow a partial Gaussian distribution. \\
       We take the fourier power spectra along the DF (originates from the umbral boundary and propagates higher up into the superpenumbral region) at three different points along the DF as shown in Figure 7. The power spectra in Figure 7(a) is at the umbral boundary, (b) is at penumbral region and (c) is at superpenumbral region. The running waves travels through the boundary of umbra and penumbra and eventually disappears at the penumbral boundary. In the power spectra we see oscillations in DF varying with time. To measure the oscillation period, we use the FFT (Fast Fourier Transform) function using IDL. As shown in Figure 7, the frequency at the umbral boundary is $\approx$ 4.3 mHz and based on the frequency of the power spectra, the oscillation periods at the same point is 235 s. The results show that the oscillations within umbra are a typical 3-min oscillation with period T = 235 s, and the
      oscillations of the same within penumbra are a  typical 5-min oscillation with period T = 345 s and T=408 s. The period of umbra oscillations is almost one half of that of the penumbral oscillations. The corresponding periods of umbral and penumbral oscillations measured by using FFT function are T$\_$umbra = 235 s and T$\_$penumbra = 408 s respectively. These oscillations propagates higher up and further develops into shocks as shown in Figure 8 suggesting  oscillations are upward propagating shock waves. Each panel in Figure 8 corresponds to three different positions along the marked DF in Figure  It is quite possible that the dynamic fibrils are usually driven in union with the sunspot oscillations \citep{2014ApJ...789..108C}. The Oscillations propagates higher up and develops into shock which propels the fibrils quite a several times. Thus, we speculate that the shocks driven DFs are related to the sunspot oscillations. To be more precise, there exists a form of propagation of shock wave fronts that starts from the sunspot center and reaches the bottom of the fibrils and propagates higher up along the DF \citep{2014ApJ...789..108C}. Hence, it is clear that the DFs are physically related to the sunspot oscillations and the form of shock fronts propagates from the sunspot center. Figure 8 shows the temporal-spectral variations (t-$\lambda$ plots) indicating  the oscillations at three different positions along the DF (connecting the sunspot to the fibril ) to infact be the upward propagating shock waves. \citet{2006ApJ...647L..73H} and \citet{2007ApJ...655..624D} simulated upward propagating  magnetoacoustic shock waves and created the dynamic fibrils observational behaviour. They found that the highly  dynamic chromospheric shock waves cause significant upward  and downward motion of the upper chromosphere. The transition of an upward propagating shock through the chromosphere produces a form of sawtooth or N-shape (as seen in Figure 8 ), which marks an impression of blueshift followed by a gradual drift towards the redshift and then a sudden appearance of  a blueshift \citep{2006ApJ...647L..73H,2009A&A...494..269V}. The N-shape pattern preponderate in all the three t-$\lambda$ plots constructed at different locations along the DF and hence confirms the idea of shock driven fibrils. The correlation coefficient between the maximum velocity and deceleration of a sample of 40 sunspot dynamic fibrils is 0.757. This proves that the oscillations in sunspots are the magnetoacoustic shock waves that propagates upward     (e.g.,\citealt{1986ApJ...301..992L,2006ApJ...640.1153C,2014ApJ...789..108C,2014ApJ...786..137T,2014ApJ...787...58Y}).\\

   \section{Discussion and conclusion}\label{sec:discussion}

   The above found correlations could probably be a significant signature for the jets that are periodic, are perhaps driven by the waves that normally propagates from photosphere into the chromosphere and hence steepen into shocks \citet{2013ApJ...776...56R}. The maximum velocities of these jet like features that we have found in our data are usually between 10$-$30 km s$^{-1}$ and is found to be similar in comparison to \citet{2006ApJ...647L..73H,2007ApJ...655..624D}. The chromospheric waves are usually created by the usual convective flows and also due to the oscillations found in photosphere and also in the convection zone. These disturbances further propagates into the chromosphere upwards and thereby shock and becomes a driving force for the plasma in chromospheric region and hence results in the formation of the DFs.
    The periods that are longer than the period of local acoustic cut off is blocked by the chromosphere. This cut off period is completely dependent on the magnetic filed line inclination with respect to the vertical  \citep{1990LNP...367..211S,2004Natur.430..536D}. A single shock could possibly drive the DFs and is also the reason to understand the correlation between the deceleration and the maximum velocity. These DFs are the direct consequence of the upward propagating chromospheric oscillations/waves that are produced in the convection zone or photosphere as a result of global p-mode oscillations and also the convective flows. The waves propagates into the chromosphere and passes through upto the region where H$\alpha$ is formed, thus forming shocks thereby while its upward propagation. The fact on which we stress our interpretation is that the DFs follow a parabolic path with a good positive correlation between the properties like the velocity and the deceleration for all of the detected DFs. We found that velocity and deceleration are directly proportional with higher the velocity, the higher is the deceleration. The velocity that \citet{1995ApJ...450..411S} reports for their event of mottles is in the range of 10$-$30 km s$^{-1}$ and we found it to be rather similar to what we found in our BBSO/NST observations for our event of DFs. The same has been seen by \citet{2006ApJ...647L..73H} and is evident that the formation of DFs is due to the shocks in the chromosphere driven by the photospheric convective flows and oscillations as well. Hence, we conclude that the driving mechanism is same. \citet{2007ApJ...666.1277H}, \citet{2011ApJ...743..142H} simulation also explains that the usually long period waves propagates into the solar atmosphere along the magnetic field lines that are inclined. The most essential conclusion is with the explanation that the driving force for the DFs are the so called magnetoacoustic shocks that are induced by the p-mode oscillations and in addition leakage of convective flows into the chromosphere. Both the observations \citealt{1992ASIC..375..261L} and simulations \citealt{2010ApJ...722..888B} provides enough evidence for the waves that propagates upwards and also the shocks in the atmosphere of the sunspots. This phenomena appears as a dark feature in the H$\alpha$ images and are due to the density of magnetic flux.  \citet{2004Natur.430..536D} proposed that the similarity between the observations and some modeling done earlier on the DFs provides reason for the upward and downward movement of the upper chromosphere which is due to the shock waves in the chromosphere. \citet{2017ApJ...838....2Z} proves that the oscillations which appears to be like surge or light walls which lies above the light bridges are also caused by shocks . Our findings indicates that, in ARs,  most of these jet like features are caused by the shocks in chromosphere.

\begin{acknowledgements}

T.G.P would like to thank Patrick Antolin for his valuable help on CRISPEX.  This work is supported by the Grant 11427901, 11773038, 11373040, 11373044, 11273034, 11303048,11178005, AGS-0847126 and NSFC-1142830911427901 and also supported partly by  State Key Laboratory for Space Weather, Center for Space Science and Applied Research, Chinese Academy of Sciences. BBSO operation is supported by NJIT, US NSF AGS-1250818 and NASA NNX13AG14G. T.G.P thanks the financial support from CAS-TWAS Presidents PhD  fellowship - 2014. We thank referee for the review and for useful comments and suggestions.

\end{acknowledgements}


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

\section{Unsupervised CL-WSD Method}
\label{sec:methodology}
\input{methodology}

\section{Experiments Setup}
\label{sec:setup}
\input{setup}

\section{Results and Discussion}
\label{sec:experiments}
\input{results}

\section{Conclusion}
\label{sec:conclusion}
We study the application of word embedding-based methods in unsupervised Cross Language Word Sense Disambiguation (CL-WSD) when translating an English noun, appeared in a paragraph, to Persian. Our semantic approach uses embedding of the candidate translations as well as translated context terms to calculate the semantic similarity of each translation to its context. The proposed approach outperforms both the CO-Graph system---a state-of-the-art system in unsupervised CL-WSD---as well as the standard baseline.

We however observe fundamental limitations of the methods based exclusively on context as bag of words. Despite this fact, the current work offers a possible solution for all languages/domains with scarce knowledge-based or parallel corpora resources, by exploiting the use of a monolingual corpus together with a simple bilingual lexicon.

\bibliographystyle{abbrv}
\vspace{-0.4cm}


\section{Introduction}

\begin{table*}
\label{tab_elements}
\begin{minipage}{115mm}
\centering
\caption{Physical data and orbital elements of the  components of the triple system (45) Eugenia $^{(1)}$}
\end{minipage}
\begin{tabular}{ccccccc}
\hline
Body		&Orbits		&$a$		&$e$		&$i$ $^{(2)}$		&Radius (km)		&Mass (kg)  		       \\
\hline
Eugenia		&Sun		&$2.72$ au	&$0.083$	&$18.2^{\circ}$		&$108.5$$^{(3)}$  	&$5.63\times10^{18}$  	       \\\\
Princesse	&Eugenia	&$610.8$ km	&$0.069$	&$18^{\circ}$		&$2.5$  		&$2.5\times10^{14}$   	       \\\\
Petit-Prince	&Eugenia	&$1164.5$ km	&$0.006$	&$9^{\circ}$		&$3.5$  		&$2.5\times10^{14}$  	       \\\\
\hline
\multicolumn{7}{l}{(1) \cite{marchis} and JPL's Horizons system, for the epoch Julian date 2452980.0 JD.}\\
\multicolumn{7}{l}{(2) Inclinations of Petit-Prince and Princesse with respect to the equator of Eugenia and }\\
\multicolumn{7}{l}{inclination of Eugenia relative to the heliocentric ecliptic.}\\
\multicolumn{7}{l}{(3) Equivalent radius \citep{marchis}.}\\
\end{tabular}
\end{table*}

To know the composition and the dynamics of the asteroids may help us
to understand the process of formation of the Solar System and, so, the formation and composition of our own planet.
Besides, it is also known that asteroids can represent a threat to our planet.
These are some of the reasons why space missions aimed to visit asteroids in the Solar System is a very popular topic in current astrodynamics and astronomy studies.

Examples of some successful missions are the OSIRIS-REx (NASA), launched in September, 2016, and that must return to the Earth with a sample of the
NEA (101955) Bennu \citep{b15}, or the Hayabusa mission (JAXA), launched in May, 2003, which explored the NEA (25413) Itokawa \citep{b10} and returned to
Earth in June, 2010, with samples from the surface of the asteroid \citep{fujita}.

Multiple systems of asteroids are interesting targets for space missions, since they increase the range of possible scientific investigations.
The ASTER mission  - the First Brazilian Deep Space Mission \citep{sukhanov}, was concepted taking this advantage into account.
In that sense, searching for stable orbits around those bodies is very important, both in terms of science and engineering.
Regions of stable orbits may be places full of dust and/or small pieces of materials, as well as possible locations for other members of the system.
In a mission, they may also indicate good regions to place a spacecraft to observe the system.

The primary target of the ASTER mission is the NEA 2001 SN263, which is a triple system of asteroids.
The announcement of this mission has motivated studies aimed to characterize regions of stability of this system.
\cite{araujo2012} and \cite{araujo2015}, characterized stable and unstable regions around the components of this triple system,
for the prograde and retrograde cases, through numerical integrations of the gravitational N-body problem.
\cite{prado2014} mapped orbits for a small body orbiting the asteroid 2001 SN263. He considered the disturbance
received from all the perturbing forces individually. This study used a new concept for mapping orbits that
shows the relative importance of each force for a given orbit in the system. Such information helped to make a decision about which forces
need to be included in the model for a given accuracy and nominal orbit.

Thus, the present paper uses a combination of both methods presented in  \cite{araujo2012}, \cite{araujo2015} and \cite{prado2014},
in order to search for direct and retrograde stable orbits in the (45) Eugenia triple system of asteroids.
(45) Eugenia is a triple system of asteroids of the MBA (Main Belt Asteroids).
It is composed by the central body Eugenia and by the satellites Petit-Prince and Princesse. Physical characteristics and orbital elements of this system
are presented at Tab. $\ref{tab_elements}$. The physical and orbital data for Eugenia and its satellites were obtained from
\cite{marchis}. \cite{Beauvalet2014} updadet those data with minor corrections in the order of the error bar of the first work.
Although those corrections must be important from the observational point of view, they do not significantly affect our long-term analysis of the
stability of the particles within the system.

In particular, in this paper we searched for regions where direct orbits are unstable, but retrograde orbits are stable.
Those orbits are very good for a potential mission.
The probe can benefit of the stability of the orbit to minimize station-keeping maneuvers and, at the same time,
to travel in regions that are expected to be free of dust, so reducing the risk of collisions with natural debris.
This idea was proposed before in the literature for the system 2001 SN263 \cite{araujo2015}, and is now extended to the (45) Eugenia system.

The main reason for this extension is that those two systems are completely different from each other in terms of physical characteristics.
The 2001 SN263 is much smaller in sizes and distances, with the bodies very close to each other.
It means that the possible orbits for the spatial exploration of the system are strongly perturbed and stable orbits are very rare in many important regions.
The (45) Eugenia has components located at much larger distances,
so the general mappings of the stable orbits are very much different. Even the Kozai effects are reduced, allowing inclined orbits
in many situations where they did not appear in the system 2001 SN263. It is easy to quantity the differences, in terms
of the perturbation received by a small body in both systems, using the integral of the perturbing forces for one orbital period \citep{prado2013}.

The approach used here to find stable orbits is similar to the one used in the previous
studies related to the 2001 SN263 system, and it is described in Sec. \ref{sec_method}

\section{Method}
\label{sec_method}

The method adopted consisted in dividing the regions around the triple system
in four distinct internal regions: two regions around the two small satellites of the system and the regions between the orbits
of the two satellites (see Fig. \ref{fig_encounter}).

Regions $3$ and $4$ around the satellites Princesse and Petit-Prince were defined by the Hill's radius \citep{murray},
considering an approximation given by the two-body problem between Eugenia and each one of the satellites
separately. We found a region of about $29$ km where the gravitational perturbation from Petit-Prince is dominant over the
other bodies of the system. For Princesse this region has a radius of about $15$ km. These values have guided us to define how the particles
would be distributed within the system. This is better described in subsections \ref{sec_r1} - \ref{sec_r4}

The N-body gravitational problem was integrated in time using the Gauss-Radau numerical integrator \citep{everhart} for a time-span of 2 years.
Numerical integrations were performed considering a system composed by seven bodies: the Sun, the planets Earth, Mars and Jupiter,
the three components of the asteroid system and thousand of particles randomly distributed around these components,
including planar and inclined prograde and retrograde orbits.

Regarding the shape of the body, following \cite{araujo2012} and \cite{araujo2015}, we considered as a first approximation the oblateness of
Eugenia through the $J_2$ value.
According to \cite{marchis} the $J_2$ value for Eugenia is equal to $0.060 \pm 0.002$ . \cite{Beauvalet2014} fit a dynamical model to
simulate the observational data and to determine constraints in the dynamical parameter. From this fit they found a
$J_2$ value equal to $0.0589 \pm 0.0004$. Based on those values and on their error bar, we considered in the integrations $J_2=0.06$.
The difference is less than $2\%$, which do not significantly affect our long-term analysis of stability of particles within the system.

The results are expressed by plots showing the percentage of
particles that survived, for each set of initial conditions, as a function of the semi-major axis, eccentricity and inclination
of the initial orbit. The regions where all the particles survived along the time span of 2 years were called stable regions. On the other
hand, regions where no particles survived along this time span were called unstable regions.

In section \ref{sec_results} we present the results of the numerical integrations for each region and also their implication
in a space mission, comparing with the results previously found for the triple system 2001 SN263 in the context
of the ASTER mission.

\begin{figure}
\begin{center}
\includegraphics[scale=0.32]{sistema__figura.png}
\caption{Representation of the triple system 45 Eugenia and the regions of influence of each component.}
\label{fig_encounter}
\end{center}
\end{figure}

\section{Results}
\label{sec_results}

\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[scale=0.85]{fig_bert1.png}}
\subfigure[]{\includegraphics[scale=0.85]{fig_bert2.png}}
\subfigure[]{\includegraphics[scale=0.85]{fig_bert3.png}}
\caption{Perturbation Integral $(m/s)$ as a function of the semi-major axis of the orbit of a small body within the triple system,
considering perturbations coming from a) Princesse and Petit Prince for circular orbits, b) Eugenia and Petit-Prince and c)
Eugenia and Princesse.}
\label{fig_bert}
\end{center}
\end{figure}

A first look at both systems of asteroids indicates that they are quite different in physical terms.
The system (45) Eugenia has a central body that is much more massive with respect to its satellite bodies,
when compared with the 2001 SN263 system. In this system the three bodies have a more uniform distribution of masses among the bodies.
The distances are also different, with the bodies of the system 2001 SN263 being much closer to each other than the bodies that composes
the (45) Eugenia system. This combination of factors affect
very much the stability of the orbits, with important consequences in selecting possible orbits for a potential mission having this
triple system as a target.
Those situations will be explored in detail when showing the results obtained here.

A previous general study of those differences can be done using the integral of the perturbing forces acting in this system
(\cite{prado2013}, \cite{carvalho2014}, \cite{lara2016}, \cite{olivprado2014}, \cite{oliveiraetal2014}, \cite{sanchesetal2016}, \cite{sanchesetal2014},
\cite{santosetal2015}, \cite{shortetal2016}).

Basically, an integral of the perturbing forces acting in a small body within a given system is made over one orbital period.
The idea is to give a first approximation of the differences of the perturbation level in both systems.
\cite{prado2014} made this study for the 2001 SN263 system. This paper also showed that the perturbation coming from the two smaller members of
the system is several orders of magnitude larger than the other perturbations, so it is valid to make a first analysis of the perturbations
using only those forces. \cite{prado2014} describes well the level of perturbation, showing two peaks near the orbits of the smaller bodies.
This result is expected, since the perturbation increases when the small body gets closer to the perturbing bodies,
but this figure quantify the level of perturbation. It is noted that the magnitude of the perturbation is of the order of 0.025 m/s in the regions
between the orbits of the two smaller bodies.

The present paper makes a similar analysis for the (45) Eugenia, which results are  shown in Figs \ref{fig_bert}. Since the idea is to measure the effects
of the perturbations with respect to the gravity field of the main body, which is the most important comparison, a new type of integral is used.
This index is calculated by the integration of the perturbing forces acting in a specific orbit divided by the gravity field of the main body, over
one orbital period of the given orbit. The index used in \citep{prado2014} was calculated to compare different forces in the same system, and for orbits
around the main body, so the effects of dividing the perturbation by the gravity field of the main body are small. But, in the present case, where
the idea is to compare two different systems, as well as orbits around the main body and the two smaller bodies, it is very important to quantify the level
of perturbations compared to the gravity field of the main body.

Figure \ref{fig_bert}a considers orbits around the main body of the system, Eugenia, perturbed by the two natural satellites of the system, Princesse and
Petit Prince. It is noted the same two
peaks near the orbits of the perturbing bodies showed in 2001 SN263 system, as expected, but the magnitude of the perturbation levels in all the
regions are smaller. It is of the order of $10^{-8}~m/s$ in the regions between both smaller bodies, and even smaller after the orbit of the exterior natural satellite.
The peaks are not shown to keep a better scale, to see more detail, but they are of the order of $10^{-6}~m/s$. The same type of integral index was calculated for the 2001 SN263,
since \citep{prado2014} uses a slight different index, showing values of the order of $10^{-6}~m/s$ for
the minimum points and near $10^{-5}~m/s$ for the peaks. It clearly indicates the weaker effects of the perturbations, in the (45) Eugenia system,
from the two companion bodies, for orbits around the main body of the system. This is a mathematical indicative of much less perturbed orbits,
which will have several consequences on the choice of the locations to place a spacecraft to observe the system. It also shows the importance of using a
scalar index to quantify the level of perturbation received in a given orbit, so it is not necessary to limit the study to the expected general
behavior of the perturbation level, without a quantitative analysis of the forces.

The integral of the perturbing forces acting in the system give us an estimation on the stability of the system. It quantify the perturbations due to the
bodies of the system, indicating regions where particles or a small body  would be more or less perturbed. Nevertheless, it did not take into
account other dynamical effects such as mean motion resonances, the Kozai resonances for high-inclined orbits or the known increase of
stability of retrograde orbits. The results presented in Fig. $2$ show, as expected, that the perturbation from the satellites is small in the triple
system Eugenia. Thus, we expect the regions around the components of this system to be stable. In order to verify if this system is subject to the
perturbations cited above and how they differ from what is predicted by the integral method, we performed numerical simulations of particles distributed
within the system and verified their stability as a function of their initial orbit, as follows.

\subsection{Region 1}
\label{sec_r1}

\begin{figure*}
\mbox{%
\subfigure[]{\includegraphics[height=2.9cm]{I0R1.png}}\qquad
\subfigure[]{\includegraphics[height=2.9cm]{I30R1.png}}
\subfigure[]{\includegraphics[height=2.9cm]{I60R1.png}}}

\mbox{%
\subfigure[]{\includegraphics[height=2.9cm]{I90R1.png}}\qquad
\subfigure[]{\includegraphics[height=2.9cm]{I120R1.png}}
\subfigure[]{\includegraphics[height=2.9cm]{I150R1.png}}}

\mbox{%
\subfigure[]{\includegraphics[height=2.8cm]{I180R1.png}}}
\caption{Diagrams $(a \times e)$ of stability of region $1$ for the triple system (45) Eugenia, for a time span of $2$ years.
In this region the particles orbit Eugenia with orbital inclination: a) $I=0^{\circ}$, b) $I=30^{\circ}$, c) $I=60^{\circ}$,
d) $I=90^{\circ}$, e) $I=120^{\circ}$, f) $I=150^{\circ}$, g) $I=180^{\circ}$.
The color-coded scale indicate the percentage of particles that survives in this region, going from $0\%$ (instability) until $100\%$ (stability).
As explained in the text, for a particle with pericentric distance smaller than $243~km$ (on the left of the green lines) a model considering a more
realistic shape of Eugenia must be adopted.}
\label{fig_histogram1}
\end{figure*}

Region $1$ was defined as the region around Eugenia going from $110$ km ($1.01$ radius of Eugenia) to $595$ km
(limited by the orbit of Princesse, considering its region of influence).

The particles were distributed within this region orbiting Eugenia with the following orbital elements:
$110~km~\leq~a~\leq~595~km$, taken every $5~km$, $0\leq e \leq 0.5$, taken every $0.05$, $0^{\circ}~\leq~I~\leq~180^{\circ}$, taken every $30^{\circ}$.
For each combination $(a \times e)$, we considered 100 particles with random values of $f$, $\omega$ and $\Omega$, where $a$ is the semi-major axis,
$e$ is the eccentricity, $I$ is the inclination, $f$ is the true-anomaly, $\omega$ is the argument of the pericenter and $\Omega$ is the longitude of the ascending node.
This combination of initial conditions resulted in a total of approximately $108,000$ particles distributed within region $1$.

The stability and instability are defined by the number of particles that survived (no ejections or collisions) throughout the numerical integrations.
The collision was defined by the physical radius of the bodies. For region $1$, an ejection was considered every time that the relative distance
between a given particle and Eugenia was greater than $595$ km (the limit of the region $1$).

The results for region $1$ is presented in the diagrams of Fig. \ref{fig_histogram1}.
They show that approximately the whole region is stable. The instability appears only for the orbits really close
to Eugenia or Princesse, or for those orbits with high eccentricities such that the crossing of orbits become possible, or for those orbits such
that the pericenter of the orbit is inside Eugenia, leading to collisions.

They also show that there are no effects from the inclination of the orbit in the survival of the particles.
It is opposite to what was found in the system 2001 SN263, where the inclination played an important role in the evolution of
the orbits of the particles. \cite{araujo2012} showed, in Fig. $6$, that the Kozai effect reduces the regions of stable orbits
with the increase of the inclination, with a maximum effect at 60 degrees of inclination.
Studying retrograde orbits, \cite{araujo2015} showed, in Fig. $8$, that there are locations in the region $1$ with
stable retrograde and unstable direct orbits. Those regions are excellent locations to place a space probe willing
to observe the central body, as explained before.

Therefore, the (45) Eugenia system does not have this option for
selecting orbits to observe the main body of the system. It happens because the third-body perturbation
acting in a small body within this system is much smaller in the (45) Eugenia system, when compared to the 2001 SN263 system. This is due to two factors.
The third-body perturbation depends on the mass ratio between the perturbing and the central body. This ratio is
(mass of Princesse, the closest perturbing body)/(mass of Eugenia) = $(2.51 \times 10^{14})/(5.63 \times 10^{18}) = 4.46 \times 10^{-5}$ in the case of the (45)
Eugenia system and $1.06 \times 10^{-2}$ for the 2001 SN263 system, a very large difference. Another point is that both companion asteroids are located
very far from the main body, more than 600 km for the closest one, as shown in Tab $1$. This distance is about $3.8$ km in the 2001 SN263 system.
It means that a combination of those two facts make the Kozai effects to be much smaller and not able to
modify significantly the evolution of planar and inclined orbits.

The particles in region $1$ are those subjected to suffer the effects due to the irregular shape of Eugenia.
\cite{jiang} explored these effects on the two satellites Princesse and Petit-Prince.
Nevertheless, \cite{chanut} has shown how the relative error on the gravitational potential depends on the radial distance of a particle relative to
a central irregular shaped body. In their work the gravitational potential was determined for the asteroid 433 Eros modeled using the polyhedral model.
They compared the results obtained with a simpler model where Eros were considered a point of mass.
Their results show an error smaller than $10\%$ at a distance of about $1.6$ of the physical semimajor axis of the body.
Applying this result to Eugenia, considering the physical semimajor axis of Eugenia as $152~km$ \citep{jiang},  we estimate that
particles within the distance of about $243~km$ are those that may suffer the effects of the irregular shape of Eugenia.
Thus, for a particle with pericentric distance smaller than $243~km$ (on the left of the green lines shown in the plot available in Fig. \ref{fig_histogram1})
a model considering a more realistic shape of Eugenia must be adopted. For the other regions the differences are not significant.

\subsection{Region 2}
\label{sec_r2}

Region $2$ is the region between the orbits of the two satellites, and its was filled with particles orbiting Eugenia.
This region goes from $625.8$ km (limited by the orbit of Princesse, considering
its region of influence) to $1136.5$ km (limited by the orbit of Petit-Prince, considering
its region of influence).
The particles were distributed within this region orbiting Eugenia with
$630\leq a \leq 1135$ km, taken every $5$ km. The other elements were taken as before for region $1$.
This combination of initial conditions resulted in a total of approximately $112,000$ particles distributed within region $2$.

The results for region $2$ are presented in the diagrams of Fig. \ref{fig_histogram2}.
As found for region $1$, we see that approximately the whole region $2$ is stable. The instability appears only for the orbits that
cross the orbit of Petit-Prince.
The results for this region also show that, once again, there are no effects from the inclination in the survival of the particles,
as occurred in region 1. It is one more time different from what happened for the system 2001 SN263, as detailed in Fig. $7$
of \cite{araujo2012}, where the inclination has a strong effect in the evolution of the orbits of the particles,
with a maximum at 75 degrees. The study of retrograde orbits made in \cite{araujo2015}
showed a very large region of stable/retrograde and unstable/direct orbits (Fig. $8$).
It means that the (45) Eugenia system does not have this option for selecting orbits for a space mission also in this region.
The reason is the same explained before. The weak Kozai effect due to the smaller mass ratio and larger distances of the third perturbing body.

\begin{figure*}
\mbox{%
\subfigure[]{\includegraphics[height=2.9cm]{I0R2.png}}\qquad
\subfigure[]{\includegraphics[height=2.9cm]{I30R2.png}}
\subfigure[]{\includegraphics[height=2.9cm]{I60R2.png}}}
\mbox{%
\subfigure[]{\includegraphics[height=2.9cm]{I90R2.png}}\qquad
\subfigure[]{\includegraphics[height=2.9cm]{I120R2.png}}
\subfigure[]{\includegraphics[height=2.9cm]{I150R2.png}}}
\mbox{%
\subfigure[]{\includegraphics[height=2.9cm]{I180R2.png}}}
\caption{Diagrams $(a \times e)$ of stability of  region $2$ for the triple system (45) Eugenia, for a time span of $2$ years. In this region the particles orbit Eugenia with orbital inclinations: a) $I=0^{\circ}$, b) $I=30^{\circ}$, c) $I=60^{\circ}$, d) $I=90^{\circ}$, e) $I=120^{\circ}$, f) $I=150^{\circ}$, g) $I=180^{\circ}$.
The color-coded scale indicates the percentage of particles that survives in this region, going from $0\%$ (instability) until $100\%$ (stability).}
\label{fig_histogram2}
\end{figure*}

\subsection{Region 3}
\label{sec_r3}

Region $3$ was defined as the region around Princesse, going from $2.6$ km ($1.04$ radius of Princesse) until $15$ km (Hill's radius of Princesse with respect to Eugenia)

The particles were distributed within this region orbiting Princesse with $630\leq a \leq 1135$ km, taken every $5$ km.
The other elements were taken as before.
This combination of initial conditions resulted in a total of approximately $14,300$ particles distributed within region $3$.

The results for region $3$ is presented in the diagrams of Fig. \ref{fig_histogram3}.
The results for this region show that there are strong effects of the inclination in the survival of the particles.
The Kozai effect is now present because the third body perturbation is now very strong since the mass ratio between
the main body and the perturbing body is also $4.46 \times 10^{-5}$. Therefore, orbits with higher inclination are more stable than
the planar ones, when out of the region where there is the kozai resonance action.
Even though, the results are still different from the ones obtained for the system 2001 SN263, shown in \cite{araujo2012}.
For the system 2001 SN263 no stable orbits were found around the closer smaller body, due to the very strong perturbation coming
from the central body. Even retrograde stable orbits do not exist, as shown in \cite{araujo2015}.
It means that the (45) Eugenia system has excellent natural orbits for a space mission to observe the smaller body closer to the central body,
opposite to the 2001 SN263 system, which requires orbital control to observe this body.

It is clear, from Fig. \ref{fig_histogram3}, that retrograde circular orbits with semi-major radius from $7$ km to $10$ km,
depending on the initial eccentricity, are excellent
for an exploration mission to this system, because they are stable while the equivalent direct orbits (same size and shape) are unstable

\begin{figure*}
\mbox{%
\subfigure[]{\includegraphics[height=2.9cm]{I0R3.png}}\qquad
\subfigure[]{\includegraphics[height=2.9cm]{I30R3.png}}
\subfigure[]{\includegraphics[height=2.9cm]{I60R3.png}}}
\mbox{%
\subfigure[]{\includegraphics[height=2.9cm]{I120R3.png}}\qquad
\subfigure[]{\includegraphics[height=2.9cm]{I150R3.png}}
\subfigure[]{\includegraphics[height=2.9cm]{I180R3.png}}}
\caption{Diagrams $(a \times e)$ of stability of  region $3$ for the triple system (45) Eugenia, for a time span of $2$ years.
In this region the particles orbit Princesse with orbital inclinations: a) $I=0^{\circ}$, b) $I=30^{\circ}$, c) $I=60^{\circ}$,
d) $I=120^{\circ}$, e) $I=150^{\circ}$, f) $I=180^{\circ}$.
The color-coded scale indicates the percentage of particles that survives in this region, going from $0\%$ (instability) until $100\%$ (stability).
The diagram for $I=90^{\circ}$ is not presented, since all the particles in this region with this inclination were lost in $2$ years (collisions or ejections).}
\label{fig_histogram3}
\end{figure*}

\subsection{Region 4}
\label{sec_r4}

Region $4$ was defined as the region around Petit-Prince, going from $3.6$ km ($1.03$ radius of Petit-Prince) until
$28$ km (Hill's radius of Petit-Prince with respect to Eugenia)

The particles were distributed within this region orbiting Petit-Prince with $3.6 \leq a \leq 28$ km, taken every $5$ km.
The other elements were taken as before.
This combination of initial conditions resulted in a total of approximately $27,500$ particles distributed within region $4$.

The results for this region show smaller effects from the inclination in the survival of the particles, when compared to region $3$.
The Kozai effect is present, but it is very small. Although the third body still has a larger mass, the distance between the bodies are much larger,
compared to region $3$.
The results for the case of the system 2001 SN263, shown in Fig. $8$ of the
paper by \cite{araujo2012}, indicate much larger effects in the most distant body.
It happens due to the much smaller distances involved in the system 2001 SN263, which contributes to higher third-body perturbations.
The study of the retrograde orbits showed larger regions of stable orbits, as occurred for the system 2001 SN263, as shown in \cite{araujo2015}, Fig. $8$.
It means that the (45) Eugenia system has excellent locations to place a spacecraft willing to observe the smaller body which is far away
from to the central body.

\begin{figure*}
\mbox{%
\subfigure[]{\includegraphics[height=2.9cm]{I0R4.png}}\qquad
\subfigure[]{\includegraphics[height=2.9cm]{I30R4.png}}
\subfigure[]{\includegraphics[height=2.9cm]{I120R4.png}}}
\mbox{%
\subfigure[]{\includegraphics[height=2.9cm]{I150R4.png}}\qquad
\subfigure[]{\includegraphics[height=2.9cm]{I180R4.png}}}
\caption{Diagrams $(a \times e)$ of stability of  region $4$ for the triple system (45) Eugenia, for a time span of $2$ years.
In this region the particles orbit Petit-Prince with orbital inclinations: a) $I=0^{\circ}$, b) $I=30^{\circ}$, c) $I=120^{\circ}$,
d) $I=150^{\circ}$, e) $I=180^{\circ}$.
The color-coded scale indicates the percentage of particles that survives in this region, going from $0\%$ (instability) until $100\%$ (stability).
The diagrams for $I=60^{\circ}$ $I=90^{\circ}$ are not presented, since all the particles in this region with this inclination were lost in $2$ years (collisions or ejections).}
\label{fig_histogram}
\end{figure*}

The results for regions $1$, $2$, $3$ and $4$ showed the existence of stable orbits in the regions around the three bodies of the system. Those results are different from what
was obtained when studying the 2001 SN263 system, where no stable orbits were found near the smallest body of the system.
Several other important differences were found. In particular, there are no places around the main body where there are stable retrograde
orbits and that are unstable in the prograde orbits. It means that this strategy to choose orbits is
not available to observe the main body. On the other hand, situations like that are found around the smaller bodies of the system.
In general, the lower perturbation level that occur in all
the regions of the (45) Eugenia system makes this system much easier to place a spacecraft, from the point of view of longer
duration missions requiring less orbital control. This is a key point in deep space missions, due to the technical
difficulties of carrying more fuel and making frequent orbital maneuvers at this distance.

It is clear, from Fig. \ref{fig_histogram}, that retrograde circular orbits with semi-major radius from $14$ km to $18$ km,
depending on the initial eccentricity, are excellent
for an exploration mission to this system, because they are stable and the direct orbits of the same size are unstable.

\section{Conclusions}

A detailed study of the (45) Eugenia system was made, with the goal of finding stable orbits in the several regions of this triple asteroid system.
Those stable orbits are important to indicate regions of possible accumulation of dust, or even good locations to looking for a new member of
the family of the asteroids system. In that sense, stable regions were found, as a function of the semi-major axis, eccentricity and inclination of the orbits,
around all the bodies that composes the system.

In general, the results showed that the (45) Eugenia system is much less perturbed than the triple system 2001 SN263, which was studied before.
The integral of the perturbing forces quantified this fact in numbers. In particular, the Kozai effect of destroying inclined orbits was not
present in orbits around the central body, due to the small mass ratio and large distance of the third-body perturbation.
However, this effect is present in orbits around the two satellites.

Those stable orbits are also important in terms of astronautical applications, because regions of the space that has unstable/direct and
stable/retrograde orbits are good locations to place a spacecraft, since they give stability for the orbit and a low risk of collisions with particles.
A new type of integral index was defined, taking into account the perturbations coming from the other bodies of the systems compared to the body that the
spacecraft is orbiting. It showed a comparison of the perturbation levels, indicating that orbits around the main body of the (45) Eugenia system are
about $10^4$ times less perturbed than the equivalent orbits on the 2001 SN263 system. Regarding the companion bodies,
the internal body of the (45) Eugenia system is also about $10^4$ times less perturbed than the equivalent body of the 2001 SN263 system.
For the external bodies, the levels of perturbations are similar for both systems.

The results showed that orbits to locate a spacecraft using the difference in the stability
for direct and retrograde orbits do not exist in the (45) Eugenia system for the regions near the central body and between the orbits of
the two smaller bodies, as happened for the 2001 SN263 system. But, on the opposite side, such regions were found around the two smaller
bodies of the system, covering regions where no stable orbits were found for the 2001 SN263 system. It means that, in terms of observing the
smaller bodies of the system, the triple asteroid (45) Eugenia offers better locations to place a spacecraft, with cleaned stable retrograde orbits.
The numerical results showed that orbits with semi-major radius from $7$ km to $10$ km are very good choices for orbits around Princesse and
orbits with semi-major radius from $14$ km to $18$ km are very good choices for orbits around Petit-Prince.

\section{Acknowledgments}
This work was funded by CNPq Procs. 150378/2015-7, Proc. 301338/2016-7 , and by FAPESP Procs. 2016/14665-2 and 2011/08171-3.
This support is gratefully acknowledged.

\renewcommand{\refname}{REFERENCES}