\section{Introduction}
\seclabel{intro}

Despite not being exposed to explicit syntactic supervision, neural language models (LMs), such as recurrent neural networks, are able to generate fluent and natural sentences, suggesting that they induce syntactic knowledge about the language to some extent.
However, it is still under debate whether such induced knowledge about grammar is robust enough to deal with syntactically challenging constructions such as long-distance subject-verb agreement.
So far, the results for RNN language models (RNN-LMs) trained only with raw text are overall negative;
prior work has reported low performance on the challenging test cases \cite{marvin-linzen:2018:EMNLP} even with the massive size of the data and model \cite{van-schijndel-EtAl:2019:EMNLP1},
or argue the necessity of an architectural change to track the syntactic structure explicitly \cite{wilcox-etal-2019-structural,kuncoro-EtAl:2018:Long}.
Here the task is to evaluate whether a model assigns a higher likelihood on a grammatically correct sentence (\ref{ex:correct_orc}) over an incorrect sentence (\ref{ex:incorrect_orc}) that is minimally different from the original one \cite{Q16-1037}.

\vspace{2pt}\begin{exe}
 \ex\label{ex:orc}
 \begin{xlist}
  \ex[]{\label{ex:correct_orc}The author that the guards like \underline{laughs}.}
  \ex[*]{\label{ex:incorrect_orc}{The author that the guards like \underline{laugh}.}}
 \end{xlist}
\end{exe}
\vspace{2pt}

In this paper, to obtain a new insight into the syntactic abilities of neural LMs, in particular RNN-LMs, we perform a series of experiments under a different condition from the prior work.
Specifically, we extensively analyze the performance of the models that are exposed to explicit negative examples.
In this work, negative examples are the sentences or tokens that are grammatically incorrect, such as (\ref{ex:incorrect_orc}) above.

Since these negative examples provide a direct learning signal on the task at test time it may not be very surprising if the task performance goes up.
We acknowledge this, and argue that our motivation for this setup is to deepen understanding, in particular the limitation or the capacity of the current architectures, which we expect can be reached with such strong supervision.
Another motivation is engineering:
we could exploit negative examples in different ways, and establishing a better way will be of practical importance toward building an LM or generator that can be robust on particular linguistic constructions.

The first research question we pursue is about this latter point:
what is a better method to utilize negative examples that help LMs to acquire robustness on the target syntactic constructions?
Regarding this point, we find that adding additional token-level loss trying to guarantee a margin between log-probabilities for the correct and incorrect words (e.g., $\log p(\textrm{laughs} | h)$ and $\log p(\textrm{laugh} | h)$ for (\ref{ex:correct_orc})) is superior to the alternatives.
On the test set of \citet{marvin-linzen:2018:EMNLP}, we show that LSTM language models (LSTM-LMs) trained by this loss reach near perfect level on most syntactic constructions for which we create negative examples, with only a slight increase of perplexity about 1.0 point.

Past work conceptually similar to us is \citet{enguehard-etal-2017-exploring}, which, while not directly exploiting negative examples, trains an LM with additional explicit supervision signals to the evaluation task.
They hypothesize that LSTMs do have enough capacity to acquire robust syntactic abilities but the learning signals given by the raw text are weak, and show that multi-task learning with a binary classification task to predict the upcoming verb form (singular or plural) helps models aware of the target syntax (subject-verb agreement).
Our experiments basically confirm and strengthen this argument, with even stronger learning signals from negative examples, and we argue this allows us to evaluate the true capacity of the current architectures.
In our experiments (Section~\secref{exp}), we show that our margin loss achieves higher syntactic performance than their multi-task learning.

Another relevant work on the capacity of LSTM-LMs is \citet{kuncoro-etal-2019-scalable}, which shows that by distilling from syntactic LMs \cite{dyer-EtAl:2016:N16-1}, LSTM-LMs can improve their robustness on various agreement phenomena.
We show that our LMs with the margin loss outperform theirs in most of the aspects, further strengthening the argument about a stronger capacity of LSTM-LMs.

The latter part of this paper is a detailed analysis of the trained models and introduced losses.
Our second question is about the true {\it limitation} of LSTM-LMs:
are there still any syntactic constructions that the models cannot handle robustly even with our direct learning signals?
This question can be seen as a fine-grained one raised by \citet{enguehard-etal-2017-exploring} with a stronger tool and improved evaluation metric.
Among tested constructions, we find that syntactic agreement across an object relative clause (RC) is challenging.
To inspect whether this is due to the architectural limitation, we train another LM on a dataset, on which we unnaturally augment sentences involving object RCs.
Since it is known that object RCs are relatively rare compared to subject RCs \cite{hale-2001-probabilistic},
frequency may be the main reason for the lower performance.
Interestingly, even when increasing the number of sentences with an object RC by eight times (more than twice of sentences with a subject RC), the accuracy does not reach the same level as agreement across a subject RC.
This result suggests an inherent difficulty in tracking a syntactic state across an object RC for sequential neural architectures.

We finally provide an ablation study to understand the encoded linguistic knowledge in the models learned with the help of our method.
We experiment under reduced supervision at two different levels:
(1) at a lexical level, by not giving negative examples on verbs that appear in the test set;
(2) at a construction level, by not giving negative examples about a particular construction, e.g., verbs after a subject RC.
We observe no huge score drops by both.
This suggests that our learning signals at a lexical level (negative words) strengthen the abstract syntactic knowledge about the target constructions, and also that the models can generalize the knowledge acquired by negative examples to similar constructions for which negative examples are not explicitly given.
The result also implies that negative examples do not have to be complete and can be noisy, which will be appealing from an engineering perspective.

\section{Target Task and Setup}
The most common evaluation metric of an LM is perplexity.
Although neural LMs achieve impressive perplexity \cite{merity2018regularizing}, it is an average score across all tokens and does not inform the models' behaviors on linguistically challenging structures, which are rare in the corpus.
This is the primary motivation to separately evaluate the models' syntactic robustness by a different task.

\subsection{Syntactic evaluation task}
\seclabel{task}
As introduced in Section~\secref{intro}, the task for a model is to assign a higher probability to the grammatical sentence over the ungrammatical one, given a pair of minimally different sentences at a critical position affecting the grammaticality.
For example, (\ref{ex:correct_orc}) and (\ref{ex:incorrect_orc}) only differ at a final verb form, and to assign a higher probability to (\ref{ex:correct_orc}), models need to be aware of the agreement dependency between {\it author} and {\it laughs} over an RC.

\paragraph{\citet{marvin-linzen:2018:EMNLP} test set}
While initial work \cite{Q16-1037,N18-1108} has collected test examples from naturally occurring sentences, this approach suffers from the coverage issue, as syntactically challenging examples are relatively rare.
We use the test set compiled by \citet{marvin-linzen:2018:EMNLP}, which consists of synthetic examples (in English) created by a fixed vocabulary and a grammar.
This approach allows us to collect varieties of sentences with complex structures.

The test set is divided by the syntactic constructions appearing in each example.
Many constructions are different types of subject-verb agreement, including local agreement on different sentential positions (\ref{ex:simple}), and non-local agreement across different types of phrases.
Intervening phrases include prepositional phrases, subject RCs, object RCs, and coordinated verb phrases (\ref{ex:long}).
(\ref{ex:orc}) is an example of agreement across an object RC.
\vspace{2pt}\begin{exe}
 \ex\label{ex:simple}The senators \underline{smile/*smiles}.
\end{exe}
\begin{exe}
 \ex\label{ex:long}The senators like to watch television shows and \underline{are/*is} twenty three years old.
\end{exe}
\vspace{2pt}
Previous work has shown that non-local agreement is particularly challenging for sequential neural models \cite{marvin-linzen:2018:EMNLP}.

The other patterns are reflexive anaphora dependencies between a noun and a reflexive pronoun (\ref{ex:reflexive}), and
on negative polarity items (NPIs), such as {\it ever}, which requires a preceding negation word (e.g., {\it no} and {\it none}) at an appropriate scope (\ref{ex:npi}):
\vspace{2pt}\begin{exe}
 \ex\label{ex:reflexive}The authors hurt \underline{themselves/*himself}.
\end{exe}
\begin{exe}
 \ex\label{ex:npi}\underline{No/*Most} authors have {\it ever} been popular.
\end{exe}
\vspace{2pt}

Note that NPI examples differ from the others in that the context determining the grammaticality of the target word (No/*Most) does not precede it.
Rather, the grammaticality is determined by the following context.
As we discuss in Section~\secref{method}, this property makes it difficult to apply training with negative examples for NPIs for most of the methods studied in this work.

All examples above (\ref{ex:orc}--\ref{ex:npi}) are actual test sentences, and we can see that since they are synthetic some may sound somewhat unnatural.
The main argument behind using this dataset is that even not very natural, they are still strictly grammatical, and an LM equipped with robust syntactic abilities should be able to handle them as a human would do.

We use the original test set used in \citet{marvin-linzen:2018:EMNLP}.\footnote{We use the ``EMNLP2018'' templates in {\color{darkblue}https://github.com/BeckyMarvin/LM\_syneval}.}
See the supplementary materials of this for the lexical items and example sentences in each construction.

\subsection{Language models}
\seclabel{lm}
\paragraph{Training data}

Following the practice, we train LMs on the dataset not directly relevant to the test set.
Throughout the paper, we use an English Wikipedia corpus assembled by \citet{N18-1108}, which has been used as training data for the present task \cite{marvin-linzen:2018:EMNLP,kuncoro-etal-2019-scalable}, consisting of 80M/10M/10M tokens for training/dev/test sets.
It is tokenized and rare words are replaced by a single unknown token, amounting to the vocabulary size of 50,000.

\paragraph{Baseline LSTM-LM}
Since our focus in this paper is an additional loss exploiting negative examples (Section~\secref{method}), we fix the baseline LM throughout the experiments.
Our baseline is a three-layer LSTM-LM with 1,150 hidden units at internal layers trained with the standard cross-entropy loss.
Word embeddings are 400-dimensional, and input and output embeddings are tied \cite{DBLP:journals/corr/InanKS16}.
Deviating from some prior work \cite{marvin-linzen:2018:EMNLP,van-schijndel-EtAl:2019:EMNLP1},
we train LMs at sentence level as in sequence-to-sequence models \cite{sutskever2014sequence}.
This setting has been employed in some previous work \cite{kuncoro-EtAl:2018:Long,kuncoro-etal-2019-scalable}.\footnote{
On the other hand, the LSTM-LM of \citet{marvin-linzen:2018:EMNLP}, which is prepared by \citet{N18-1108}, is trained at document level through truncated backpropagation through time (BPTT) \cite{conf/icassp/MikolovKBCK11}.
Since our training regime is more akin to the task setting of syntactic evaluation, it may provide some advantage at test time.
}

Parameters are optimized by SGD.
For regularization, we apply dropout on word embeddings and outputs of every layer of LSTMs, with weight decay of 1.2e-6, and anneal the learning rate by 0.5 if the validation perplexity does not improve successively, checking every 5,000 mini-batches.
Mini-batch size, dropout weight, and initial learning rate are tuned by perplexity on the dev set of Wikipedia dataset.\footnote{Following values are found: mini-batch size: 128; initial learnin rate: 20.0; dropout weight on the word embedding layer and each output layer of LSTM: 0.1.}
Note that we tune these values for the baseline LSTM-LM and fix them across the experiments.

The size of our three-layer LM is the same as the state-of-the-art LSTM-LM at document-level \cite{merity2018regularizing}.
\citet{marvin-linzen:2018:EMNLP}'s LSTM-LM is two-layer with 650 hidden units and word embeddings.
Comparing two, since the word embeddings of our models are smaller (400 vs.~650) the total model sizes are comparable (40M for ours vs.~39M for theirs).
Nonetheless, we will see in the first experiment that our carefully tuned three-layer model achieves much higher syntactic performance than their model (Section~\secref{exp}), being a stronger baseline to our extensions, which we introduce next.

\section{Learning with Negative Examples}
\seclabel{method}

Now we describe four additional losses for exploiting negative examples.
The first two are existing ones, proposed for a similar purpose or under a different motivation.
As far as we know, the latter two have not appeared in past work.\footnote{
The loss for large-margin language models \cite{huang-EtAl:2018:EMNLP3} is similar to our sentence-level margin loss.
Whereas their formulation is more akin to the standard large-margin setting, aiming to learn a reranking model, our margin loss is simpler, just comparing two log-likelihoods of predefined positive and negative sentences.
}

We note that we create negative examples by modifying the original Wikipedia training sentences, not sentences in the test set.
As a running example, let us consider the case where sentence (\ref{ex:correct_batch}) exists in a mini-batch, from which we create a negative example (\ref{ex:incorrect_batch}).
\vspace{2pt}\begin{exe}
 \ex\label{ex:batch}
 \begin{xlist}
  \ex[]{\label{ex:correct_batch}An industrial park with several companies \underline{is} located in the close vicinity.}
  \ex[*]{\label{ex:incorrect_batch}An industrial park with several companies \underline{are} located in the close vicinity.}
 \end{xlist}
\end{exe}
\paragraph{Notations}
By a {\it target} word, we mean a word for which we create a negative example (e.g., {\it is}).
We distinguish two types of negative examples:
a {\it negative token} and a {\it negative sentence}; the former means a single incorrect word (e.g., {\it are}), while the latter means an entire ungrammatical sentence.

\subsection{Negative Example Losses}

\paragraph{Binary-classification loss}
This is proposed by \citet{enguehard-etal-2017-exploring} to complement a weak inductive bias in LSTM-LMs for learning syntax.
It is multi-task learning across the cross-entropy loss ($L_{lm}$) and an additional loss ($L_{add}$):
\begin{equation}\eqlabel{totalloss}
 L = L_{lm} + \beta L_{add},
\end{equation}
where $\beta$ is a relative weight for $L_{add}$.
Given outputs of LSTMs, a linear and binary softmax layers predict whether the next token is singular or plural.
$L_{add}$ is a loss for this classification, only defined for the contexts preceding a target token $x_{i}$:
\begin{equation*}

  L_{add} = \sum_{x_{1:i} \in \mathbf{h^*}} -\log p(\textrm{num}(x_{i}) | x_{1:i-1}),
\end{equation*}
where $x_{1:i} = x_1 \cdots x_{i}$ is a prefix sequence and $\mathbf{h^*}$ is a set of all prefixes ending with a target word (e.g., {\it An industrial park with several companies is}) in the training data.
$\textrm{num}(x) \in \{ \textrm{singular, plural} \}$ is a function returning the number of $x$.
In practice, for each mini-batch for $L_{lm}$, we calculate $L_{add}$ for the same set of sentences and add these two to obtain a total loss for updating parameters.

As we mentioned in Section~\secref{intro}, this loss does not exploit negative examples explicitly;
essentially a model is only informed of a key position (target word) that determines the grammaticality.
This is rather an indirect learning signal, and we expect that it does not outperform the other approaches.

\paragraph{Unlikelihood loss}
This is recently proposed \cite{welleck2019neural} for resolving the {\it repetition} issue, a known problem for neural text generators \cite{holtzman2019curious}.
Aiming at learning a model that can suppress repetition, they introduce an unlikelihood loss, which is an additional loss at a token level and explicitly penalizes choosing words previously appeared in the current context.

We customize their loss for negative tokens $x_i^*$ (e.g., {\it are} in (\ref{ex:incorrect_batch})).
Since this loss is added at token-level, instead of Eq.~\eqref{totalloss} the total loss is $L_{lm}$, which we modify as:
\begin{align}
 \sum_{\mathbf x \in D} \sum_{x_i \in \mathbf x} &-\log p(x_i|x_{1:i-1}) + \sum_{x_i^* \in \textrm{neg}_t(x_i)} g(x_i^*), \nonumber \\
 g(x_i^*) &= - \alpha \log (1 -p(x_i^* | x_{1:i-1})), \nonumber
\end{align}
where $\textrm{neg}_t(\cdot)$ returns negative tokens for a target $x_i$.\footnote{
Empty for non-target tokens.
It may return multiple tokens sometimes, e.g., themselves$\rightarrow$\{himself, herself\}.
}
$\alpha$ controls the weight.
$\mathbf x$ is a sentence in the training data $D$.
The unlikelihood loss strengthens the signal to penalize undesirable words in a context by explicitly reducing the likelihood of negative tokens $x_i^*$.
This is a more direct learning signal than the binary classification loss.

\paragraph{Sentence-level margin loss}
We propose a different loss, in which the likelihoods for correct and incorrect sentences are more tightly coupled.
As in the binary classification loss, the total loss is given by Eq.~\eqref{totalloss}.
We consider the following loss for $L_{add}$:
\begin{equation*}
 \sum_{\mathbf x \in D} \sum_{\mathbf x_j^* \in \textrm{neg}_s(\mathbf x)} \max(0, \delta - (\log p(\mathbf x) - \log p(\mathbf x_j^*))),
\end{equation*}
where $\delta$ is a margin value between the log-likelihood of original sentence $\mathbf x$ and negative sentences $\{ \mathbf x_j^* \}$.
$\textrm{neg}_s(\cdot)$ returns a set of negative sentences by modifying the original one.
Note that we change only one token for each $\mathbf x_j^*$, and thus may obtain multiple negative sentences from one $\mathbf x$ when it contains multiple target tokens (e.g., {\it she \underline{leaves} there but \underline{comes} back ...}).\footnote{
In principle, one can cumulate this loss within a single mini-batch for $L_{lm}$ as we do for the binary-classification loss.
However, obtaining $L_{add}$ needs to run an LM entirely on negative sentences as well, which demands a lot of GPU memories.
We avoid this by separating mini-batches for $L_{lm}$ and $L_{add}$.
We precompute all possible pairs of ($\mathbf x$, $\mathbf x_j^*$) and create a mini-batch by sampling from them.
We make the batch size for $L_{add}$ (the number of pairs) as the half of that for $L_{lm}$, to make the number of sentences contained in both kinds of batches equal.
Finally, in each epoch, we only sample at most the half mini-batches of those for $L_{lm}$ to reduce the total amount of training time.
}

Comparing to the unlikelihood loss, not only decreasing the likelihood of a negative example, this loss tries to guarantee a certain difference between the two likelihoods.
The learning signal of this loss seems stronger in this sense;
however, the token-level supervision is missing, which may provide a more direct signal to learn a clear contrast between correct and incorrect words.
This is an empirical problem we pursue in the experiments.

\paragraph{Token-level margin loss}
Our final loss is a combination of the previous two, by replacing $g(x_i)$ in the unlikelihood loss by a margin loss:
\begin{align}
 g(x_i^*) = \max(0, \delta-&(\log p(x_i|x_{1:i-1}) \nonumber \\
                 &- \log p(x_i^*|x_{1:i-1})). \nonumber
\end{align}
We will see that this loss is the most advantageous in the experiments (Section~\secref{exp}).

\subsection{Parameters}
Each method employs a few additional hyperparameters ($\beta$ for the binary classification loss, $\alpha$ for the unlikelihood loss, and $\delta$ for the margin losses).
We preliminary select $\beta$ and $\alpha$ from $\{1,10,100,1000\}$ that achieve the best average syntactic performance and find $\beta=1$ and $\alpha=1000$.
For the two margin losses, we fix $\beta=1.0$ and $\alpha=1.0$ and only see the effects of margin value $\delta$.

\begin{table*}
\centering
\scalebox{0.74}{
\begin{tabular}{lrrrrrrr}
\toprule
 & \multicolumn{2}{c}{LSTM-LM} & \multicolumn{2}{c}{Additional margin loss ($\delta=10$)} & \multicolumn{2}{c}{Additional loss {\small ($\alpha=1000, \beta=1$)}} & \multicolumn{1}{c}{Distilled} \\
 \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} \cmidrule{8-8}
 & \multicolumn{1}{c}{M\&L18}&\multicolumn{1}{c}{Ours}&\multicolumn{1}{c}{Sentence-level}&\multicolumn{1}{c}{Token-level}&\multicolumn{1}{c}{Binary-pred.}&\multicolumn{1}{c}{Unlike.} & \multicolumn{1}{c}{K19} \\
\midrule
\textsc{Agreement}: &&&&&&\\
Simple & 94.0&98.1 ($\pm$1.3)&\textbf{100.0} ($\pm$0.0)&\textbf{100.0} ($\pm$0.0)&99.1 ($\pm$1.2)&99.7 ($\pm$0.6)&\textbf{100.0} ($\pm$0.0) \\
In a sent. complement & 99.0&96.1 ($\pm$2.0)&95.8 ($\pm$0.7)&\textbf{99.3} ($\pm$0.4)&96.9 ($\pm$2.4)&92.7 ($\pm$3.1)&98.0 ($\pm$2.0) \\
Short VP coordination & 90.0&93.6 ($\pm$3.0)&\textbf{100.0} ($\pm$0.0)&99.4 ($\pm$1.1)&93.8 ($\pm$3.3)&95.6 ($\pm$3.0)&99.0 ($\pm$2.0) \\
Long VP coordination & 61.0&82.2 ($\pm$3.4)&94.5 ($\pm$1.0)&\textbf{99.0} ($\pm$0.8)&83.9 ($\pm$3.2)&90.0 ($\pm$2.4)&80.0 ($\pm$2.0) \\
Across a PP & 57.0&92.6 ($\pm$1.4)&\textbf{98.8} ($\pm$0.4)&98.6 ($\pm$0.3)&92.7 ($\pm$1.3)&95.2 ($\pm$1.2)&91.0 ($\pm$3.0) \\
Across a SRC & 56.0&91.5 ($\pm$3.4)&99.6 ($\pm$0.4)&\textbf{99.8} ($\pm$0.2)&91.9 ($\pm$2.5)&97.1 ($\pm$0.7)&90.0 ($\pm$2.0) \\
Across an ORC & 50.0&84.5 ($\pm$3.1)&93.5 ($\pm$4.0)&\textbf{93.7} ($\pm$2.0)&86.3 ($\pm$3.2)&88.7 ($\pm$4.1)&84.0 ($\pm$3.0) \\
Across an ORC (no that) & 52.0&75.7 ($\pm$3.3)&86.7 ($\pm$4.2)&\textbf{89.4} ($\pm$2.7)&78.6 ($\pm$4.0)&86.4 ($\pm$3.5)&77.0 ($\pm$2.0) \\
In an ORC & 84.0&84.3 ($\pm$5.5)&99.8 ($\pm$0.2)&\textbf{99.9} ($\pm$0.1)&89.3 ($\pm$6.2)&92.4 ($\pm$3.5)&92.0 ($\pm$4.0) \\
In an ORC (no that) & 71.0&81.8 ($\pm$2.3)&97.0 ($\pm$1.0)&\textbf{98.6} ($\pm$0.9)&83.0 ($\pm$5.1)&88.9 ($\pm$2.4)&92.0 ($\pm$2.0) \\
\midrule
\textsc{Reflexive}: &&&&&\\
Simple & 83.0&94.1 ($\pm$1.9)&99.4 ($\pm$1.1)&\textbf{99.9} ($\pm$0.2)&91.8 ($\pm$2.9)&98.0 ($\pm$1.1)&91.0 ($\pm$4.0) \\
In a sent. complement & 86.0&80.8 ($\pm$1.7)&\textbf{99.2} ($\pm$0.6)&97.9 ($\pm$0.8)&79.0 ($\pm$3.1)&92.6 ($\pm$2.9)&82.0 ($\pm$3.0) \\
Across an ORC & 55.0&74.9 ($\pm$5.0)&72.8 ($\pm$2.4)&73.9 ($\pm$1.3)&72.3 ($\pm$3.0)&\textbf{78.9} ($\pm$8.6)&67.0 ($\pm$3.0) \\
\midrule
\textsc{NPI}: &&&&&\\
Simple & 40.0&\textbf{99.2} ($\pm$0.7)&98.7 ($\pm$1.6)&97.7 ($\pm$2.0)&98.0 ($\pm$3.1)&98.2 ($\pm$1.2)&94.0 ($\pm$4.0) \\
Across an ORC & 41.0&63.5 ($\pm$15.0)&56.8 ($\pm$6.0)&64.1 ($\pm$13.8)&64.5 ($\pm$14.0)&48.5 ($\pm$6.4)&\textbf{91.0} ($\pm$7.0) \\
\midrule
Perplexity & 78.6&\textbf{49.5} ($\pm$0.2)&56.4 ($\pm$0.5)&50.4 ($\pm$0.6)&49.6 ($\pm$0.3)&50.3 ($\pm$0.2)&56.7 ($\pm$0.2) \\
\bottomrule
\end{tabular}
}
 \caption{Comparison of syntactic dependency evaluation accuracies across different types of dependencies and perplexities.
 Numbers in parentheses are standard deviations.
 M\&L18 is the result of two-layer LSTM-LM in \citet{marvin-linzen:2018:EMNLP}.
 K19 is the result of distilled two-layer LSTM-LM from RNNGs \cite{kuncoro-etal-2019-scalable}.
 VP: verb phrase; PP: prepositional phrase; SRC: subject relative clause; and ORC: object-relative clause.
 Margin values are set to 10, which works better according to Figure~\figref{margin}.
 Perplexity values are calculated on the test set of the Wikipedia dataset.
 The values of M\&L18 and K19 are copied from \citet{kuncoro-etal-2019-scalable}.
 }\tablabel{main}
\end{table*}

\begin{figure*}[t]
 \centering
 \scalebox{0.97}{
 \includegraphics[width=\linewidth]{figures/margin_vs_avg.pdf}}
 \caption{Margin value vs.~macro average accuracy over the same type of constructions,
 or perplexity, with standard deviation for the sentence and token-level margin losses.
 $\delta=0$ is the baseline LSTM-LM without additional loss.
 }\figlabel{margin}
\end{figure*}

\subsection{Scope of Negative Examples}
\seclabel{scope}
Since our goal is to understand to what extent LMs can be sensitive to the target syntactic constructions by giving explicit supervision via negative examples, we only prepare negative examples on the constructions that are directly tested at evaluation.
Specifically, we mark the following words in the training data, and create negative examples:
\begin{description}
 \item[Present verb] To create negative examples on subject-verb agreement, we mark all present verbs and change their numbers.\footnote{
            We use Stanford tagger \cite{toutanova2003feature} to find the present verbs.
            We change the number of verbs tagged by VBZ or VBP using \texttt{inflect.py} ({\color{darkblue}https://pypi.org/project/inflect/}).}
 \item[Reflexive pronoun] We also create negative examples on reflexive anaphora, by flipping between \{\textit{themselves}\}$\leftrightarrow$\{\textit{himself, herself}\}.
\end{description}
These two are both related to the syntactic number of a target word.
For binary classification we regard both as a target word, apart from the original work that only deals with subject-verb agreement \cite{enguehard-etal-2017-exploring}.
We use a single common linear layer for both constructions.

In this work, we do not create negative examples for NPIs.
This is mainly for technical reasons.
Among four losses, only the sentence-level margin loss can correctly handle negative examples for NPIs, essentially because other losses are token-level.
For NPIs, left contexts do not have information to decide the grammaticality of the target token (a quantifier; no, most, etc.) (Section~\secref{task}).
Instead, in this work, we use NPI test cases as a proxy to see possible negative (or positive) impacts as compensation for specially targeting some constructions.
We will see that in particular for our margin losses, such negative effects are very small.

\section{Experiments on Additional Losses}
\seclabel{exp}

We first see the overall performance of baseline LSTM-LMs as well as the effects of additional losses.
Throughout the experiments, for each setting, we train five models from different random seeds and report the average score and standard deviation.
The code is available at {\color{darkblue}https://github.com/aistairc/lm\_syntax\_negative}.

\paragraph{Naive LSTM-LM performs well}
The main accuracy comparison across target constructions for different settings is presented in Table~\tabref{main}.
We first notice that our baseline LSTM-LM (Section~\secref{lm}) performs much better than \citet{marvin-linzen:2018:EMNLP}'s LM.
A similar observation is recently made by \citet{kuncoro-etal-2019-scalable}.\footnote{We omit the comparison but the scores are overall similar.}
This suggests that the original work underestimates the true syntactic ability induced by LSTM-LMs.
The table also shows the results by their distilled LSTM-LM from RNNGs (Section~\secref{intro}).

\begin{figure*}[t]
 \centering
 \scalebox{0.97}{
 \includegraphics[width=\linewidth]{figures/boost_orc.pdf}}
 \caption{Accuracies on ``Across an ORC'' (with and without complementizer ``that'') by models trained on augmented data with additional sentences containing an object RC.
 Margin is set to 10.
 X-axis denotes the total number of object RCs in the training data.
 0.37M roughly equals the number of subject RCs in the original data.
 ``animate only'' is a subset of examples (see body).
 Error bars are standard deviations across 5 different runs.
 }\figlabel{orc}
\end{figure*}

\paragraph{Higher margin value is effective}
For the two types of margin loss, which margin value should we use?
Figure~\figref{margin} reports average accuracies within the same types of constructions.
For both token and sentence-levels, the task performance increases along $\delta$, but a too large value (15) causes a negative effect, in particular on reflexive anaphora.
Increases (degradations) of perplexity are observed in both methods but this effect is much smaller for the token-level loss.
In the following experiments, we fix the margin value to 10 for both, which achieves the best syntactic performance.

\paragraph{Which additional loss works better?}
We see a clear tendency that our token-level margin loss achieves overall better performance.
Unlikelihood loss does not work unless we choose a huge weight parameter ($\alpha=1000$), but it does not outperform ours, with a similar value of perplexity.
The improvements by binary-classification loss are smaller, indicating that the signals are weaker than other methods with explicit negative examples.
Sentence-level margin loss is conceptually advantageous in that it can deal with any type of sentence-level grammaticality including NPIs.
We see that it is overall competitive with token-level margin loss but suffers from a larger increase of perplexity (4.9 points), which is observed even with smaller margin values (Figure~\figref{margin}).
Understanding the cause of this degradation as well as alleviating it is an important future direction.

\section{Limitations of LSTM-LMs}
\seclabel{orc}
In Table~\tabref{main}, the accuracies on dependencies across an object RC are relatively low.
The central question in this experiment is whether this low performance is due to the limitation of current architectures, or other factors such as frequency.
We base our discussion on the contrast between object (\ref{ex:author_orc}) and subject (\ref{ex:author_src}) RCs:
\vspace{2pt}\begin{exe}
 \ex\label{ex:author_orc}The authors (that) the chef likes \underline{laugh}.
\end{exe}
\begin{exe}
 \ex\label{ex:author_src}The authors that like the chef \underline{laugh}.
\end{exe}
\vspace{2pt}
Importantly, the accuracies for a subject RC are more stable, reaching 99.8\% with the token-level margin loss, although the content words used in the examples are common.\footnote{\label{foot:animate}
Precisely, they are not the same.
Examples of object RCs are divided into two categories by the animacy of the main subject ({\it animate} or not), while subject RCs only contain animate cases.
If we select only animate examples from object RCs the vocabularies for both RCs are the same,
remaining only differences in word order and inflection, as in (\ref{ex:author_orc}, \ref{ex:author_src}).
}

It is known that object RCs are less frequent than subject RCs \cite{hale-2001-probabilistic,Levy2008-LEVESC}, and it could be the case that the use of negative examples still does not fully alleviate this factor.
Here, to understand the true limitation of the current LSTM architecture, we try to eliminate such other factors as much as possible under a controlled experiment.

\paragraph{Setup}
We first inspect the frequencies of object and subject RCs in the training data, by parsing them with the state-of-the-art Berkeley neural parser \cite{kitaev-klein:2018:Long}.
In total, while subject RCs occur 373,186 times, object RCs only occur 106,558 times.
We create three additional training datasets by adding sentences involving object RCs to the original Wikipedia corpus (Section~\secref{lm}).
To this end, we randomly pick up 30 million sentences from Wikipedia (not overlapped to any sentences in the original corpus), parse by the same parser, and filter sentences containing an object RC, amounting to 680,000 sentences.
We create augmented training sets by adding a subset, or all of these sentences to the original training sentences.
Among the test cases about object RCs we only report accuracies on subject-verb agreement, on which the portion for subject RCs also exists.
This allows us to compare the difficulties of two types of RCs for the present models.
We also evaluate on ``animate only'' subset, which has a correspondence to the test cases for subject RCs with only differences in word order and inflection (like (\ref{ex:author_orc}) and (\ref{ex:author_src}); see footnote \ref{foot:animate}).
Of particular interest to us is the accuracy on these animate cases.
We expect that the main reason for lower performance for object RCs is due to frequency, and with our augmentation the accuracy will reach the same level as that for subject RCs.

\begin{figure}[t]
 \centering
 \scalebox{0.9}{
 \includegraphics[width=\linewidth]{figures/exclude_token_pattern.pdf}}
 \caption{An ablation study to see the performance of models trained with reduced explicit negative examples (token-level and construction-level).
 One color represents the same models across plots, except the last bar (construction-level), which is different for each plot.}\figlabel{ablation}
\end{figure}

\paragraph{Results}
However, for both all and animate cases, accuracies are below those for subject RCs (Figure~\figref{orc}).
Although we see improvements from the original score (93.7), the highest average accuracy by the token-level margin loss on the ``animate'' subset is 97.1 (``with that''), not beyond 99\%.
This result indicates some architectural limitations of LSTM-LMs in handling object RCs robustly at a near perfect level.
Answering why the accuracy does not reach (almost) 100\%, perhaps with other empirical properties or inductive biases \cite{khandelwal-etal-2018-sharp,ravfogel-etal-2019-studying} is future work.

\section{Do models generalize explicit supervision, or just memorize it?}

One distinguishing property of our margin loss, in particular token-level loss, is that it is highly lexical, making a contrast explicitly between correct and incorrect words.
This direct signal may make models acquire very specialized knowledge about each target word, not very generalizable one across similar words and occurring contexts.
In this section, to get insights into the transferability of syntactic knowledge induced by our margin losses, we provide an ablation study by removing certain negative examples during training.

\paragraph{Setup}
We perform two kinds of ablation.
For token-level ablation (\textsc{-Token}), we avoid creating negative examples for all verbs that appear as a target verb\footnote{{\it swim, smile, laugh, enjoy, hate, bring, interest, like, write, admire, love, know}, and {\it is}.} in the test set.
Another is construction-level (\textsc{-Pattern}), by removing all negative examples occurring in a particular syntactic pattern.
We ablate a single construction at a time for \textsc{-Pattern}, from four non-local subject-verb dependencies (across a prepositional phrase (PP), subject RC, object RC, and long verb phrase (VP)).\footnote{We identify all these cases from the parsed training data, which we prepared for the analysis in Section~\secref{orc}.}
We hypothesize that models are less affected by token-level ablation, as knowledge transfer across words appearing in similar contexts is promoted by language modeling objective.
We expect that construction-level supervision would be necessary to induce robust syntactic knowledge, as perhaps different phrases, e.g., a PP and a VP, are processed differently.

\begin{table}
\centering
\scalebox{0.78}{
\begin{tabular}{lccc}
\toprule
& \multicolumn{3}{c}{Second verb (V1 and \textbf{V2})} \\
Models & All verbs & like & other verbs \\
\midrule
LSTM-LM & 82.2 ($\pm$3.4)&13.0 ($\pm$12.2)&89.9 ($\pm$3.6) \\
Margin (token) & 99.0 ($\pm$0.8)&94.0 ($\pm$6.5)&99.6 ($\pm$0.5) \\
~~ \textsc{-Token} & 90.8 ($\pm$3.3)&51.0 ($\pm$29.9)&95.2 ($\pm$2.6) \\
~~ \textsc{-Pattern} & 90.1 ($\pm$4.6)&50.0 ($\pm$30.6)&94.6 ($\pm$2.2) \\
\bottomrule
\end{tabular}
}
 \caption{Accuracies on long VP coordinations by the models with/without ablations.
 ``All verbs'' scores are overall accuracies.
 ``like'' scores are accuracies on examples on which the second verb (target verb) is {\it like}.
 }\tablabel{vp_coord_first}
\end{table}
\begin{table}[t]
\centering
\scalebox{0.78}{
\begin{tabular}{lcc}
\toprule
& \multicolumn{2}{c}{First verb (\textbf{V1} and V2)} \\
Models & likes & other verbs \\
\midrule
LSTM-LM & 61.5 ($\pm$20.0)& 93.5 ($\pm$3.4) \\
Margin (token) & 97.0 ($\pm$4.5)&99.9 ($\pm$0.1) \\
~~ \textsc{-Token} &63.5 ($\pm$18.5)& 99.2 ($\pm$1.1) \\
~~ \textsc{-Pattern} &67.0 ($\pm$21.2)& 98.0 ($\pm$1.4) \\
\bottomrule
\end{tabular}
}
 \caption{Further analysis of accuracies on the ``other verbs'' cases of Table~\tabref{vp_coord_first}.
 Among these cases, the second column (``likes'') shows accuracies on examples where the first verb (not target) is {\it likes}.
 }\tablabel{second_vp}
\end{table}

\paragraph{Results}
Figure~\figref{ablation} is the main results.
Across models, we restrict the evaluation on four non-local dependency constructions, which we select as ablation candidates as well.
For a model with \textsc{-Pattern}, we evaluate only on examples of construction ablated in training (see caption).
To our surprise, both \textsc{-Token} and \textsc{-Pattern} have similar effects, except ``Across an ORC'', on which the degradation by \textsc{-Pattern} is larger.
This may be related to the inherent difficulty of object RCs for LSTM-LMs that we verified in Section~\secref{orc}.
For such particularly challenging constructions, models may need explicit supervision signals.
We observe lesser score degradation by ablating prepositional phrases and subject RCs.
This suggests that, for example, the syntactic knowledge strengthened for prepositional phrases with negative examples could be exploited to learn the syntactic patterns about subject RCs, even when direct learning signals on subject RCs are missing.

We see approximately 10.0 points score degradation on long VP coordination by both ablations.
Does this mean that long VPs are particularly hard in terms of transferability?
We find that the main reasons for this drop, relative to other cases, are rather technical, essentially due to the target verbs used in the test cases.
See Table~\tabref{vp_coord_first,second_vp}, which show that failed cases for the ablated models are often characterized by the existence of either {\it like} or {\it likes}.
Excluding these cases (``other verbs'' in Table~\tabref{second_vp}), the accuracies reach 99.2 and 98.0 by \textsc{-Token} and \textsc{-Pattern}, respectively.
These verbs do not appear as a target verb in the test cases of other tested constructions.
This result suggests that the transferability of syntactic knowledge to a particular word may depend on some characteristics of that word.
We conjecture that the reason for weak transferability to {\it likes} and {\it like} is that they are polysemous;
e.g., in the corpus, {\it like} is much more often used as a preposition and being used as a present tense verb is rare.
This type of issue due to frequency may be one reason for lessening the transferability.
In other words, {\it like} can be seen as a challenging verb to learn its usage only from the corpus, and our margin loss helps for such cases.

\section{Discussion and Conclusion}
Our results with explicit negative examples are overall positive.
We have demonstrated that models exposed to these examples at training time in an appropriate way will be capable of handling the targeted constructions at near perfect level except a few cases.
We found that our new token-level margin loss is superior to the other approaches and the remaining challenging cases are dependencies across an object relative clause.

Object relative clauses are known to be harder for a human as well, and our results may indicate some similarities in the sentence processing behaviors by a human and RNN, though other studies also find some dissimilarities between them \cite{linzen2018distinct,wilcox-etal:2019-what-syntactic-structures}.
The difficulty of object relative clauses for RNN-LMs has also been observed in the prior work \cite{marvin-linzen:2018:EMNLP,van-schijndel-EtAl:2019:EMNLP1}.
A new insight provided by our study is that this difficulty holds even after alleviating the frequency effects by augmenting the target structures along with direct supervision signals.
This indicates that RNNs might inherently suffer from some memory limitation like a human subject, for which the difficulty of particular constructions, including center-embedded object relative clauses, are known to be incurred due to memory limitation \cite{gibson1998linguistic,demberg-keller-2008} rather than purely frequencies of the phenomena.
In terms of language acquisition, the supervision provided in our approach can be seen as direct negative evidence \cite{MARCUS199353}.
Since human learners are known to acquire syntax without such direct feedback we do not claim that our proposed learning method itself is cognitively plausible.

One limitation of our approach is that the scope of negative examples has to be predetermined and fixed.
Alleviating this restriction is an important future direction.
Though it is challenging, we believe that our final analysis for transferability, which indicates that the negative examples do not have to be complete and can be noisy, suggests a possibility of a mechanism to induce negative examples themselves during training, perhaps relying on other linguistic cues or external knowledge.

\section*{Acknowledgements}
We would like to thank Naho Orita and the members of Computational Psycholinguistics Tokyo for their valuable suggestions and comments.
This paper is based on results obtained from projects commissioned by the New Energy and Industrial Technology Development Organization (NEDO).


\section{Introduction}
Studies on photoionization processes and other photo-induced breakup reactions in atoms and molecules have always helped deepening our understanding of the structure and dynamics of matter on a microscopic scale. Due to well-defined transfers of energy and momentum during photoabsorption, such processes can provide detailed information on the underlying interaction mechanism as well as the involved atomic system. Accurate and comprehensive tests of corresponding theoretical models are nowadays enabled  by kinematically complete experiments \cite{Ullrich}.

Within the multitude of different photoionization processes, electron-electron correlations often play an integral part. A well known example is resonant photoionzation, where photoexcitation creates an autoionizing state which subsequently leads to ionization via Auger decay.
This mechanism can also be generalized to systems consisting of two (or more) atoms. In this case, the resonant excitation of one atom leads -- via
interatomic electron-electron correlations -- to a radiationless energy transfer to a neighbour atom, resulting in its ionization. This decay mechanism is a two-center version of the (usually intraatomic) Auger effect and also well-known as interatomic Coulombic decay (ICD) \cite{ICD,ICDres,ICDrev}. It has been experimentally observed in a variety of systems, e.g. in noble gas dimers \cite{dimers}, clusters \cite{ICDcluster} and water molecules \cite{ICDwater} after photoabsorption. ICD can also be triggered by electron impact which has been demonstrated in various noble-gas systems for a range of impact energies \cite{Grieves,Dorn, Lanzhou, Lanzhou2}.

Stabilization through ICD is included in resonant two-center photoionization (2CPI) which represents
the generalization of resonant photoionization to diatomic systems \cite{2CPI1,2CPI2}. The process describes the ionization of an
atom $A$ via photoexcitation of a neighbouring atom $B$ of different species. This way, an autoionizing two-center state is formed which can decay via a two-center Auger decay, that is through ICD. 2CPI has theoretically been studied in systems of two atoms whose internuclear separation was assumed to be fixed. Under suitable conditions it was shown to dominate over direct (i.e. single-center) photoionization
of atom $A$ by several orders of magnitude. As specific example, a diatomic system of Li and He was considered \cite{2CPI1,2CPI2}. Related theoretical studies (see also \cite{Perina}) have treated two-center single-photon
double ionization \cite{Alex}, strong-field multiphoton ionization \cite{Jacqui} and electron-impact ionization \cite{gruell}
in heteroatomic systems as well as two-photon single-ionization in homoatomic systems \cite{Mueller2011}.

Experiments on two-center processes are mostly performed on van-der-Waals molecules or clusters which
are characterized by small binding energies and large bond lengths. In particular, experimental studies on
2CPI were carried out on He-Ne dimers \cite{2CPIexp,Mhamdi18} and on Ne-Ar clusters \cite{Hergenhahn} using synchrotron radiation to induce
the 1s-3p transition in He at about 23 eV and the 2p-3s transition in Ne at about 17 eV, respectively. Both experiments found strong enhancements of the photoelectron yield as compared with the direct ionization channels of Ne or Ar, respectively.

In these experiments, the internuclear distance between
the two atoms is not fixed but varies due to the vibrational
molecular motion. Therefore, the original theory of 2CPI \cite{2CPI1}
needs to be further developed to account for the nuclear
motion as well. Of particular interest is the influence of the
latter on the total cross section of 2CPI and its relative
strength in comparison with the direct photoionization.
It is worth mentioning that angular distributions of electrons
emitted via 2CPI from He-Ne dimers have been calculated in
\cite{Mhamdi18}, including their vibrational level structure (for
a similar consideration of ICD in He$_2$, see \cite{Mhamdi20}).
Besides, we note that 2CPI has recently been studied in slow
atomic collisions where the internuclear distance is not fixed
either but subject to the relative atomic motion which was
assumed to follow a straight-line trajectory \cite{Voitkiv}.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.47\textwidth]{Fig1.eps}
\end{center}
\vspace{-0.5cm}
\caption{Scheme of resonant two-center photoionization 2CPI. A two-center autoionizing state is created via photoexcitation of atom $B$ (left). Subsequently, the latter transfers the excitation energy radiationlessly to atom $A$, leading to its ionization via two-center Auger decay (right).}
\label{figure1}
\end{figure}
In the present paper, we expand the theoretical description of 2CPI by incorporating the nuclear motion in
weakly bound heteroatomic dimers. This way,corresponding treatments for ICD \cite{Moiseyev, Jabbari} are extended
by fully including the photoexcitation step. While 2CPI in van-der-Waals molecules is a rather complex process,
we try to keep its treatment as analytical as possible in order to facilitate the physical
interpretation of the results. Our approach allows us to examine the changes, which originate from the
nuclear motion, to the predictions on 2CPI for fixed internuclear separations \cite{2CPI1,2CPI2}. Besides, a comparison
between 2CPI and the direct photoionization in a dimer is drawn and a formula for the ratio of the
corresponding cross sections is obtained.

Our paper is organized as follows: In Sec. II we present our theoretical considerations of 2CPI.
Applying the Born-Oppenheimer approximation, we first treat the (fast) electronic transitions
involved in 2CPI at fixed internuclear distance \cite{Fussnote2} and, afterwards, include the (slow) nuclear motion
in the potential formed by the Coulomb and van-der-Waals interactions between the atoms $A$ and $B$.
An expression for the 2CPI cross section in a weakly bound dimer is derived which contains transition
matrix elements between the relevant molecular vibrational states. In Sec. III we apply our theory to 2CPI
in Li-He dimers. The influence of the vibrational nuclear motion is demonstrated and its physical implications
are discussed. Concluding remarks are given in Sec. IV.

Atomic units (a.u.) will be used throughout unless
otherwise stated. The Bohr radius is denoted by $a_{0}$.

\section{Theory of resonant two-center photoionization}
\label{sec:theory}
As in \cite{2CPI1}, we describe the process of 2CPI in a simple system of two atomic centers $A$ and $B$. The geometry of the system, in which both atoms are initially in their ground states, is characterized by the linking vector ${\bf R}$ between the nuclei. Within the Born-Oppenheimer picture, electronic and nuclear processes can be treated separately. The fast electronic transitions can be calculated for fixed interatomic distances $R$, whereas the slow nuclear motion can be considered afterwards.
We note that $R$ has to be sufficiently large (covering at least several Bohr radii) in order to apply the approximations used in the following.
\subsection{Electronic transitions}
\label{sec:general}
To begin with, we calculate the process treating the two atoms individually and incorporate molecular interactions in the following sections.
The nuclei which are assumed to be at rest carry the corresponding charges $Z_{A}$ and $Z_{B}$.
We take the position of $Z_{A}$ as the origin and denote the coordinates of the  nucleus $Z_{B}$, the electron associated with $A$ and the electron associated with $B$ by $\bf R$, $\bf r$ and ${\bf r'}={\bf R}+{\boldsymbol \xi}$ respectively, where $\boldsymbol \xi$ is the position of the electron of atom $B$ relative to the nucleus $Z_{B}$.

The considered process involves two steps, the photoexcitation of atom $B$ and the subsequent radiationless energy transfer leading to the ionization of $A$. Therefore, we need to define initial, intermediate and final configurations of the two electrons associated with $A$ and $B$, which are illustrated in Fig. \ref{figure1}. Within this basic consideration, we restrict the expressions to one active electron on each center included in the transition. Therefore, we obtain the following states:
(I) The initial state $\Phi_{g,g}=\varphi_{g}({\bf r})\chi_{g}({\boldsymbol \xi})$, with total energy $E_{g,g}=\varepsilon_{g}+\epsilon_{g}$  has both electrons of $A$ and $B$ in their ground state. (II) In the intermediate state $\Phi_{g,e}=\varphi_{g}({\bf r})\chi_{e}(\boldsymbol \xi)$
with total energy $E_{g,e}=\varepsilon_{g}+\epsilon_{e}$, the electron of $B$ has been excited by the electromagnetic field and the electron in $A$ remains in its ground state. (III) The final state  $\Phi_{{\bf k},g}=\varphi_{{\bf k}}({\bf r})\chi_{g}(\boldsymbol \xi)$
with total energy $E_{{\bf k},g}=\varepsilon_{k}+\epsilon_{g}$ and $\varepsilon_{k}=\frac{k^{2}}{2}$  consists of the electron from $A$, which  has been emitted into the continuum with asymptotic momentum $\bf k$ and the electron of $B$, which has returned to its ground state.

In order to ionize atom $A$ in a two-center process including atom $B$, the ionization potential $I_{A}=|\varepsilon_{g}|$ has to be smaller than the energy difference $\omega_{B}=\epsilon_{e}-\epsilon_{g}$ of the electronic transition in atom $B$. The probability amplitude is calculated via the time-dependent perturbation theory,
\begin{eqnarray}
\label{eq: S2}
S^{(2)}&=&-\int\limits_{-\infty}^{\infty}dt\, \mathcal{V}^{AB}({\bf k},{\bf R})e^{-i(E_{g,e}-E_{k,g})t}\nonumber\\&&\times\int\limits_{-\infty}^{t}dt'\, \mathcal{W}^{B}(\omega)e^{-i(E_{g,g}+\omega-E_{g,e})t'}.
\end{eqnarray}
The matrix elements are given by
\begin{eqnarray}
\mathcal{V}^{AB}({\bf k},{\bf R})&=&\langle \Phi_{{\bf k},g}|\hat{V}_{AB}({\bf R})|\Phi_{g,e}\rangle\\
\mathcal{W}^{B}(\omega)&=&\langle\Phi_{g,e}|\hat{W}_{B}(\omega)|\Phi_{g,g}\rangle .
\end{eqnarray}
The operator $ \hat{W}_{B}e^{-i\omega t }$ induces the photoexcitation of atom $B$.
The external field interacting with the atomic centers is described as a classical field with linear polarization. Applying dipole approximation simplifies the employed vector potential:
\begin{equation}
 {\bf A}(t)= {\bf A}_{0}\cos(\omega t).
 \end{equation}
 The direction of polarization is set along the $z$-axis, ${\bf A}_{0}= A_{0}{\bf e}_{z}$. Hence, the interaction $\hat{W}_{\text{B}}$ reads:
  \begin{equation}
  \hat{W}_{B}=\frac{{\bf A}_{0}\cdot{\bf \hat{p}}_{{\bf \xi}}}{2c}.
  \end{equation}
 The interatomic interaction $\hat{V}_{AB}$ causes the energy transfer between $A$ and $B$ resulting in the ionization of $A$. The excited state in atom $B$ decays in a dipole-allowed transition. We neglect retardation effects and obtain
\begin{equation}
\label{eq: VAB}
\hat{V}_{AB}({\bf R})=\frac{{\bf r}\cdot{\boldsymbol \xi}}{R^{3}}-\frac{3({\bf r}\cdot{\bf R})({\boldsymbol \xi}\cdot{\bf R})}{R^{5}}.
\end{equation}
Integration of Eq. \eqref{eq: S2} leads to
\begin{eqnarray}
\label{eq: S2-delta}
S^{(2)}=-2\pi i \delta(E_{{\bf k},g}-E_{g,g}-\omega)\frac{\mathcal{V}^{AB}({\bf k},{\bf R})\mathcal{W}^{B}(\omega)}{\Delta+\frac{i}{2}\Gamma},
\end{eqnarray}
where $\Delta=\epsilon_{g}+\omega-\epsilon_{e}$ is the energy detuning. The resonance width $\Gamma$ accounts for the instability of the excited state of atom $B$ and includes the radiative width $\Gamma_{\text{rad}}$ and the two-center Auger width $\Gamma_{\text{aug}}$  \cite{Fussnote}
\begin{equation}
\Gamma=\Gamma_{\text{rad}}+\Gamma_{\text{aug}}.
\end{equation}
The radiative decay rate $\Gamma_{\text{rad}}$ describes the decay of the excited state of atom $B$ via spontaneous emission of a photon. The two-center Auger width $\Gamma_{\text{aug}}$ describes radiationless transfer of the transition energy in atom $B$ to atom $A$, leading to its ionization. Accordingly, the widths are obtained from the equations
 \begin{align}
  \label{eq: Gamma}
  \Gamma_{\text{rad}}&=\frac{4 \omega_{B}^{3}}{3 c^{3}}\left|\left\langle\chi_{g}\left|{\boldsymbol \xi}\right|\chi_{e}\right\rangle\right|^{2}\\
  \label{eq: Gammaaug}
  \Gamma_{\text{aug}}({\bf R})&=\int \frac{d^{3}k'}{(2\pi)^{2}}|\mathcal{V}^{AB}({\bf k'},{\bf R})|^{2}\delta(\varepsilon_{k}'+\epsilon_{g}-\varepsilon_{g}-\epsilon_{e}).
  \end{align}
  \subsection{Ratio of cross sections at fixed ${\bf R}$}
In order to compare the two-center photoionization, we also calculate the direct photoionization of atom $A$. The corresponding transition amplitude is given by
\begin{equation}
\label{eq: S1}
S^{(1)}=-i\int\limits_{-\infty}^{\infty}dt\,\langle \varphi_{{\bf k}}({\bm r})|\hat{W}_{A}(\omega)| \varphi_{g}({\bm r})\rangle e^{-i(\varepsilon_{g}+\omega-\varepsilon_{k})t}.
\end{equation}
Here, the operator $\hat{W}_{A}=\frac{{\bf A_{0}\cdot\hat{p}}_{r}}{2c}$ describes the coupling of the electron in atom $A$ to the electromagnetic field.
The direct photoionization process can interfere with the indirect process leading to the same final state. The resulting transition amplitude thus consists of two individual amplitudes
\begin{eqnarray}
 S^{(12)}=S^{(1)}+S^{(2)}.
 \end{eqnarray}
We then obtain the transition cross sections $\sigma^{(1)}$, $\sigma^{(2)}$ and $\sigma^{(12)}$ for the one-center, two-center and total photoionization processes, respectively; by virtue of
\begin{eqnarray}
\label{eq: sigmaallg}
\sigma^{(M)}=\frac{1}{\tau j}\int\frac{d^{3}k}{(2\pi)^{3}}|S^{(M)}|^{2},
\end{eqnarray}
where $\tau$ denotes the interaction time and $j=\frac{\omega A_{0}^{2}}{8\pi c}$ is the incident photon flux.
The cross sections for two-center and direct one-center photoionization can be related. Fixing ${\bf R}=R{\bf e}_{z}$, one obtains:
\begin{align}
\sigma^{(2)}&=\frac{1}{j}\frac{4}{R^{6}}\int \frac{d^{3}k}{(2\pi)^{3}}\left|\left\langle\varphi_{{\bf k}}\chi_{g}\left|z \xi_{z}\right|\varphi_{g}\chi_{e}\right\rangle\right|^{2}\times\nonumber\\
& 2\pi\delta(E_{{\bf k},g}-E_{g,g}-\omega)\frac{\left|\left\langle\chi_{e}\left|\frac{1}{2c}{\bf\hat{A}\cdot\hat{p}}_{\xi}\right|\chi_{g}\right\rangle\right|^{2}}{\Delta^{2}+\frac{1}{4}\Gamma^{2}}\nonumber\\
&\approx \sigma^{(1)}\frac{\Gamma_{\text{rad}}^{2}}{\Delta^{2}+\frac{1}{4}\Gamma^{2}}\left(\frac{3c^{3}}{4R^{3}\omega^{3}}\right)^{2},
\end{align}
On resonance, the ratio reads:
\begin{equation}
\label{eq: Abschaetzung_atom}
\frac{\sigma^{(2)}}{\sigma^{(1)}}\approx \frac{\Gamma_{\text{rad}}^{2}}{\Gamma^{2}}\left(\frac{3c^{3}}{2R^{3}\omega^{3}}\right)^{2}.
\end{equation}
$\Gamma_{\text{rad}}$ and $\Gamma_{\text{aug}}(R)$ can be calculated from Eqs. \eqref{eq: Gamma} and \eqref{eq: Gammaaug} or obtained from literature.
\subsection{Inclusion of nuclear motion}
\subsubsection{Energy shift}
\label{energy shift}
The presence of an atom in the proximity of another atomic center leads to an appearance of an interaction between them, described by a static potential. As a consequence, the energy of the system depends on their separation ${\bf R}$.
We therefore include interactions of Coulombic and exchanging nature as first order perturbations. Also van-der-Waals terms have to be included in the picture as they generally become more relevant than first order shifts for larger distances \cite{Patil}.
The Coulomb correction to the energy is given by
\begin{eqnarray}
V_{\text{coul}}(R)&=&\int d^{3}r d^{3}\xi \Big(\frac{\beta_{1}}{R}-\frac{\beta_{2}}{|{\bf R}-{\bf r}|}-\frac{\beta_{3}}{|{\bf R}+{\boldsymbol \xi}|}\nonumber\\
&&+\frac{\beta_{4}}{|{\bf R}-{\bf r}+{\boldsymbol \xi}|}\Big)|\Phi({\bf r},{\boldsymbol \xi})|^{2},
\end{eqnarray}
where the spatial dependencies are shown in Fig. \ref{figure2}.
Here, $\beta_{m}$ are coefficients containing the charges of particles involved in the respective interaction.\\
The exchange energy takes the indistinguishability of the electrons into account and is given by
\begin{align}
V_{\text{exc}}(R)&=-\frac{1}{2}\int d^{3}r d^{3}\xi \Big(\frac{\beta_{1}}{R}-\frac{\beta_{2}}{|{\bf R}-{\bf r}|}-\frac{\beta_{3}}{|{\bf R}+{\boldsymbol \xi}|}\nonumber\\
&+\frac{\beta_{4}}{|{\bf R}-{\bf r}+{\boldsymbol \xi}|}\Big)\Phi({\bf r},{\bf R}+{\boldsymbol \xi})\Phi^{*}({\bf R}+{\boldsymbol \xi},{\bf r}).
\end{align}
 The factor of $\frac{1}{2}$ is due to the spin. In general, the van-der-Waals interaction can be calculated as a second order correction of the form
\begin{eqnarray}
V_{\text{vdW}}(R)=\sum_{n'\neq m}\frac{|\langle \Phi_{m}|\hat{V}_{AB}|\Phi_{n'}\rangle |^{2}}{E_{n'}^{(0)}-E_{m}^{(0)}}.
\end{eqnarray}
However, when possible, we use literature values for better accuracy , where higher order terms ($\frac{C_{8}}{R^{8}}$,$\frac{C_{10}}{R^{10}}$) can be included. Consequently, the van-der-Waals interaction can be described by
\begin{equation}
V_{\text{vdW}}(R)=-f_{6}(R)\frac{C_{6}}{R^{6}}-f_{8}(R)\frac{C_{8}}{R^{8}}-f_{10}(R)\frac{C_{10}}{R^{10}}.
\end{equation}
Literature values require the restriction to one orientation of the internuclear vector ${\bf R}$.
The Van-der-Waals coefficients are multiplied by damping factors \cite{Patil}
\begin{equation}
 f_{2L+4}(R)=1-e^{-R/\overline{a}}\sum_{n=0}^{2L+7}\frac{1}{n\!}\left(\frac{R}{\overline{a}}\right)^{n}.
 \end{equation}
Here, $\overline{a}$ is defined by the ionization energies of both constituents $A$ and $B$:
 $\overline{a}=\frac{1}{4}\left(\frac{1}{\sqrt{2 \varepsilon_{g}}}+\frac{1}{\sqrt{2 \epsilon_{g}}}\right)$.
The damping functions are introduced to include the effects of charge overlap \cite{damping}.

\begin{figure}
\begin{center}
\includegraphics[width=0.47\textwidth]{Fig2.eps}
\end{center}
\caption{Scheme of the spatial dependencies considered in the calculation of the energy shift.}
\label{figure2}
\end{figure}
The total static interaction energy between the atoms is then given by the sum of these three interactions
\begin{equation}
V_{\text{int}}({\bf R})=V_{\text{coul}}+V_{\text{exc}}+V_{\text{vdW}}.
\end{equation}
\subsubsection{Vibrational States}
\label{sec:vibrational states}

The interaction of neighbouring atoms $A$ and $B$ can lead to a complex of these two constituents which cannot be treated individually anymore.
 This bound system is then characterized not only by its motion as a whole but also by its relative motion. The rotation and vibration of diatomic molecules  influences the binding energy on different energy scales. Since the rotational motion creates levels with a spacing much smaller than the levels due to vibrational motion, we will restrict our calculations to vibrational effects solely \cite{Bransden}.
We consider vibrational levels and wave functions for the three electronic states described in section \ref{sec:general}.

(I) $\Psi_{g,g}=\psi_{\text{i}}(R,\nu_{\text{i}})$ refers to the vibrational wave function for the  intial system including both atoms $A$ and $B$ in their electronic ground states. The corresponding vibrational energy reads $E_{\text{i}}(\nu_{\text{i}})=-D_{\text{i}}+E_{\text{vib}}(\nu_{\text{i}})$, where $D_{\text{i}}$ is the maximum depth of the corresponding interaction potential curve.
(II) $\Psi_{g,e}=\psi_{\text{a}}(R,\nu_{\text{a}})$ denotes the vibrational wave function for the intermediate system of an autoionizing state with $A$ in its electronic ground state and $B$ excited via photoexcitation. The vibration leads to the energy shift $E_{\text{a}}(\nu_{\text{a}})=-D_{\text{a}}+E_{\text{vib}}(\nu_{\text{a}})$.
(III) $\Psi_{{\bf k},g}=\psi_{\text{f}}(R,\nu_{\text{f}})$ describes the final  vibrational wave function  with $A$ being ionized and $B$ back in its electronic ground state. The vibration leads to the energy shift $E_{\rm{f}}(\nu_{\rm{f}})=-D_{\rm{f}}+E_{\rm{vib}}(\nu_{\rm{f}})$.\\
We calculate the vibrational wave functions and the energy levels from the potential energy shift.
In order to avoid a numerical computation of the vibrational energies and states, we fit a Morse potential to our potential curves in order to find the necessary parameters. The potential reads
\begin{eqnarray}
\label{eq:Morse}
V_{\rm{Morse}}(R)=D\left[\left(1-e^{-\alpha (R-R_{\rm{eq}})}\right)^{2}-1\right],
\end{eqnarray}
where $D$ is the depth of the potential, $R_{\text{eq}}$ is the equillibrium distance and $\alpha$ describes the width of the potential.
Employing this fit, we obtain the vibrational energy levels
\begin{equation}
\label{eq:Morse_Energie}
E_{\rm{vib}}(\nu)=\sqrt{\frac{2D\alpha^{2}}{\mu}}\left[\left(\nu+\frac{1}{2}\right)-\frac{1}{\kappa}\left(\nu+\frac{1}{2}\right)^{2}\right]
 \end{equation}
 with $\kappa=\sqrt{\frac{8\mu D}{\alpha^{2}}}$ and the reduced mass of the system $\mu$.
The parameters obtained by the fit are then inserted in the vibrational wave function \cite{Gallas}:
\begin{eqnarray}
\psi(R,\nu)=\sqrt{\frac{\alpha b (\nu!)}{\Gamma(\kappa-\nu)}}e^{-z/2}z^{(\kappa-2\nu-1)/2}\mathcal{L}_{\nu}^{(\kappa-2\nu-1)}(z)
\end{eqnarray}
with the associated Laguerre polynomials $\mathcal{L}_{\nu}^{(\kappa-2\nu-1)}(z)$, $b=\kappa-2\nu-1$ and $z=\kappa e^{-\alpha (R-R_{\rm{eq}})}$.

The vibrational wave functions contribute to the cross section by the insertion of Franck-Condon factors, which describe the overlap of the vibrational nuclear wave functions  associated with the two electronic states involved in a transition. The factor therefore favours transitions with large overlaps, involving no substantial changes of the nuclear motion in terms of position and kinetic energy \cite{LaForge}. Photoionization, photoexcitation and the energy transfer via dipole-dipole interaction do not limit the possible combinations of vibrational states. Every matrix element describing an electronic transition contains an integration over $R$ as well. The general form of a transition matrix element is given by:
	    \begin{equation}
	 \mathcal{M}=   \langle \Phi_{2}^{el}({\bf r},{\boldsymbol \xi})\Psi_{2}(R,\nu_{2})|\hat{V}({\bf R},{\bf r},{\boldsymbol \xi})|\Phi_{1}^{el}({\bf r},{\boldsymbol \xi})\Psi_{1}(R,\nu_{1})\rangle,
	    \end{equation}
	     where $\hat{V}$ is an arbitrary operator, $i$ and $f$ refer to the inital and final states of the transition. The matrix element depends on the vibrational levels $\nu$.

We do not apply the Condon approximation in our calculation since $\hat{V}_{AB}$ depends heavily on $R$. Therefore, the corresponding Franck-Condon factor cannot be separated from the matrix element consisting of electronic and vibrational transitions.
The vibrational levels $\nu$ for the intermediate and final state  give rise to a multitude of transitions. Hence, we have to sum coherently over all possible intermediate states in the transition amplitude $S^{(2)}_{\text{mol}}$ and incoherently over all possible final states within the associated two-center cross section.

Considering the vibrational wave functions, the previously position-dependent  $\Gamma_{\text{aug}}(R)$ is transformed into $\bar{\Gamma}_{\text{aug}}$. For this purpose we average $\Gamma_{\text{aug}}(R)$ over the probability density $\left|\psi_{a}(R,\nu_{a})\right|^{2}$ of the intermediate vibrational state for every vibrational level $\nu_{a}$ \cite{Fussnote3}. For the radiative width the atomic value from Eq. \eqref{eq: Gamma} is inserted because $\Gamma_{\text{rad}}$ remains practically unaltered in the presence of the neighbour atom.

Including vibrational levels, vibrational wave functions and decay rates, the explicit expressions accounting for the molecular effects read
\begin{widetext}
\begin{align}
S^{(2)}_{\text{mol}}=-2\pi i\sum_{\nu_{\text{a}}}\delta(E_{{\bf k},g}+E_{\text{f}}(\nu_{\text{f}})-E_{g,g}-E_{\text{i}}(\nu_{\text{i}})-\omega)\frac{\langle \Phi_{{\bf k},g}\Psi_{{\bf k},g}|\hat{V}_{AB}({\bf R})|\Phi_{g,e}\Psi_{g,e}\rangle \langle\Phi_{g,e}\Psi_{g,e}|\hat{W}_{B}(\omega)|\Phi_{g,g}\Psi_{g,g}\rangle}{\epsilon_{g}+E_{\text{i}}(\nu_{\text{i}})+\omega-\epsilon_{e}-E_{\text{a}}(\nu_{\text{a}})+\frac{i}{2}(\Gamma_{\text{rad}}+\bar{\Gamma}_{\text{aug}})}
\end{align}
\begin{align}
S^{(1)}_{\text{mol}}=-2\pi i \delta(E_{{\bf k},g}+E_{\text{f}}(\nu_{\text{f}})-E_{g,g}-E_{\text{i}}(\nu_{\text{i}})-\omega) \langle\Phi_{{\bf k},g}\Psi_{{\bf k},g}|\hat{W}_{A}(\omega)|\Phi_{g,g}\Psi_{g,g}\rangle
\end{align}
\end{widetext}
 Since $\hat{W}_{A}$ and $\hat{W}_{B}$ do not depend on $R$, the corresponding Franck-Condon factors can be separated from the electronic transition matrix elements. Therefore,
\begin{align}
F_{\text{i},\text{a}}(\nu_{\text{i}},\nu_{\text{a}})&=\int dR \, \psi_{\text{i}}(R,\nu_{\text{i}})\psi^{*}_{\text{a}}(R,\nu_{\text{a}})\\
F_{\text{i},\text{f}}(\nu_{\text{i}},\nu_{\text{f}})&=\int dR \, \psi_{\text{i}}(R,\nu_{\text{i}})\psi^{*}_{\text{f}}(R,\nu_{\text{f}}).
\end{align}
For $F_{\text{a},\text{f}}(\nu_{a},\nu_{f})$ however, one has to incorporate the $R$-dependency of $\hat{V}_{AB}$. The factor in question can be simplified to $1/R^{3}$, see Eq. \eqref{eq: VAB}. This yields
\begin{equation}
F_{\text{a},\text{f}}(\nu_{\text{a}},\nu_{\text{f}})=\int dR\,\frac{ \psi_{\text{a}}(R,\nu_{\text{a}})\psi^{*}_{\text{f}}(R,\nu_{\text{f}})}{R^{3}}.\\
\end{equation}
\section{Results and Discussion}
\label{Results}
Let us apply the theory established in Sec. \ref{sec:theory} to the photoionization of a weakly bound van-der-Waals molecule, which is chosen as a $^{7}$Li$^{4}$He dimer. The internuclear axis is taken along the $z$-axis, which also serves as the quantization axis.
We choose effective electronic states $\chi$ and $\varphi$ in order to achieve a reasonable comparableness to the atomic species under consideration.
Neither Li nor He are single-electron atoms. To describe He as a two-electron system we employ symmetrized
superpositions of product states formed by hydrogen-like wave functions
with an effective charge $Z_{\rm He}$. The He ground state
correspondingly reads $\chi_{1s}({\boldsymbol \xi}_{1}) \chi_{1s}({\boldsymbol \xi}_2)$,
and the excited state is given by $\frac{1}{\sqrt{2}}[\chi_{1s}({\boldsymbol \xi}_{1})
\chi_{2p_{0}}({\boldsymbol \xi}_{2})+\chi_{2p_{0}}({\boldsymbol \xi}_{1}) \chi_{1s}({\boldsymbol \xi}_{2})]$. Within such an approach one can calculate the matrix elements including helium by using only a one-particle wave function $\chi({\bf \xi})$ and then multiplying by  a factor of 2 contributing for the superposition. Note, that this also applies to the two-center Auger decay width $\Gamma_{\text{aug}}$ for both the atomic and molecular consideration.
For the valence electron of the alkali atom lithium in its ground state, we choose the Bates-Darngaard wave function as in \cite{Patil}, so that
\begin{eqnarray}
\varphi_{g}({\bf r})=\frac{1}{\sqrt{8\pi}}\frac{1}{\Gamma(a+1)}\left(\frac{2}{a}\right)^{a+\frac{1}{2}}r^{a-1}\left(1+\frac{v}{r}\right)e^{-r/a}
\end{eqnarray}
 with $a=\frac{1}{\sqrt{2|\varepsilon_{g}|}}$ depending on the binding energy and $v=-\frac{1}{2}a^{2}(a-1)$. Note that $a$ corresponds to an effective nuclear charge $Z_{\text{Li}}\approx 1.259$ for $\varepsilon_{g}= 5.39\text{eV}$ \cite{NIST}.

After ionization, the influence of the remaining lithium ion on the emitted electron of Li is accounted for by employing a Coulomb wave $\varphi_{{\bf k}}$ \cite{LandauQM}. For this state, the same effective nuclear charge $Z_{\text{Li}}$ as in the ground state is chosen.
In accordance with Eq. \eqref{eq: sigmaallg}, this continuum state is normalized to a quantization volume of unity.\\
Within this system, we consider the dipole-allowed transition $(1s-2p_{0})$ for the photo-excitation of atom $B$. The transition energy defines the effective nuclear charge for both states. With an excitation energy $\omega_{\text{He}}=\text{21.218 \text{eV}}$ \cite{NIST}, we find $Z_{\text{He}}=1.435$ .
We have now captured some basic features of the Li-He dimer.

\subsubsection{Results at fixed nuclei}
Before incorporating molecular effects into our description of the system  \cite{Friedrich}, let us first calculate the cross sections for the ionization of lithium by disregarding the molecular structure and rather considering Li and He as two unbound atoms which are separated by a fixed internuclear distance $R$ on which the two-center cross section strongly depends \cite{2CPI1,2CPI2}.

Due to the atomic approach, the two-center Auger decay width $\Gamma_{\text{aug}}$ in Eq. \eqref{eq: Gammaaug} is very sensitive to $R$. Both the decay width $\Gamma_{\text{aug}}$ and the square of the matrix element of the radiationless energy transfer include the dependency on $\frac{1}{R^{6}}$. On resonance, the $R$-dependence of the ratio $\eta_{\text{atomic}}^{(2)}(\omega,R)=\sigma^{(2)}_{\text{atomic}}/\sigma^{(1)}_{\text{atomic}}$ is depicted in Fig. \ref{fig:ohne_Kernbewegung} (a). For relatively small distances $R$, the ratio of cross sections shows a  steep incline for increasing  spacings. For sufficiently large separations, $\Gamma_{\text{aug}}$ becomes small in comparison to the constant radiative decay width  $\Gamma_{\text{rad}}$. Consequently, the saturation of the total decay width $\Gamma$ sets in quickly. The dependence on $R$ of the matrix element mentioned above, however, does not exhibit saturation. Therefore, the cross section ratio depicted in Fig.\ref{fig:ohne_Kernbewegung} (a) shows the mentioned $\frac{1}{R^{6}}$ behaviour after the maximum is reached.

For fixed internuclear separations $R=20\,a_{0}$ and $R=30\,a_{0}$, the ratio of cross sections is plotted against the incident photon energy $\omega$ in Fig. \ref{fig:ohne_Kernbewegung} (b). The plot shows a single peak where the photon energy is resonant with the electronic $1s$-$2p$ transition in He. The peak height demonstrates that 2CPI can strongly dominate over the direct photoionization.
\begin{figure}[t]
\includegraphics[width=0.49\textwidth]{Fig3.eps}
\caption{Ratios $\eta^{(2)}_{\text{atomic}}(\omega,R)=\frac{\sigma^{(2)}_{\text{atomic}}}{\sigma^{(1)}_{\text{atomic}}}$ in a Li-He system, ignoring molecular effects. (a) $\eta^{(2)}_{\text{atomic}}(\omega,R)$ plotted against the internuclear distance $R$ on resonance. (b) $\eta^{(2)}_{\text{atomic}}(\omega,R)$ plotted against the incident photon energy $\omega$ for fixed internuclear distance $R=20a_{0}$ (solid) and $R=30a_{0}$ (dashed). }
\label{fig:ohne_Kernbewegung}
\end{figure}
\subsubsection{Effects of nuclear motion}
Until now, we have not taken any effects of the molecular bond in Li-He into account.
In neighbourhood of each other, the influence of the two atoms leads to a potential surface creating a possibly bound state.
We now consider these potential energy curves for the three electronic states participating in the process depending on the internuclear distance $R$.

Note that for the initial state, where both atoms are in their ground state, the direction of $R$ does not play a role. For the intermediate state, however, where helium is excited, the interaction depends on the direction of ${\bf R}$. Here, ${\bf R}=R{\bf e}_{z}$ is considered, as mentioned above.

For the chosen geometry, the van-der-Waals coefficients can be given for inital, intermediate and final electronic state of Li-He.
(I): For the initial state, we find $C_{6}=22.5 \text{ a.u.}$, $C_{8}=1.06\times10^{3}\text{ a.u.}$ \cite{Patil}.
(II): The intermediate state is described using, $C_{6}=6123\text{ a.u.}$, $C_{8}=7.85\times10^{5}\text{ a.u.}$, $C_{10}=1.02\times10^{8}\text{ a.u.}$ \cite{vdW_2p}.
(III): For the ionic system of the final state, $C_{6}=0.298\text{ a.u.}$, $C_{8}=1.98\text{ a.u.}$ are employed. This way, we obtain the relevant terms as calculated in \cite{Soldan}.

\begin{figure}
\includegraphics[width=0.49\textwidth]{Fig4a.eps}
\caption{Potential energy curves (solid) of the electronic states involved in the two-center photoionization process of Li-He and some corresponding vibrational levels (dashed). (a) Initial state.  The only supported vibrational level $\nu_{\text{\text{i}}}=0$ exhibits a binding energy  that is too small to be graphically resolved here. (b) Intermediate state with the lowest-lying vibrational levels.
In total, 11 vibrational levels ($0\le\nu_a\le 10$) associated with bound
molecular states were obtained; all of them are included in our numerical
calculations.
(c) Final state including 6  bound vibrational levels ($0\le\nu_{\text{f}}\le 5$). As an absolute scale, the potential energy $E$ relative to the ground state energy at infinite internuclear distance is depicted on the right side. }
\label{fig:energyshift}
\end{figure}

In Fig. \ref{fig:energyshift}, the interaction potential curves of the relevant electronic states of Li-He are depicted. For the initial state $\Phi_{g,g}$ and the intermediate state $\Phi_{g,e}$ the potential minima are close to each other. The final state $\Phi_{{\bf k},g}$, however, features its minimum at a much smaller distance. The depths of the potential for the final and intermediate state are comparable in magnitude.

The initial state, however, exhibits a much more shallow potential curve.
This leads to the initial state supporting only one vibrational level. The binding energy of merely $6\text{mK}=1.9\times 10^{-8}\text{a.u.}$ \cite{Friedrich} is not depicted graphically here since it is too close to the plot axis. The vibrational level expands widely and
spans over a large range of interatomic distances $R$ allowing for a very large bond length. In the results from  \cite{Friedrich}, the equillibrium distance $R_{\text{eq}}\approx 11.3 a_{0}$, representing the minimum of the potential, is much smaller than the mean distance $R \approx 55a_{0}$.
Literature values for the potential curves can be found for the initial and final state \cite{Patil, Soldan}and show good agreement with our calculations (within about $10\%$). For the intermediate states, no literature values beyond the van-der-Waals coefficients were available.

From the potential curves depicted in Fig.\ref{fig:energyshift}, the parameters for a Morse potential [see Eq. \eqref{eq:Morse} and Table \ref{tab: Morse}] are extracted.
\begin{table}[h]
\begin{center}
\begin{tabular}{|c||c|c|c| }
&ground state & intermediate state & final state\\
 \hline
$R_{\text{eq}}$&11.9&11.0&3.5\\
$D$& $5.7\times10^{-6}$&$2.9\times 10^{-3}$&$3.5\times 10^{-3}$ \\
$\alpha$&0.43&0.44&0.80
\end{tabular}
 \caption{Fitted Morse parameters (in a.u.) for the interaction potential curves. }\label{tab: Morse}
\end{center}\end{table}
 The resulting vibrational energy levels are calculated  according to Eq. \eqref{eq:Morse_Energie} and also shown in the figure. The inclusion of molecular effects leads to a shift concerning the energy of the Li-He system  for every electronic state. Compared to the atomic consideration, the resonance energies of the system are blueshifted \cite{Hergenhahn}.

The vibrational wave functions are calculated from the obtained Morse potential as described in Sec. \ref{sec:vibrational states}.
Since the positions of the potential minima for the initial and intermediate state are quite close to each other, their corresponding vibrational wave functions yield a significant overlap for all intermediate vibrational levels $\nu_{\text{\text{a}}}$. The potential minimum of the final state, however, is strongly shifted to smaller interatomic distances $R$. An overlap with the other sets of states $\lbrace \psi_{\text{a}}(R,\nu_{\text{a}})\rbrace$, $\lbrace \psi_{\text{i}}(R,\nu_{\text{i}})\rbrace$ is only present for high vibrational levels $\nu_{\text{\text{f}}}$. This results in a strong dependence of the Franck-Condon factors on the combination of vibrational levels.

For the calculation including molecular effects, the decay widths have to be evaluated. Using the wave functions mentioned in Sec. \ref{sec:general}, $\Gamma_{\text{rad}}^{\text{He}}=6.62\times10^{-8}\,a.u.$
 is obtained. This value differs from the literature value $\Gamma_{\text{rad}}^{\text{He}}=4.35\times 10^{-8}\,a.u.$ \cite{NIST} by a factor of about $1.5$. However, since the wave functions are inserted in all matrix elements of both processes, we use our value for selfconsistency. The two-center Auger decay rate $\Gamma_{\text{aug}}(R)$  in Eq. \eqref{eq: Gammaaug} depends on the internuclear distance $R$ and therefore is influenced by the nuclear motion. The two-center Auger decay is only allowed for distances greater than $\approx 4.3a_{0}$, where our potential curves for intermediate and final states cross \cite{Scheit}.

Before showing our numerical results, the analytical approach provides the possibility to obtain an approximate formula displaying the influence of the nuclear motion on the ratio of cross sections $\eta^{(2)}(\omega)=\frac{\sigma^{(2)}(\omega)}{\sigma^{(1)}(\omega)}$. With the restriction to one combination of vibrational levels $\nu_{\text{a}}$ and $\nu_{\text{f}}$, an approximation of $\eta^{(2)}$ can be given:
\begin{equation}
\label{eq: FC_Vereinfachung}
\frac{\sigma^{(2)}}{\sigma^{(1)}}\Bigg|_{\omega=\omega_{\text{res}}}\approx \frac{9c^{6}}{\omega_{\text{res}}^{6}R_{\text{eq}}^{6}}\frac{\Gamma_{\text{rad}}^{2}}{\Gamma^{2}}\left(\frac{R_{\text{eq}}^{6}\left(F_{\text{i,a}}F_{\text{a,f}}\right)^{2}}{\left(F_{\text{i,f}}\right)^{2}}\right),
\end{equation}
where $\Gamma$ already includes the two-center Auger width $\bar{\Gamma}_{\text{aug}}$ incorporating the nuclear motion. Note, that the resonant energy  $\omega_{\text{res}}$ depends on the vibrational energy levels involved. Ignoring molecular effects, $\omega_{\text{res}}=\omega_{B}$ meets the criterion for  resonance.

A comparison of Eq. \eqref{eq: FC_Vereinfachung} with Eq. \eqref{eq: Abschaetzung_atom}, which was obtained in an atomic picture, reveals the influence of molecular effects on 2CPI. Considering the photoionization processes including molecular effects, the   formula in Eq. \eqref{eq: Abschaetzung_atom} is modified by multiplying the term in parentheses including the Franck-Condon overlaps. $R_{\text{eq}}$ has been inserted in the parentheses of Eq. \eqref{eq: FC_Vereinfachung} to make the ratio of Franck-Condon factors dimensionless, since $F_{\text{a,f}}$ includes the $R$-dependence of the electron-electron interaction. Note, that the ratios  for Li-He calculated from an atomic point of view are still dependent on $R$. Besides, the expression in Eq. \eqref{eq: Abschaetzung_atom} differs by a factor of $4$ due to the fact that here, we consider helium as a two-electron system.

It should be mentioned that the full calculation of the 2CPI cross section including vibrational effects calls for the summation over $\nu_{\text{a}}$ of the transition amplitudes as well as for the summation over $\nu_{\text{f}}$ in the cross section. Therefore, Eq. \eqref{eq: FC_Vereinfachung} can only indicate estimated values for particular vibrational transitions. As we will see below, for the strongest transitions $\nu_{\rm{f}}=5$ the values computed from Eq. \eqref{eq: FC_Vereinfachung} are only slightly reduced as compared with $\eta_{\text{atomic}}^{(2)}(\omega_{B},R)$. The approximated values are depicted by crosses for individual $\nu_{\rm{a}}$ in Fig. \ref{fig:P2/P1} a).

With these preparations we now proceed to the numerical computation of the ratios $\eta^{(2)}(\omega)=\frac{\sigma^{(2)}(\omega)}{\sigma^{(1)}(\omega)}$ and $\eta^{(12)}(\omega)=\frac{\sigma^{(12)}(\omega)}{\sigma^{(1)}(\omega)}$. We include the vibrational levels $\nu_{\text{a}}=0,1,...,10$ and $\nu_{\text{f}}=0,1,..,5$.
While the direct photoionization process exhibits a cross section smoothly depending on the incident photon energy $\omega$, the two-center process is a  resonant one \cite{2CPI1,2CPI2,2CPIexp}. Including the various combinations of vibrational levels, the single atomic resonance peak  of Fig. \ref{fig:ohne_Kernbewegung} b) is expected to fan out with respect to the photon energy $\omega$.
\begin{figure}
\includegraphics[width=0.49\textwidth]{Fig5.eps}\\
\caption{Ratios of cross sections in a Li-He dimer including the vibrational levels $\nu_{\text{a}}\in\left\lbrace 0,1,2,..,10\right\rbrace $ and $\nu_{\text{f}}\in\left\lbrace 0,1,2,3,4,5\right\rbrace $. a) $\eta^{(2)}(\omega)$ and values from the approximate expression in Eq \eqref{eq: FC_Vereinfachung} (crosses). b) $\eta^{(12)}(\omega)$ containing interference of direct and two-center photoionization}
\label{fig:P2/P1}
\end{figure}
The ratios depicted in Fig. \ref{fig:P2/P1} show a multiplet of peaks caused by the various vibrational levels of the intermediate molecular state. They are responsible for the shift of the particular resonant energy. The vibrational level of the final state does not affect the resonance condition and therefore causes no additional splitting. The peak on the most lefthand side can be associated with  vibrational level $\nu_{\text{a}}=0$ and $\nu_{\text{a}}$ is ascending from left to right for the other peaks.

The ratio $\eta^{(12)}$ considering the total cross section in Fig. \ref{fig:P2/P1} b) shows a sequence of Fano profiles (alternating  in phase between adjacent vibrational levels $\nu_{\text{a}}$ and $\nu_{\text{a}}+1$). The cross section $\sigma^{(12)}$ includes both the direct and indirect photoionization processes which allows for interference and therefore leads to the Fano profiles including minima and maxima.
However, the Fano profile of a given resonance is strongly modified due to the summation over the final vibrational level $\nu_{\text{f}}$. The characteristics of the minimum differ for each level $\nu_{\text{f}}$. Considering $\eta^{(12)}$ for each $\nu_{\rm{f}}$ individually, the peak position remains the same while the position and depth of the minimum changes. As a result, some Fano minima are cancelled by the overlap of all $\nu_{\text{f}}$. A relatively deep minimum can only be seen in between $\nu_{\text{a}}=0$ and $\nu_{\text{a}}=1$ and next to $\nu_{\text{f}}=10$, where the individual components overlap in order to form a minimum.
This modification represents a strong contrast to the atomic calculation \cite{2CPI1}, where the Fano profile is much more pronounced.
Off resonance, the ratio eventually falls off to $1$. The ratio $\eta^{(12)}$ shown in Fig. 5b) involves incoherent summations
over the final vibrational states, separately in the numerator and the
denominator. In order to obtain a pure Fano profile, a fixed value of
$\nu_{\text{f}}$ must be considered instead. For $\nu_{\text{f}}=5$, representing the strongest
contribution to a peak, the Fano parameters q are given in Table \ref{tab:Fano}.
While the sign of the Fano parameters alternates, their absolute value,
roughly speaking, increases from small to large $nu_{\text{a}}$. Being on the order
of $10^{3}$, they resemble the corresponding values obtained for fixed nuclei
(see caption of Table \ref{tab:Fano}).

\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c||c|c||c|c| }
$\nu_{\text{a}}$& $q$&$\nu_{\text{a}}$&$q$&$\nu_{\text{a}}$& $q$\\
 \hline
0& $1.0\times 10^{3}$&4 & $1.7\times 10^{3}$&8 & $3.3\times 10^{3}$    \\
1 & $-5.6\times 10^{2}$ &5 & $-1.7\times 10^{3}$ & 9 & $-3.8\times 10^{3}$   \\
2 & $1.2\times 10^{3}$ &6 & $2.3\times 10^{3}$   &10 & $5.0\times 10^{3}$  \\
3 & $-1.1\times 10^{3}$  &7 & $-2.5\times 10^{3}$&&  \\
\end{tabular}
 \caption{Approximated Fano parameters $q$ for each vibrational state $\nu_{\text{a}}$ for $\nu_{\text{f}}=5$. For fixed nuclei and $R=10 \text{a.u.}$ ($R=20 \text{a.u.}$), one finds $q \approx 3.2\times 10^{3}$ ($q \approx 3.6\times 10^{3}$).}\label{tab:Fano}
\end{center}\end{table}

\begin{figure}[h]
\includegraphics[width=0.49\textwidth]{Fig6.eps}
\caption{Ratio  $\bar{\eta}^{(2)}=\frac{\bar{\sigma}^{(2)}}{\bar{\sigma}^{(1)}}$ for a Li-He dimer.}
\label{fig:Mittelung}
\end{figure}

The orders of magnitude of the peak heights shown  in Fig. \ref{fig:P2/P1}  a) are comparable to those in Fig \ref{fig:ohne_Kernbewegung} b) and Ref. \cite{2CPI1}. Therefore, the inclusion of molecular effects does not change the magnitude of the resonant enhancement significantly. We note that from left to right, i.e. from
$\nu_{\text{a}}=0$ to $\nu_{\text{a}}=10$, the peaks become more narrow
and have a tendency to reach larger maximum values.

In order to further compare the calculations for fixed nuclei with the calculations involving molecular effects, the total decay widths of the peaks can be considered in table \ref{tab: 1}. They have been obtained by analyzing the data from the numerical calculations. We see that the total width decreases with the increase of the vibrational level of the intermediate state $\nu_{\text{a}}$. This is because the internuclear distance tends to grow for higher lying vibrational excitations, as illustrated in Fig. \ref{fig:energyshift}. The same trend of decreasing ICD widths has also been observed in He-Ne dimers \cite{2CPIexp,Mhamdi18}.
\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c||c|c||c|c| }
$\nu_{\text{a}}$& \textbf{$\Gamma$}(a.u.)&$\nu_{\text{a}}$& \textbf{$\Gamma$}(a.u.)&$\nu_{\text{a}}$& \textbf{$\Gamma$}(a.u.)\\
 \hline
0& $5.3\times 10^{-7}$&4 & $3.8\times 10^{-7}$&8 & $2.2\times 10^{-7}$    \\
1 & $5.0\times 10^{-7}$ &5 & $3.4\times 10^{-7}$ & 9 & $1.8\times 10^{-7}$   \\
2 & $4.6\times 10^{-7}$ &6 & $3.0\times 10^{-7}$   &10 & $1.4\times 10^{-7}$  \\
3 & $4.2\times 10^{-7}$  &7 & $2.6\times 10^{-7}$&&  \\
\end{tabular}
 \caption{Total decay widths of the intermediate states in a Li-He dimer,
depending on the vibrational level.
For fixed nuclei the decay width amounts to
$\Gamma = 6.91\times 10^{-7}$\,a.u., ($\Gamma = 7.59\times 10^{-8}$\,a.u.) at $R=10$\,a.u. ($R=20$\,a.u.).}\label{tab: 1}
\end{center}\end{table}

Another relevant quantity is the integral resonance strength. It can be roughly estimated by the product of peak height and width.
The corresponding value for $\nu_{\text{a}} = 0$ approximately is $\sigma^{(2)}\Gamma\big|_{\nu_{a}=0}\approx 1.34\times10^{-7}\text{ a.u.}$., for example,
whereas the {\it relative} integral resonance strength amounts to
$\eta^{(2)}\Gamma\vert_{\nu_{\text{a}}=0}\approx 0.65 $. Summing the relative integral resonance strength over all peaks, we obtain a total value of about $15 \text{a.u.}.$
It is comparable to values of $\approx 12 \text{a.u.}$ which are obtained in the atomic calculations $\eta_{\text{atomic}}(\omega_{B},R)*\Gamma(R)$ for a nuclear distance around $R_{\text{eq}}$. Hence, loosely speaking, the atomic resonance height is redistributed over the vibrational multiplet of resonance lines in the molecular case. This interpretation fits to our previous approximation in Eq. \eqref{eq: FC_Vereinfachung}.

In an experimental setup, the frequency of the external field is not exactly defined, but exhibits a width $\Delta\omega$. In order to account for this, we perform a corresponding averaging over all frequencies  by employing a Gaussian distribution with $\varsigma=\frac{\Delta \omega}{2\sqrt{2\ln{2}}}$. The averaged cross section reads
\begin{equation}
\bar{\sigma}^{(M)}(\omega)=\int d\omega' \sigma^{(M)}(\omega')\frac{1}{\sqrt{2\pi\varsigma^{2}}}e^{-\frac{(\omega'-\omega)^{2}}{2\varsigma^{2}}}.
\end{equation}
Applying an  FWHM energy width $\Delta \omega=1.7\text{meV}$ of the incident photon beam as used in  \cite{2CPIexp} we obtain averaged ratios depicted in Fig. \ref{fig:Mittelung}.
Due to the finite width of the field frequency, the peaks are reduced in height but also broadened. As a consequence, the integral resonance strength remains comparable to the one deduced from Fig. \ref{fig:P2/P1}. For example. when calculating the value again for $\nu_{\text{a}}=0$, we obtain the same order of magnitude as given above.

Before proceeding to the conclusions, we note that
the inclusion of molecular rotations, which have been disregarded
in our treatment, is expected to lead to an additional fine structure
of the resonance lines in Fig. \ref{fig:P2/P1}. Each of them would be split into a multiplet of lines which are associated with different rovibrational
transitions from the ground to the intermediate state \cite{Fussnote_rot}. The resulting
line splitting within such a multiplet would be very narrow,
since the energy differences between rotational states belonging to
the same vibrational level is very small. An experimental resolution
of this substructure would therefore be challenging.

 \section{Conclusion and Outlook}
Resonant two-center photoionization in weakly bound systems has been studied, including the nuclear motion due to molecular vibrations. Earlier treatments of the process assumed spatially fixed nuclei and obtained, accordingly, a 'local' 2CPI cross section which strongly depends on the internuclear distance $R$ and exhibits a single resonance peak.
When treating the system in a molecular way instead, the $R$-dependent interaction potential is incorporated in a 'global' way by defining the vibrational states of the nuclear motion.

An approximate expression highlighting the molecular effects
was established (see  Eq. \eqref{eq: FC_Vereinfachung}) which shows that the changes
to the calculation assuming fixed $R$ can be attributed to
Franck-Condon overlaps and vibrational energy levels. As
a consequence, the two-center cross section splits up into
a multitude of peaks at various resonance energies which are
determined by the electronic and vibrational transition energies.
In the studied case of Li-He dimers, the cross section of each peak
is reduced as compared to the value obtained for fixed (equilibrium)
internuclear distance. Summing over all peaks, however,
yields a total value similar to the `local' resonant cross sections
which are obtained when a fixed internuclear separation $R$
close to the equilibrium distance is taken.

The latter observation indicates that theoretical treatments of
other two-center processes assuming fixed internuclear distances
can be expected to provide meaningful predictions for weakly bound molecular
systems, as well. Corresponding recent studies have considered,
for example, electron-impact ionization \cite{gruell} and photo double ionization \cite{Alex},
resonant electron scattering and recombination \cite{Eckey_Jacob}, and electron capture \cite{Sisourat}.
The influence of nuclear motion in a van-der-Waals molecule on (some of)
these two-center processes will be the subject of a forthcoming study.

\section*{Acknowledgement}
This work has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Grant No. 349581371 (MU 3149/4-1 and VO 1278/4-1)