# Lower Perplexity is Not Always Human-Like

Tatsuki Kuribayashi1,2, Yohei Oseki3,4, Takumi Ito1,2,
Ryo Yoshida3, Masayuki Asahara5, Kentaro Inui1,4
1Tohoku University 2Langsmith Inc. 3University of Tokyo 4RIKEN 5NINJAL
{kuribayashi, takumi.ito.c4<EMAIL_ADDRESS>,
{oseki<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In computational psycholinguistics, various language models have been
evaluated against human reading behavior (e.g., eye movement) to build human-
like computational models. However, most previous efforts have focused almost
exclusively on English, despite the recent trend towards linguistic universal
within the general community. In order to fill the gap, this paper
investigates whether the established results in computational
psycholinguistics can be generalized across languages. Specifically, we re-
examine an established generalization —the lower perplexity a language model
has, the more human-like the language model is— in Japanese with typologically
different structures from English. Our experiments demonstrate that this
established generalization exhibits a surprising lack of universality; namely,
lower perplexity is not always human-like. Moreover, this discrepancy between
English and Japanese is further explored from the perspective of (non-)uniform
information density. Overall, our results suggest that a cross-lingual
evaluation will be necessary to construct human-like computational models.

##  1 Introduction

It is well known that the probability of a word in context (i.e., surprisal)
impacts its processing difficulty in incremental human language comprehension
Hale (2001); Demberg and Keller (2008); Levy (2008); Smith and Levy (2013).
Building on this basis, researchers have compared a variety of language models
(LMs) in terms of how well their surprisal correlates with human reading
behavior Roark et al. (2009); Frank and Bod (2011); Fossum and Levy (2012);
Hale et al. (2018); Goodkind and Bicknell (2018); Aurnhammer and Frank (2019);
Merkx and Frank (2020); Wilcox et al. (2020). Such investigations could
provide insights into the development of a general computational model of
human language processing. For example, recent studies reported that LMs with
better performance for next-word prediction could also better predict the
human reading behavior (i.e. more human-like) Fossum and Levy (2012); Goodkind
and Bicknell (2018); Wilcox et al. (2020).

Figure 1: Gaze duration from human subjects and surprisal from language models
for the Japanese sentence “Yononakaniwa samazamana hitoga irutoiu kotoga yoku
wakatta.” (I understood well that there are all kinds of people in the world.)

In this paper, we re-examine whether the recent findings on human-like
computational models can be generalized across languages. Despite the
community’s ongoing search for a language-independent model Bender (2011),
existing studies have focused almost exclusively on the English language.
Having said that, broad-coverage cross-linguistic evaluation of the existing
reports is prohibitively difficult. In fact, data on human reading behavior
(e.g., eye movement) is available only in limited languages. As an initial
foray, this study focuses on the Japanese language as a representative of
languages that have typologically different characteristics from the English
language. If the observation is different between English and Japanese, the
current findings on English data might lack a universality across languages.

We specifically revisit the recent report—the lower perplexity a LM has, the
more human-like the LM is—in the English and Japanese languages Fossum and
Levy (2012); Goodkind and Bicknell (2018); Wilcox et al. (2020). In addition
to the importance of cross-linguistic evaluation, the report itself is worth
investigating. Recent studies in the machine learning field have reported that
more parameters, training data, and computation cost can result in better PPL
Kaplan et al. (2020); Brown et al. (2020). Our investigation has implications
for whether a human-like model might exist beyond such improvements.

More concretely, over three dozens of LMs were trained for each language, with
variants in their architecture, training data size, and the number of
parameter updates. Then, the surprisals computed by each LM were compared to
human eye movement data (Figure 1). The analysis of the relationship between
PPL and the psychometric predictive power revealed substantively different
trends between the Japanese and English LMs. In Japanese, a lower PPL of a LM
does not indicate better performance for modeling reading behavior. By
contrast, in English, there was a clear relationship between the two metrics
as reported in the prior studies.

This opens a remaining and important question: why are English and Japanese
different in this aspect? We discuss the differing results between English and
Japanese from the perspective of the uniform information density hypothesis
Genzel and Charniak (2002); Levy (2005); Jaeger and Levy (2007). We find that
the processing difficulty (i.e., gaze duration) of segments is less uniformly
distributed within a Japanese sentence. Given this, the discrepancy of the
results between English and Japanese might stem from a mismatch between the
information uniformity of the target language and the LM’s training objective.
We demonstrate that tuning Japanese LMs to this training objective collapses
the human-like nonuniformity of the processing difficulty observed in Japanese
subjects. Our code is made publicly
available.111https://github.com/kuribayashi4/surprisal_reading_time_en_ja

##  2 Related work

###  2.1 Human sentence processing and LMs

What factor determines the incremental difficulty of human language
processing? At present, surprisal theory Hale (2001); Levy (2008) has been
widely adopted in the field of computational psycholinguistics. This theory
suggests that the processing difficulty of a segment is determined by how
predictable the segment is in its preceding context
($-\log{p(\mathrm{segment}|\mathrm{preceding\ context})}$).

Existing studies have compared various computational models by checking the
effectiveness of their surprisals in modeling human reading behavior Hale
(2001); Roark et al. (2009); Frank and Bod (2011); Fossum and Levy (2012);
Hale et al. (2018); Goodkind and Bicknell (2018); Merkx and Frank (2020);
Wilcox et al. (2020). Data such as eye movement Kennedy et al. (2003) and
brain activity Frank et al. (2015); Brennan et al. (2016) are used as measures
of human reading behavior. For example, using eye movement data, Frank and Bod
(2011) compared the surprisals from phrase-structure grammars (PSGs) with
those from a non-hierarchical, sequential model, tentatively concluding that
human sentence processing was insensitive to hierarchical structures since
non-hierarchical models displayed better psychological predictive power than
PSGs. Recently, researchers reported that surprisals from LMs with low PPL
correlate well with human reading behaviors Fossum and Levy (2012); Goodkind
and Bicknell (2018); Aurnhammer and Frank (2019); Wilcox et al. (2020).

The work most closely related to this study is Wilcox et al. (2020). They
examined the relationship between PPL, psychometric predictive power, and
syntactic knowledge in LMs using a variety of models, including modern neural
LMs Radrof et al. (2018). They found a tight relationship between PPL and
psychometric predictive power in the English corpora. This study investigates
whether this relationship can be generalized across languages.

###  2.2 Reading behavior in Japanese

In comparison to English speakers, Japanese speakers display different
patterns in sentence processing. For example, an anti-locality effect (the
more modifiers a word has in its preceding context, the easier the word is to
process) has typically been observed in head-final languages, including
Japanese Konieczny (2000). Such differences between the languages are assumed
to be more or less due to their different sentence structures. Recently, eye
movement data for naturally occurring Japanese texts have recently become
available Asahara et al. (2016) and was extensively annotated with various
linguistic properties Asahara and Kato (2017); Asahara (2017, 2018).

##  3 Methods

This section describes the settings of LMs, eye movement data, and evaluation
metrics.

###  3.1 Language models

A variety of sentence-level, left-to-right sequential LMs was used.

#### Training data of English LMs:

We used the WikiText-103 dataset to train the English LMs. Based on the
reports that subword-level English LMs exhibits superior psychometric
predictive power Wilcox et al. (2020), input texts were divided into subwords
by a byte-pair encoding (BPE) Sennrich et al. (2016).222Implemented in
SentencePiece Kudo and Richardson (2018). We set character coverage to
0.9995，and vocabulary size to 32,000 in English. In Japanese, the vocabulary
size is 100,000, reflecting its rich morphemes. The training data consist of
approximately 4M sentences (114M subwords units).

#### Training data of Japanese LMs:

We used news articles and the Japanese part of Wikipedia to train the Japanese
LMs. Input texts were first segmented into morphemes by MeCab Kudo (2006) with
unidic dictionary, and then further divided into subwords by
BPE.22footnotemark: 2 The training data consist of approximately 5M sentences
(146M subwords units).

#### Architectures:

The following four variants of LMs were used: Transformer-large (Trans-lg)
Vaswani et al. (2017), Transformer-small (Trans-sm), LSTM (LSTM) Hochreiter
and Schmidhuber (1997), and N-gram LMs (N-gram).333The neural LMs were trained
with the fairseq toolkit Ott et al. (2019). N-gram LMs were trained using
KenLM https://github.com/kpu/kenlm. The parameter size was almost the same for
Trans-sm and LSTM. With respect to the N-gram models, 3-gram, 4-gram, and
5-gram LMs were used. Appendix A shows the hyperparameters of the neural LMs.

#### Training data size:

For each neural LM architecture (Trans-lg, Trans-sm, and LSTM), three variants
were trained using different training data sizes: lg (full training data), md
(1/10 training data), and sm (1/100 training data). The N-gram LMs were
trained on lg datasets.

#### Number of updates:

The parameters of each neural LM were saved at four different points during
training: 100, 1K, 10K, and 100K parameter updates.

To summarize, 39 LM training settings were attained for each language (3
architectures $\times$ 3 data size $\times$ 4 parameter updates $=$ 36 neural
LMs, plus 3 N-gram LMs). In addition, our experiments use three LMs trained
using different random seeds for each neural LM training configure; hence, 111
LMs (36 neural LMs $\times$ 3 seeds, plus 3 N-gram LMs) were tested for each
language.

Corpus | #articles | #sents. | #segments |  #data points (used) |  #subjects per article |  Avg. GD per segment |  Avg. #subwords per segment
---|---|---|---|---|---|---|---
Dundee Corpus | 20 | 2,478 | 51,501 | 214,955 | 10 | 227.1 | 1.3
BCCWJ-EyeTrack | 20 | 218 | 1,643 | 6,009 | 12 | 361.6 | 3.4
Table 1: Statistics of the corpora used for evaluating the psychometric
predictive power of LMs. “#articles” and “#sents.” are the number of articles
and sentences in each corpus. “#segments” denotes the number of segments
annotated with human reading time in each corpus. “#data points” is the number
of reading time annotations used in our experiments. Each segment has the
reading time annotations from multiple subjects (#subjects per article). “Avg.
GD per segment” is the averaged gaze duration per segment. “Avg. #subwords per
segment” denotes the averaged number of subwords consisting of each segment.

###  3.2 Eye movement data

#### English:

The Dundee Corpus Kennedy et al. (2003), which contains gaze duration
annotation for each word, was used. Following Smith and Levy (2013), the
first-pass gaze duration was analyzed. Then, following Goodkind and Bicknell
(2018), the data points that met any of the following criteria were excluded:

  * •

data points with zero gaze duration or that beyond three standard deviations

  * •

segments with punctuation or numeric characters

  * •

segments whose next segment has punctuation or numeric characters

  * •

first or last segment in a line

In total, the analysis included 214,955 data points in the corpus.

#### Japanese:

The BCCWJ-EyeTrack Asahara et al. (2016), which contains gaze duration
annotation for each phrasal unit, was used. Note that the phrasal unit (i.e.,
bunsetsu) consists of at least one content morpheme and its postpositional
function morphemes. Henceforth, an English word and a Japanese phrasal unit
are referred to as a “segment.” The same exclusion criteria as the Dundee
Corpus was applied to the BCCWJ-EyeTrack data.444Strictly speaking, the
exclusion criteria was slightly different between Japanese and English data.
In the Japanese data, we included the segments whose next segment had
punctuation or a numeric character, as there is no spillover effect in
Japanese (see Section 3.3) In total, the analysis included 6,009 data points
in the corpus. Note that the BCCWJ-EyeTrack data was deliberately designed to
address language-specific issues in Japanese such as the lack of segmentation
spaces in Japanese texts Asahara et al. (2016).

#### Statistics:

Table 1 shows the statistics of the Dundee Corpus and BCCWJ-EyeTrack data. The
BCCWJ-EyeTrack has more than 10 times a smaller number of data points than the
Dundee Corpus. Notably, the segment annotated with eye movement information
differs between English and Japanese. On average, a Japanese segment consists
of 3.4 subwords, while an English segment consists of 1.3 subwords. Smith and
Levy (2013) theoretically proved that the more fragments a word is divided
into when computing its surprisal, the better the calculated surprisal
approximates the human cognitive effort if the human language processing is
highly incremental. Thus, we tentatively consider that this difference did not
make a negative impact on the results using the Japanese data.

###  3.3 Evaluation metrics

#### Perplexity (PPL):

PPL, the inverse geometric mean of next-word probabilities $p(w_{i}|w_{<i})$
in a text that consists of $N$ signals ($w_{1},w_{2},\cdots,w_{N}$), is a
typical evaluation metric for unidirectional LMs (Eq. 1):

PPL $\displaystyle=\prod_{i=0}^{N}p(w_{i}|w_{<i})^{-\frac{1}{N}}\>\>.$ (1)

Figure 2: Relationship between PPL (X-axis) and psychometric predictive power,
i.e., $\Delta$LogLik (Y-axis) in the English and Japanese languages. Each
point corresponds to each LM. A low score on the X-axis indicates the high
linguistic accuracy of the model. The PPL was calculated on the eye movement
data, and the LMs with PPL more than 106 were excluded from the figure. A high
score on the Y-axis indicates that the model has a high psychometric
predictive power. Note that the X-axis is on a log scale.

Low PPL indicates that the model can accurately predict the upcoming signal
based on its preceding context. The training objective of LMs works to
minimize the PPL computed by the model. In the experiments, the PPL of a LM is
evaluated with the texts in the eye movement data, which do not overlap with
the training data. A model with low PPL is also called a linguistically
accurate model Frank and Bod (2011).

#### Psychometric predictive power:

The surprisal measure, a negative logarithmic probability of a segment in
context ($-\log{p({\rm segment}|{\rm preceding\ context})}$), is a widely used
information-theoretic complexity metric. Intuitively, a model is considered to
have high psychometric predictive power (i.e., psychological accuracy) if the
surprisals of segments computed by the model have trends similar to the human
subject’s cognitive load (e.g., measured by gaze duration). Following the
existing studies Goodkind and Bicknell (2018); Merkx and Frank (2020); Wilcox
et al. (2020), the psychometric predictive power of a model was measured by
comparing surprisal from the model and gaze duration from human subjects.

While LMs process a text subword-by-subword, gaze duration is annotated in a
larger segment. Following the study using subwords Wilcox et al. (2020), the
surprisal of each segment was calculated using the joint probability of its
constituent subwords. Formally, given a text consisting of $N$ subwords
$w_{1:N}=(w_{1},w_{2},\cdots,w_{N})$, surprisal $I(\cdot)$ of a segment
$s_{k}=(w_{l},w_{l+1},\cdots,w_{m})$, where $1\leq l\leq m\leq N$, was
calculated as follows:

$\displaystyle\begin{split}I(s_{k})&=-\log p(w_{l},\cdots,w_{m}|w_{<l})\\\
&=-\sum_{k=l}^{m}\log p(w_{k}|w_{1},\cdots,w_{k-1})\>\>.\end{split}$ (2)

The effect of surprisals for modeling human reading behavior was calculated
using a linear mixed-effects regression Bates et al. (2015). Specifically, the
gaze duration (GD) was modeled using the following formula:

$\displaystyle\begin{split}\texttt{GD}&\sim\texttt{surprisal}+\texttt{surprisal\\_prev\\_1}\\\
&+\texttt{surprisal\\_prev\\_2}+\texttt{freq}*\texttt{length}\\\
&+\texttt{freq\\_prev\\_1}*\texttt{length\\_prev\\_1}\\\
&+\texttt{screenN}+\texttt{lineN}+\texttt{segmentN}\\\
&+\texttt{(1|article)}+\texttt{(1|subj)}\>\>.\end{split}$ (3)

The regression model includes baseline factors (e.g., frequency of a segment)
that are of no interest in the comparison of LMs. A collection of factors used
in the existing studies Asahara et al. (2016); Wilcox et al. (2020) were
initially examined and the factors that were not significant ($p>0.05$) for
gaze duration modeling both in the Dundee Corpus and BCCWJ-EyeTrack were
excluded. The frequency of a segment (freq) was calculated using the entire
training data for LMs. Appendix B shows the details of each factor in Eq. 3.

In English experiments, surprisals of preceding words (surprisal_prev_1 and
surprisal_prev_2) were included in order to handle the spillover effect (the
processing cost of a certain segment is affected by its preceding segments)
Rayner and Well (1996); Smith and Levy (2013). In Japanese experiments, the
surprisals of preceding words were not included because our preliminary
experiment showed that these factors were not significantly effective for
modeling gaze duration in the BCCWJ-EyeTrack.555The reason is probably that a
Japanese phrasal unit (i.e., bunsetsu) could be a larger unit than an English
word. All the regression models used in our experiments were converged.

Figure 3: Separate effect of model architecture, training data size, and the
number of parameter updates for LMs’ psychometric predictive power in each
language. Each point corresponds to each LM. The box shows the quartiles of
the data. The whiskers show 1.5 times interquartile range.

To isolate the effect of surprisal for gaze duration modeling, a baseline
regression model was trained without surprisal information (excluding the
surprisal, surprisal_prev_1, and surprisal_prev_2 terms from Eq. 3). Following
Wilcox et al. (2020), the mean by-segment difference of log-likelihood between
the model using surprisal values (Eq. 3) and the baseline model was
calculated. Henceforth, this metric is called $\Delta$LogLik. When surprisal
from a LM is not effective for gaze duration modeling, the $\Delta$LogLik
score becomes zero. A high $\Delta$LogLik means that the surprisal values
obtained by the LM are effective for modeling gaze duration (i.e., the LM has
a high psychometric predictive power).

##  4 Experiments

The relationship between PPL and psychometric predictive power is
investigated. Furthermore, the relationship is analyzed with respect to the
training configures of LMs (e.g., the number of parameter updates). Then, we
discuss the results from the perspective of the uniformity of information
density.

###  4.1 Psychometric predictive power and PPL

Figure 2 shows the relationship between PPL and psychometric predictive power
(i.e., $\Delta$LogLik) of LMs in each of the languages. Each point corresponds
to each LM, and a score on the Y-axis indicates the psychometric predictive
power of a LM (higher is better). The X-axis is PPL on a log scale (lower is
better).

#### Dundee Corpus:

First, the results of the data from the Dundee Corpus show a clear
relationship between PPL and psychometric predictive power; namely, lower PPL
corresponds to more psychometric predictive power, as reported by prior
studies Goodkind and Bicknell (2018); Wilcox et al. (2020). Spearman’s rank
correlation coefficient between the two metrics was $-0.87$.

#### BCCWJ-EyeTrack:

By contrast, in BCCWJ-EyeTrack, there was no clear, consistent trend between
the PPL and psychometric predictive power. While LMs with PPL over 400 show
the correlation between PPL and psychometric predictive power ($-0.68$ with
Spearman’s $\rho$), there is a positive correlation ($0.53$ with Spearman’s
$\rho$) for LMs with PPL below 400. The positive correlation means that the
more accurately the LMs can predict the upcoming word, the worse the
psychometric predictive power of the LMs is. These results demonstrate the
non-universality of the recent report across languages; lower perplexity is
not always human-like. The lstm LM trained using the md dataset with 1K
updates achieved the best psychometric predictive power. Notably, surprisal
was effective for gaze duration modeling in all the Japanese LMs.
$\Delta$logLik scores were significantly higher than zero with the chi-square
test ($p<$0.05).

###  4.2 Model architectures, data sizes, number of parameter updates

Which factor (e.g., model architecture, training data size, and the number of
parameter updates) characterizes the psychometric predictive power of LMs? Is
the collection of effective factors consistent between the two languages? This
study takes a more in-depth look at the separate effects of (i) model
architecture, (ii) training data size, and (iii) the number of parameter
updates for the psychometric predictive power.

Figure 3 summarizes the effect of each factor, where the Y-axis denotes the
psychometric predictive power. The most noticeable trend is that Japanese LMs
with a relatively fewer number of parameter updates (1K) have better
psychometric predictive power than the other Japanese LMs (bottom right part
of Figure 3), while this trend does not exist in the English LMs (top right
part). This implies that the training objective of the LMs, maximizing
$\frac{1}{N}\sum_{i=1}^{N}\log P(w_{i}|w_{<i})$, had a negative impact on the
psychometric predictive power of LMs, at least in Japanese. We discuss this
point in Section 4.3.

To quantitatively test the differences in Figure 3, a linear regression model
was trained to estimate psychometric predictive power with the factors of the
model architecture, the training data size, and the parameter update number in
each language. The training data size and the parameter update number are
represented as logarithmically transformed numerical factors. The following
trends were found: (i) ; (ii) the training data size positively affects the
performance in English alone; and (iii) the number of parameter updates
positively affects the performance only in English. There was no factor that
boosted the psychometric predictive power of LMs in both English and Japanese
languages.

###  4.3 Discussion: uniform information density

The key question is: why do Japanese and English show different trends between
PPL and psychometric predictive power? One possible interpretation connecting
our results to the uniform information density is discussed in this section.

In computational psycholinguistics, it is commonly assumed that language is
designed to enable efficient communication. This principle has been typically
investigated under the uniform information density (UID) hypothesis Genzel and
Charniak (2002); Levy (2005); Jaeger and Levy (2007). This hypothesis suggests
that speakers seek to keep the amount of information constant across the
signals (e.g., segments).

Assuming this hypothesis holds for all languages, the reasonable expectation
would be for human subjects to show a near-uniform gaze duration across
segments regardless of their native language. However, this study found that
the coefficient of variation666Coefficient of variation is
$\frac{\sigma}{\mu}$, where $\sigma$ and $\mu$ are the standard deviation and
the mean of the first-pass gaze durations in the eye movement data. in gaze
duration over the whole corpus was around 1.7 times higher in Japanese
compared to English (0.75 vs. 0.44). Specifically, in Japanese, the gaze
duration tended to speed up towards the end of sentences, whereas the duration
was near-uniform in English (Figure 4).777At least in our experimental setup,
token position within the sentence was not significantly effective for gaze
duration modeling in English sentences, whereas it was significant in Japanese
sentences. We checked the coefficient of the factor of position in sentence
segmentN using the linear regression model of GD $\sim$ sengmentN. These
observations imply that the Japanese language might have a less uniform
information density than English. This phenomenon was also investigated
through the lens of word order, where SOV languages such as Japanese are
reported to show less uniformity of information density Maurits et al. (2010).

Based on this observation, the discrepancy between English and Japanese low-
PPL LMs’ psycholinguistic predictive power could stem from a mismatch between
the LM’s training objective and the information uniformity of the target
language. The objective function, $\frac{1}{N}\sum_{i=1}^{N}\log
P(w_{i}|w_{<i})$, defines that the “ideal” is to maximize all next word
probabilities to 1.0 (a uniform goal).888PPL,
$\prod_{i=1}^{N}P(w_{i}|w_{<i})^{-\frac{1}{N}}$, is minimized when the LM
objective are maximized. That is, LMs are, in theory, trained to approach a
model satisfying the UID assumption Bloem (2016), where all surprisals from
the LM are equally, sufficiently small across the segments. Therefore, the
objective function might lead to a worse approximation of human-like surprisal
in languages that are further from the UID assumption, such as Japanese, while
it might be more compatible with English, which has a more uniform processing
difficulty across segments. This explanation would be consistent with the
observation that more tuning to the LM training objective (i.e., a lower PPL)
had a negative impact on the psycholinguistic performance of the Japanese LMs
(Section 4.2). Note the tendency of LMs to assign unreasonably high
probabilities to segments has also attracted attention from the viewpoint of
memorization capability of LMs Carlini et al. (2020). In addition, the
connection of the UID hypothesis to the modern NLP techniques has been
recently explored Meister et al. (2020); Wei et al. (2021). We further
investigate our hypothesis in Section 5.

Figure 4: Uniformity of gaze duration with respect to segment position in a
sentence. This plot is computed by the generalized additive model of GD $\sim$
segmentN. Here, segmentN is denoted as the position of a segment in a
sentence. Figure 5: Relationship between the LM’s psychometric predictive
power and the effect of the syntactic category on the surprisal computed by
each LM (left part), and that between PPL and the effect of the syntactic
category (right part). Each point corresponds to each LM. The PPL was
calculated on the eye movement data, and the LMs with PPL more than 106 was
excluded from the right part of the figure.

##  5 Probing nonuniform information density of Japanese LMs

This study hypothesized that tuning to the LM objective (i.e., uniform goal)
obscures the nonuniform trend observed in the reading behavior of Japanese
subjects. We investigated whether the nonuniformity of the processing
difficulty observed in human reading time is mirrored by LM surprisals.

#### Settings:

In a preliminary experiment, we observed that the _syntactic category_
(similar to part-of-speech) was the most dominant linguistic factor for
explaining the difference in human gaze duration in Japanese sentences (see
Appendix D). Based on this observation, we analyze the nonuniformity of
surprisals in Japanese LMs with respect to the syntactic categories.

The segments in BCCWJ-EyeTrack were classified into one of the following
syntactic categories: (a) nominal (nouns), (b) verbal (verbs), (c) modifier
(adjectives and adverbs), and (d) other entries, as follows:

Kanojo-ga akai kaban-o kat-ta She-NOM red bag-ACC buy-PAST nominal modifier
nominal verbal

As Asahara and Kato (2017) reported, verbal and modifier segments have a
shorter gaze duration than the other segments in Japanese sentences. An
analysis was conducted on how strongly the Japanese LM’s surprisals on
segments are influenced by their syntactic category. This influence can be
evaluated by examining how effectively syntactic category factors can model LM
surprisals.

In this experiment, surprisal was regarded as “simulated gaze duration” from
an “LM subject,” and the importance of syntactic category information for
modeling the simulated gaze duration (simulated_GD) was evaluated. To inspect
the effect of the syntactic category information for modeling the simulated
gaze duration, the following regression model999sentN and tokenN denote the
sentence position and the segment position in a sentence (see Appendix B).
Note that the tokenN and syntactic category exhibit low correlation (0.02 with
Pearson’s $r$). was used, including a factor defining which syntactic category
the segment falls into (syn_category):

$\displaystyle\begin{split}\texttt{simulated\\_GD}&\sim\texttt{syn\\_category}+\texttt{sentN}\\\
&+\texttt{tokenN}+\texttt{freq}*\texttt{length}\>\>.\end{split}$ (4)

From this regression model, a log-likelihood score for the simulated gaze
duration was obtained. To evaluate the separate effect of syn_category,
$\Delta$LogLik between Eq. 4 and a baseline model was calculated. The baseline
model was
$\texttt{simulated\\_GD}\sim\texttt{sentN}+\texttt{tokenN}+\texttt{freq}*\texttt{length}$.
The $\Delta$LogLik is denoted as “Effect of syntactic category.” A lower score
means that the LM lacked the property of varying processing difficulty with
respect to the syntactic category.

#### Results:

The results are shown in Figure 5. First, the higher psychometric predictive
power the LMs exhibit, the greater the effect of syntactic category on
surprisals (left part in Figure 5). This means that, depending on the
syntactic category of the segment they processed, LMs with high psychometric
predictive power computed surprisals with a more nonuniform trend. The right
part of Figure 5 shows that, as PPL decreases below a certain value
($\mathrm{PPL}\sim 400$), the Japanese LMs compute surprisals that obscure the
nonuniform trends with respect to the syntactic category of segments.101010The
correlation between PPL and the effect of syntactic category in the LMs with
PPL less than 400 was 0.45 and 0.34 with Pearson’s $r$ and Spearman’s $\rho$,
respectively. This trend supports our hypothesis that tuning to LM objectives
obscures the human-like nonuniformity of the processing difficulty. Even
though LMs that are not fully tuned to the LM objective ($\mathrm{PPL}\sim
400$) acquire human-like trends with respect to syntactic category, these
biases tend to be lost by further lowering their PPL.

Notably, we also observed that not all the types of linguistic nonuniformity
were obscured in surprisals computed by the LMs with low PPL. For example,
Appendix E shows that LMs with lower PPL compute surprisals that better
correlates with a particular syntactic factor although that factor is a less
dominant trend in human reading behavior than the syntactic category (Appendix
D).

##  6 Limitations and future works

To test the universality of the recent findings in computational
psycholinguistics across languages, the initial focus is on English and
Japanese as a pair of languages with different linguistic properties. Although
the discrepancy of the results in the two languages is discussed from the
viewpoint of the UID hypothesis, the two languages are also different in
various ways, such as writing systems, agglutinative property, case marking,
sentence structure, and pro-drop nature. To identify the difference that
relates to the human-like behaviors of LMs, experiments that include
additional languages should be conducted in the future.

In addition, the corpus size of the BCCWJ-EyeTrack data is smaller than the
Dundee Corpus. While the reading time data in the BCCWJ-EyeTrack was collected
from various human subjects, the number of the independent segments was
limited (1,643 segments, 218 sentences). Thus, whether the trends reported in
this study generalize to more diverse Japanese texts should be explored in
future work. It is hoped that this study motivates the creation of a large-
scale corpus of human reading behaviors in diverse languages.

##  7 Conclusion

This study has investigated whether the recent reports on the psychometric
predictive power of LMs can be generalized across languages. Our initial
investigation has re-examined the recent report—the lower PPL a LM has, the
more human-like the LM is—using Japanese eye movement data. Our experiments
have demonstrated a surprising lack of universality of this report; lower
perplexity is not always human-like. This discrepancy of the results between
the languages reinforces the need for the cross-lingual evaluation of the
psychometric predictive power of LMs. The discussion considers potential
factors that make the observation different across languages from the
viewpoint of the uniform information density hypothesis. We believe that this
is an important first step for seeking a language-agnostic model of human
sentence processing. Hopefully, this study encourages researchers to further
investigate the universality of human language processing across languages.

## Acknowledgements

We would like to thank the members at the Tohoku NLP Lab for their valuable
advice, particularly Ana Brassard for proofreading. This work was supported by
Grant-in-Aid for JSPS Fellows Grant Number JP20J22697, JSPS KAKENHI Grant
Number 19H04990, and JST CREST Grant Number JPMJCR20D2. This work was also
supported by the National Institute for Japanese Language and Linguistics
(NINJAL) Collaborative Research Project “Computational Psycholinguistics of
Language Processing with Large Corpora.”

## Ethical considerations

Language models with low perplexity are typically trained with a high
computational cost. Our work demonstrates that further up-scaling the models
might not be a reasonable direction of searching for human-like language
models. This could potentially contribute to reducing energy and carbon costs,
which are needed to train large-scale language models.

## References

  * Asahara (2017) Masayuki Asahara. 2017.  Between Reading Time and Information Structure.  In _Proceedings of PACLIC_ , pages 15–24.
  * Asahara (2018) Masayuki Asahara. 2018.  Between Reading Time and Clause Boundaries in Japanese - Wrap-up Effect in a Head-final Language.  In _Proceedings of PACLIC_ , pages 19–27.
  * Asahara and Kato (2017) Masayuki Asahara and Sachi Kato. 2017.  Between Reading Time and Syntactic / Semantic Categories.  In _Proceedings of IJCNLP_ , pages 404–412.
  * Asahara et al. (2016) Masayuki Asahara, Hajime Ono, and Edson T Miyamoto. 2016.  Reading-Time Annotations for “Balanced Corpus of Contemporary Written Japanese”.  In _Proceedings of COLING_ , pages 684–694.
  * Aurnhammer and Frank (2019) C Aurnhammer and S L Frank. 2019.  Comparing Gated and Simple Recurrent Neural Network Architectures as Models of Human Sentence Processing.  In _Proceedings of CogSci_ , pages 112–118.
  * Bates et al. (2015) Douglas Bates, Martin Mächler, Ben Bolker, and Steve Walker. 2015.  Fitting linear mixed-effects models using lme4.  _Journal of Statistical Software, Articles_ , 67(1):1–48.
  * Bender (2011) Emily M Bender. 2011.  On Achieving and Evaluating Language-Independence in NLP.  _Linguistic Issues in Language Technology_ , 6(3):1–26.
  * Bloem (2016) Jelke Bloem. 2016.  Testing the processing hypothesis of word order variation using a probabilistic language model.  In _Proceedings of the Workshop on Computational Linguistics for Linguistic Complexity (CL4LC)_ , pages 174–185, Osaka, Japan. The COLING 2016 Organizing Committee.
  * Brennan et al. (2016) Jonathan R Brennan, Edward P Stabler, Sarah E Van Wagenen, Wen-Ming Luh, and John T Hale. 2016.  Abstract linguistic structure correlates with temporal activity during naturalistic comprehension.  _Brain and language_ , 157:81–94.
  * Brown et al. (2020) Tom B Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel M Ziegler, Jeffrey Wu, Clemens Winter, Christopher Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei. 2020.  Language Models are Few-Shot Learners.  In _Proceedings of NeurIPS_.
  * Carlini et al. (2020) Nicholas Carlini, Florian Tramèr, Eric Wallace, Matthew Jagielski, Ariel Herbert-Voss, Katherine Lee, Adam Roberts, Tom B. Brown, Dawn Song, Úlfar Erlingsson, Alina Oprea, and Colin Raffel. 2020.  Extracting training data from large language models.  _CoRR_ , abs/2012.07805.
  * Demberg and Keller (2008) Vera Demberg and Frank Keller. 2008.  Data from eye-tracking corpora as evidence for theories of syntactic processing complexity.  _Journal of Cognition_ , 109(2):193–210.
  * Fossum and Levy (2012) Victoria Fossum and Roger Levy. 2012.  Sequential vs. Hierarchical Syntactic Models of Human Incremental Sentence Processing.  In _Proceedings of CMCL_ , pages 61–69, Montréal, Canada.
  * Frank and Bod (2011) Stefan L Frank and Rens Bod. 2011.  Insensitivity of the Human Sentence-Processing System to Hierarchical Structure.  _Psychological science_ , 22(6):829–834.
  * Frank et al. (2015) Stefan L. Frank, Leun J. Otten, Giulia Galli, and Gabriella Vigliocco. 2015.  The erp response to the amount of information conveyed by words in sentences.  _Brain and Language_ , 140:1–11.
  * Genzel and Charniak (2002) Dmitriy Genzel and Eugene Charniak. 2002.  Entropy rate constancy in text.  In _Proceedings of the 40th Annual Meeting of the Association for Computational Linguistics_ , pages 199–206, Philadelphia, Pennsylvania, USA. Association for Computational Linguistics.
  * Goodkind and Bicknell (2018) Adam Goodkind and Klinton Bicknell. 2018.  Predictive power of word surprisal for reading times is a linear function of language model quality.  In _Proceedings of CMCL2018_ , pages 10–18.
  * Hale (2001) John Hale. 2001.  A Probabilistic Earley Parser as a Psycholinguistic Model.  In _Proceedings of NAACL_ , pages 159–166.
  * Hale et al. (2018) John Hale, Chris Dyer, Adhiguna Kuncoro, and Jonathan R. Brennan. 2018.  Finding Syntax in Human Encephalography with Beam Search.  In _Proceedings of ACL_ , pages 2727–2736.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. 1997.  Long Short-Term Memory.  _Journal of Neural Computation_ , 9(8):1735–1780.
  * Jaeger and Levy (2007) T Jaeger and Roger Levy. 2007.  Speakers optimize information density through syntactic reduction.  In _Advances in Neural Information Processing Systems_ , volume 19, pages 849–856. MIT Press.
  * Kaplan et al. (2020) Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei. 2020.  Scaling Laws for Neural Language Models.  _arXiv preprint arXiv:2001.08361_.
  * Kennedy et al. (2003) Alan Kennedy, Robin Hill, and Joël Pynte. 2003.  The dundee corpus.  In _Proceedings of the 12th European conference on eye movement_.
  * Konieczny (2000) Lars Konieczny. 2000.  Locality and Parsing Complexity.  _Journal of Psycholinguistic Research_ , 29(6):627–645.
  * Kudo (2006) Taku Kudo. 2006.  MeCab: Yet Another Part-of-speech and Morphological Analyzer.  _http://mecab. sourceforge. jp_.
  * Kudo and Richardson (2018) Taku Kudo and John Richardson. 2018.  SentencePiece: A simple and language independent subword tokenizer and detokenizer for Neural Text Processing.  In _Proceedings of EMNLP_ , pages 66–71.
  * Levy (2005) Roger Levy. 2005.  _Probabilistic models of word order and syntactic discontinuity_.  stanford university.
  * Levy (2008) Roger Levy. 2008.  Expectation-based syntactic comprehension.  _Journal of Cognition_ , 106(3):1126–1177.
  * Maurits et al. (2010) Luke Maurits, Dan Navarro, and Amy Perfors. 2010.  Why are some word orders more common than others? a uniform information density account.  In _Advances in Neural Information Processing Systems_ , volume 23. Curran Associates, Inc.
  * Meister et al. (2020) Clara Meister, Ryan Cotterell, and Tim Vieira. 2020.  If beam search is the answer, what was the question? In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 2173–2185, Online. Association for Computational Linguistics.
  * Merkx and Frank (2020) Danny Merkx and Stefan L. Frank. 2020.  Human Sentence Processing: Recurrence or Attention? In _Procceding of CMCL 2021_.
  * Ott et al. (2019) Myle Ott, Sergey Edunov, Alexei Baevski, Angela Fan, Sam Gross, Nathan Ng, David Grangier, and Michael Auli. 2019.  fairseq: A fast, extensible toolkit for sequence modeling.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics (Demonstrations)_ , pages 48–53, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Radrof et al. (2018) Alec Radrof, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2018.  Language Models are Unsupervised Multitask Learners.  _Technical report, OpenAI_.
  * Rayner and Well (1996) Keith Rayner and Arnold D Well. 1996.  Effects of contextual constraint on eye movements in reading: A further examination.  _Psychonomic Bulletin & Review_, 3(4):504–509.
  * Roark et al. (2009) Brian Roark, Asaf Bachrach, Carlos Cardenas, and Christophe Pallier. 2009.  Deriving lexical and syntactic expectation-based measures for psycholinguistic modeling via incremental top-down parsing.  In _Proceedings of EMNLP_ , pages 324–333, Singapore.
  * Sennrich et al. (2016) Rico Sennrich, Barry Haddow, and Alexandra Birch. 2016.  Neural Machine Translation of Rare Words with Subword Units.  In _Proceedings of ACL_ , pages 1715–1725.
  * Smith and Levy (2013) Nathaniel J. Smith and Roger Levy. 2013.  The effect of word predictability on reading time is logarithmic.  _Journal of Cognition_ , 128(3):302–319.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention Is All You Need.  In _Proceedings of NIPS_ , pages 5998–6008.
  * Wei et al. (2021) Jason Wei, Clara Meister, and Ryan Cotterell. 2021.  A cognitive regularizer for language modeling.  _arXiv preprint arXiv:2105.07144_.
  * Wilcox et al. (2020) Ethan Gotlieb Wilcox, Jon Gauthier, Jennifer Hu, Peng Qian, and Roger Levy. 2020\.  On the Predictive Power of Neural Language Models for Human Real-Time Comprehension Behavior.  In _Proceedings of CogSci_ , pages 1707–1713.

Figure 6: Relationship between PPL (X-axis) and psychometric predictive power
(Y-axis). Each point corresponds to each LM. Low score on X-axis indicates the
high linguistic accuracy of the model. High score on Y-axis indicates that the
model has a high psychometric predictive power. Note that X-axis is on a log
scale. The shape, color, and size of each point is same as Figure 2. Figure 7:
Relationship between PPL (X-axis) and the effect of the anti-locality
(Y-axis). Each point corresponds to each LM. A low score on X-axis indicates
the high linguistic accuracy of the model. A high score on Y-axis indicates
that the surprisals computed by the corresponding model are highly biased
towards the anti-locality effect. Note that X-axis is on a log scale. The
shape, color, and size of each point is the same as Figure 2.

##  Appendix A Hyperparameters of LMs

Table 2 shows the hyperparameters of Trans-sm, Trans-lg, and LSTM,
respectively. Note that the number of parameter updates varies as described in
Section 3.

Fairseq model | architecture | transformer_lm_gpt2_small
---|---|---
adaptive softmax cut off | 50,000, 140,000
share-decoder-input-output-embed | True
embed_dim | 1,024
ffn_embed_dim | 4,096
layers | 24
heads | 16
dropout | 0.1
attention_dropout | 0.1
Optimizer | algorithm | AdamW
learning rates | 5e-4
betas | (0.9, 0.98)
weight decay | 0.01
clip norm | 0.0
Learning rate scheduler | type | inverse_sqrt
warmup updates | 4,000
warmup init lrarning rate | 1e-7
Training | batch size | 61,440 tokens
sample-break-mode | none
(a) Trans-lg.

Fairseq model | architecture | transformer_lm_gpt
---|---|---
adaptive softmax cut off | 50,000, 140,000
share-decoder-input-output-embed | True
embed_dim | 384
ffn_embed_dim | 2,048
layers | 8
heads | 6
dropout | 0.1
attention_dropout | 0.1
Optimizer | algorithm | AdamW
learning rates | 5e-4
betas | (0.9, 0.98)
weight decay | 0.01
clip norm | 0.0
Learning rate scheduler | type | inverse_sqrt
warmup updates | 4,000
warmup init lrarning rate | 1e-7
Training | batch size | 61,440 tokens
sample-break-mode | none
(b) Trans-sm.

Fairseq model | architecture | lstm_lm
---|---|---
adaptive softmax cut off | 50,000, 140,000
share-decoder-input-output-embed | True
embed_dim | 400
hiden_size | 1,024
layers | 2
dropout | 0.1
Optimizer | algorithm | AdamW
learning rates | 1e-3
betas | (0.9, 0.98)
weight decay | 0.01
clip norm | 0.0
Learning rate scheduler | type | inverse_sqrt
warmup updates | 4,000
warmup init lrarning rate | 1e-7
Training | batch size | 20,480 tokens
sample-break-mode | none
(c) LSTM.

Table 2: Hyperparameters for the LMs.

##  Appendix B Factors used in regression models

Descriptions for the factors used in our experiments are shown in Table 3. The
frequency of a segment (freq) was estimated using the full training data for
the LMs.

Factor name | Type |  Description
---|---|---
surprisal | num |  surprisal caluzulted by LMs
GD | num |  reading time (first pass time)
article | factor |  article ID
screenN | int |  screen display order
lineN | int |  the serial number of line the segment is displayed
segmentN | int |  the serial number of segment in a screen
sentN | int |  the serial number of sentence the segment belongs to
tokenN | int |  the position of segment in sentence
length | int |  number of characters
freq | num |  geometric mean of the frequencies of subword constituents in a segment
subj | factor |  participant ID
syn_category | factor |  syntactic category the segment falls into (nominal, verbal, modifier, or other)
sem_category | factor |  semantic category the segment falls into (relation, subject, action, product, or nature)
n_dependents | int |  number of dependents before the segment
Table 3: Factor names and their description. syntactic category | number of segments | Avg. gaze duration
---|---|---
nominal | 4,322 | 388.4
verbal | 1,090 | 291.0
modifier | 588 | 297.1
other | 9 | 239.3
Table 4: The statistics of the syntactic category labels in BCCWJ-EyeTrack.

##  Appendix C Results of modeling logarithmic gaze duration in BCCWJ-EyeTrack

Existing studies Asahara et al. (2016) performed experiments using the
logarithmic gaze duration because the logarithmic gaze duration more matches
the normal distribution than the raw gaze duration. Given this, we
additionally conducted experiments in Section 4, changing the target variable
from the raw gaze duration to its logarithmic gaze duration. The result with
this setting is shown in Figure 6. There was no substantial difference with
the results shown in Section 4.

##  Appendix D Preliminary experiments in Section 5

Which linguistic factor is helpful for explaining the difference in gaze
duration? We conducted experiments using linguistic annotation in the BCCWJ-
EyeTrack. Following the existing studies, we checked the separate effect of
syntactic category, semantic category Asahara and Kato (2017), and a
particular aspect of hierarchical syntactic structure (i.e., the anti-locality
effect) Asahara et al. (2016). Specifically, we used the factors,
syn_category, sem_category, and n_dependents, shown in Table 3. For each
factor, we inspect the separate effect of each factor for modeling gaze
duration. As Eq. 4, we first modeled the gaze duration using each factor
(factor_X):

$\displaystyle\begin{split}\texttt{GD}&\sim\texttt{factor\\_X}+\texttt{sentN}\\\
&+\texttt{segmentN}+\texttt{freq}*\texttt{length}\>\>.\end{split}$ (5)

Then, we calculated the $\Delta$LogLik between X and a baseline model. The
baseline model was
$\texttt{GD}\sim\texttt{sentN}+\texttt{segmentN}+\texttt{freq}*\texttt{length}$.

The $\Delta$LogLik for each collection of factors are shown in 5. We found
that syntactic category is the most influential factor for modeling gaze
duration, at least in this experiment.

linguistic property |  $\Delta$LogLik
---|---
syntactic category | 58.37
semantic category | 17.08
number of dependents | 13.84
Table 5: The separate effect of each linguistic annotation for modeling gaze
duration.

##  Appendix E Anti-locality effect in LMs

Similar to Section 5, we analyzed how strongly the surprisals from each
Japanese LM are biased towards a particular linguistic property. In this
section, we investigated the anti-locality effect in the surprisals from LMs.
The anti-locality is that the more dependents a segment has in its preceding
context, the cognitive effort of the head segment is reduced (i.e., modifiers
alleviate the processing cost of their head).

Analogous to the Section 5, we regarded surprisal as “simulated gaze duration”
from an “LM subject,” and evaluated the importance of the number of the
dependents in its preceding context (n_dependents) for modeling the simulated
gaze duration (simulated_GD). To inspect the effect of the n_dependents for
modeling the simulated gaze duration, we used the following regression model:

$\displaystyle\begin{split}\texttt{simulated\\_GD}&\sim\texttt{n\\_dependents}+\texttt{sentN}\\\
&+\texttt{tokenN}+\texttt{freq}*\texttt{length}\>\>.\end{split}$ (6)

From this regression model, we obtained a log-likelihood score for the
simulated gaze duration. To evaluate the separate effect of n_dependents, we
calculated the $\Delta$LogLik between Eq. 6 and a baseline model. The baseline
model was
$\texttt{simulated\\_GD}\sim\texttt{sentN}+\texttt{segmentN}+\texttt{freq}*\texttt{length}$.
The $\Delta$LogLik is denoted as “Effect of the anti-locality.”

The results are shown in Figure 7. There is a clear trend that the LMs with
lower PPL exhibit surprisals that are more consistent with the anti-locality
effect (Spearman’s $\rho=-0.77$ between PPL and the strength of the anti-
locality effect). This suggests that the surprisals from LMs with low PPL are
biased towards the hierarchical structure of sentences rather than the
syntactic category.


#  Infrared Properties of OGLE4 Mira variables in our Galaxy

Kyung-Won Suh<EMAIL_ADDRESS>
(Received September 16, 2022; accepted November 22, 2022)

##  1 Introduction

Most asymptotic giant branch (AGB) stars are believed to be long-period
variables (LPVs) with outer dust envelopes (e.g., Suh 2021). LPVs are
classified into small-amplitude red giants (SARGs), semiregular variables
(SRVs), and Mira variables (e.g., Soszyński et al. 2013). Though most of SRVs
and a majority of SARGs could also be in the AGB phase, all Mira variables are
known to be in the AGB phase (e.g., Höfner & Olofsson 2018).

Depending on the chemical composition of the photosphere or the outer dust
envelope, AGB stars are classified into O-rich or C-rich stars. Amorphous
silicates in the envelopes around O-rich AGB (OAGB) stars produce 10 $\mu$m
and 18 $\mu$m features in emission or absorption (e.g., Suh 1999). Dust
envelopes around C-rich AGB (CAGB) stars consist of amorphous carbon (AMC) and
SiC dust (e.g., Suh 2000). SiC dust produces the 11.3 $\mu$m emission feature.

Suh (2021) presented a catalog of 11,209 OAGB stars and 7172 CAGB stars in our
Galaxy and presented various IR two-color diagrams (2CDs) for the AGB stars.
The IR 2CDs have been useful in studying the chemical and physical properties
of various celestial objects (e.g., van der Veen & Habing 1988; Suh 2015).

Thanks to the optical gravitational lensing experiment IV (OGLE4) project, the
number of known Mira variable in our Galaxy has increased significantly.
Iwanek et al. (2022) presented a sample of 65,981 OGLE4 Mira variables in our
Galaxy, from which the majority (47,532 objects) comprises new discoveries.

In this work, we investigate IR properties of OGLE4 Mira variables in our
Galaxy. For each object, we cross-identify the Wide-field Infrared Survey
Explorer (WISE), Two-Micron All-Sky Survey (2MASS), and Infrared Astronomical
Satellite (IRAS) counterparts. We present period-magnitude and period-color
relations for the Mira variables. We also revise and update the previous
catalog of AGB stars in our Galaxy using the new sample of OGLE4 Mira
variables. For the new list of AGB stars, we cross-identify the IRAS, AKARI,
Midcourse Space Experiment (MSX), AllWISE, 2MASS, OGLE4, Gaia Data Release 3
(DR3), and American Association of Variable Star Observers (AAVSO;
international variable star index; version 2022 September 23; Watson et al.
2022) counterparts and present various IR 2CDs. We compare the observations
with theoretical models and discuss general properties of the AGB stars.

Table 1: OGLE4 Mira variables in our Galaxy

OGLE4M1 | AllWISE2 | 2MASS2 | Gaia2,3 | IRAS4 |  |
---|---|---|---|---|---|---
65,981 | 64,962 | 65,931 | 57,482 | 16,379 |  |

1OGLE4 Mira variables from Iwanek et al. (2022). 2The number of counterparts
within 5′′ (2′′ for Gaia; duplicate sources are excluded). 3 Gaia DR3 LPV
sources (Lebzelter et al. 2022). 4The number of positive counterparts of IRAS
PSC sources within 60′′ (see Section 2.7).

Figure 1: WISE-2MASS 2CDs for color-selected OAGB and CAGB stars (the upper
panel; see Section 2.4) and for OAGB and CAGB stars classified from the Gaia
RP spectra (the lowe panel; see Section 2.5) for the OGLE4 Mira variables. The
2CDs also show known AGB stars from Suh (2021) (see Table 3). The observations
are compared with theoretical models (see Section 4.1). The brown line in the
upper-panel shows the boundary line between OAGB and CAGB stars for the color-
selection method. The number of objects and the number of plotted objects (in
parenthesis) with good-quality observed data are shown for each class.  Figure
2: Number density distributions of the Gaia magnitudes at the BP band for OAGB
and CAGB stars classified from the Gaia RP spectra for the OGLE4 Mira
variables (see Section 2.5). Figure 3: Period-Magnitude relations for the
OGLE4 Mira variables from Iwanek et al. (2022) with known distances from Gaia
DR3 (see Section 2.6). Figure 4: Period-color relations for the OGLE4 Mira
variables from Iwanek et al. (2022).

##  2 OGLE4 Mira variables in our Galaxy

Iwanek et al. (2022) presented a sample of 65,981 OGLE4 Mira variables in our
Galaxy (40,356 objects in the Galactic bulge fields and 25,625 objects in the
Galactic disk).

###  2.1 Infrared Photometric Data

The IRAS point source catalog (PSC) (Beichman et al., 1988) provided
photometric data at four bands (12, 25, 60, and 100 $\mu$m) for 245,889
objects. The IRAS’s photometric sensitivity (averaged 10 sigma value at the 12
$\mu$m band) is 0.7 Jy. MSX (Egan et al. 2003) conducted the Galactic plane
survey with a higher sensitivity (about 0.17 Jy at the 8.28 $\mu$m band) at
four MIR bands (8.28, 12.13, 14.65, and 21.34 $\mu$m) for 454,091 objects
(MSXC6). AKARI (Murakami et al., 2007) provided the PSC data at two bands (9
and 18 $\mu$m; Ishihara et al. 2010) for 870,973 objects and bright-source
catalogue (BSC) data at four bands (65, 90, 140, and 160 $\mu$m; Yamamura et
al. 2010) for 427,071 objects. The field of view (FOV) pixel sizes of the
IRAS, MSX, AKARI PSC, and AKARI BSC images are
0$.\mkern-4.0mu^{\prime}$75x(4$.\mkern-4.0mu^{\prime}$5-4$.\mkern-4.0mu^{\prime}$6),
18$.\\!\\!^{\prime\prime}$3, 10′′, and 30′′, respectively.

Two-Micron All-Sky Survey (2MASS) (Cutri et al., 2003) provided PSC data at J
(1.25 $\mu$m), H (1.65 $\mu$m), and K (2.16 $\mu$m) bands for more than 470
million objects. The FOV pixel size of the 2MASS image is 2′′. WISE (Wright et
al., 2010) conducted the entire sky survey at W1 (3.4 $\mu$m), W2 (4.6
$\mu$m), W3 (12 $\mu$m), and W4 (22 $\mu$m) bands. For the four WISE bands,
the FOV pixel sizes are 2$.\\!\\!^{\prime\prime}$75,
2$.\\!\\!^{\prime\prime}$75, 2$.\\!\\!^{\prime\prime}$75, and
5$.\\!\\!^{\prime\prime}$5, and the 5$\sigma$ photometric sensitivities are
0.068, 0.098, 0.86, and 5.4 mJy. In 2013, the AllWISE source catalog released
to present more accurate positions and four-band fluxes for over 747 million
objects.

The FOV pixel size of the CCD camera system for OGLE4 is
0$.\\!\\!^{\prime\prime}$26 (see https://ogle.astrouw.edu.pl/).

###  2.2 Cross-matches

Table 1 lists the numbers of the cross-identified AllWISE, 2MASS, Gaia, and
IRAS PSC counterparts for the 65,981 OGLE4 Mira variables in our Galaxy
(Iwanek et al. 2022). For each object, we find the nearest cross-matched
source within 5′′ (2′′ for Gaia). See Section 2.7 for the IRAS counterparts.

In some cases, multiple OGLE4 objects may have the same cross-matched source
because the telescope for OGLE4 has a smaller FOV pixel size. So we have
checked all of the duplicated cross-matches and selected only the nearest
object for the one source in the AllWISE, 2MASS, Gaia, and IRAS PSC. We find
that there are 64,962 OGLE4 Mira variables that have AllWISE counterparts.

We need to confirm that the cross-matched sources show the expected properties
(brightness, colors, etc) of the original object. A diagram comparing the
fluxes of the identified source and original source is useful. In most cases,
they show a good correlation (e.g., see Figure 1 in Suh 2021), but some
objects show large deviations. A reason for the deviations could be because
they are different objects. In this work, we have reduced the number of the
deviated objects by using the multiple cross-match and comparison processes
(See Section 2.7). Because another reason for the deviations could be the
large scale variations of AGB stars or different beam sizes of the two
telescope systems, we do not exclude the remaining objects that show the large
deviations.

IR 2CDs that show a general character of the object is also useful. For OGLE4
Mira variables, the IR 2CDs using the cross-matched sources need to show that
the new objects produce expected colors similar to previously known Mira
variables or AGB stars (Section 2.4).

In this paper, we use only good-quality observational data for all wavelength
bands for the IRAS, 2MASS, WISE, AKARI, and MSX photometric data (q = 3 for
IRAS and AKARI; q = A for 2MASS; q = A or B for WISE; q = 3 or 4 for MSX).

###  2.3 Color indices

The color index is defined by

$M_{\lambda 1}-M_{\lambda 2}=-2.5\log_{10}{{F_{\lambda 1}/ZMF_{\lambda
1}}\over{F_{\lambda 2}/ZMF_{\lambda 2}}}$ (1)

where $ZMF_{\lambda i}$ is the zero magnitude flux (ZMF) at given wavelength
($\lambda i$). $ZMF_{\lambda i}$ is specified in the manual for the
corresponding telescope system (see Table 4 in Suh 2021).

Figures 1, 4, and 5 show IR 2CDs using various IR color indices. Note that we
use only good quality observational data for all wavelength bands in plotting
the IR 2CDs (see Section 2.2).

To consider the Galactic extinction processes, we plot reddening vectors for
IR 2CD using NIR data (see Suh 2021 for details of the reddening vectors).

###  2.4 The 2MASS-WISE 2CD

IR 2CDs have been useful to distinguish different characteristics of various
heavenly bodies (e.g., van der Veen & Habing 1988; Suh 2018). Though the
classification using the IR spectra would be much more reliable, the IR color-
selection method is useful because the method can make right classification
for a major portion of known OAGB or CAGB stars with IR spectroscopic data
(Suh & Hong 2017; Suh 2021). Therefore, the IR color-selection method using an
IR 2CD can be useful to distinguish between OAGB and CAGB for AGB stars
without reliable spectroscopic data.

The WISE-2MASS 2CD using W3[12]$-$W4[22] versus K[2.2]$-$W3[12] looks to be
more useful than any other 2CDs because the color-combination can distinguish
different classes of AGB stars more correctly (Suh 2018; Suh 2020; Suh 2021;
see Section 3.1) and the good quality observed data (for the two IR colors)
are available for a larger number of sample objects. The brown line on the
upper panel of Figure 1 shows the boundary line between OAGB and CAGB stars
obtained for known AGB stars (Suh 2021), which is used for the color-selection
method in this paper.

The upper panel of Figure 1 shows the WISE-2MASS 2CD for the OGLE4 Mira
variables (see Table 1) as well as the AGB stars in the catalog of Suh (2021)
(see Table 3). On this 2CD, we can find new color-selected AGB stars from the
65,981 OGLE4 Mira variables: 59,928 OAGB and 2674 CAGB stars. Because another
3379 objects do not have the good-quality observed color indices for the 2CD,
the objects can not be classified using the color-selection method.

###  2.5 Gaia RP spectra

The Gaia DR3 LPV (Lebzelter et al. 2022) data provide the low-resolution RP
spectra at visual wavelength bands to identify C-rich stars. Using the Gaia RP
spectral data for the cross-identified 57,482 OGLE4 Mira variables (see Table
1), they can be classified into 56,596 O-rich and 886 C-rich stars by using
the selection criteria (7 $<$ median_delta_wl_rp $<$ 11 and GBP $<$ 19 mag)
presented in Section 2.4 in Lebzelter et al. (2022). The lower panel of Figure
1 shows the WISE-2MASS 2CD for OAGB and CAGB stars classified from the Gaia RP
spectra for the OGLE4 Mira variables.

Figure 2 shows number density distributions of the Gaia magnitudes at the BP
band for OAGB and CAGB stars classified from the Gaia RP spectra for the OGLE4
Mira variables. The number density in the histogram denotes the raw count
divided by the total number of counts and the bin width, so that the area
under the histogram integrates to one. We find that most CAGB stars classified
from the Gaia RP spectra are bright nearby AGB stars.

We find that the classification using the Gaia RP spectra is right for a
portion of known AGB stars with IR spectroscopic data. Compared with the
color-selection method (see Section 2.4), we find that the accuracy ratio of
the RP spectral method is lower (see Section 3.1). This could be due to severe
interstellar extinction processes at the visual bands or the effect of dust
envelopes around AGB stars. The RP spectral data would be more useful for
studying nearby AGB stars with thin dust envelopes.

In this work, we do not use the Gaia RP spectral data for classifying the
OGLE4 Mira variables because we believe that the color-selection method (see
Section 2.4) is more useful for the large sample of AGB stars in our Galaxy.
However, the catalog data for all of the sample AGB stars (see Section 3) will
include the information from the available Gaia RP spectral data.

###  2.6 Period-magnitude and period-color relations

Figure 3 shows period-magnitude relations at W2[4.6] and W3[12] bands for the
OGLE4 Mira variables from Iwanek et al. (2022) with known distances from Gaia
DR3. For each object, the absolute magnitudes at the two WISE bands are
calculated from the corrected distances from Bailer-Jones et al. (2021) for
the cross-matched Gaia DR3 LPV object (Lebzelter et al. 2022). Generally, we
find that Mira variables with longer periods show brighter IR fluxes. The
large deviations from the linear relations would be mainly due to uncertainty
in the Gaia distances. As expected for Mira variables (e.g., Soszyński et al.
2013; Suh 2020), we find that the OGLE4 Mira variables occupy a single
sequence on the period-magnitude diagrams.

Figure 4 shows period-color relations for 59,928 OAGB and 2674 CAGB stars
classified from the 65,981 OGLE4 Mira variables (see Section 2.4). Generally,
we find that Mira variables with longer periods show redder IR colors. Mira
variables with longer periods are expected to have thicker dust envelopes as
they are generally more evolved AGB stars.

###  2.7 IRAS PSC sources

Because IRAS has very low angular resolution, the AllWISE or 2MASS
counterparts obtained from the positions of cross-matched AKARI or MSX sources
(with a higher angular resolution; see Section 2.1) would be much more
reliable.

For all of the 245,889 IRAS PSC sources, we have found the AKARI, MSX,
AllWISE, 2MASS, OGLE4, Gaia DR3, and AAVSO counterparts by using following
method. We cross-identify the AKARI PSC, MSX, and AKARI BSC counterpart by
finding the nearest source within 60′′ using the position given in the IRAS
PSC (version 2.1). Then, we cross-identify the AllWISE, 2MASS, and OGLE4
counterparts using the best position of the available AKARI PSC, MSX, or AKARI
BSC counterpart. Only when there is no AKARI or MSX counterpart, we use the
position of the IRAS PSC.

The similar method was used in Suh (2018) and Suh (2021). But in this work, we
use the MSX data in addition to the AKARI data to improve the accuracy of the
cross-match (see Section 2.2).

From the 65,981 OGLE4 Mira variables, we find 18,413 IRAS counterparts
(excluding duplicate sources) within 60′′. From the 18,426 sources, we find
that 16,379 sources are positive (or reliable) counterparts, for which the
AllWISE counterpart of the IRAS source is the same as the AllWISE counterpart
of the OGLE4M object (see Table 2).

##  3 A new sample of AGB stars

Suh (2021) presented a catalog AGB stars in our Galaxy in two parts: one is
based on the IRAS PSC and the other one is based on the AllWISE source
catalog. In the IRAS based catalog, there are 5908 OAGB-IRAS and 3596 CAGB-
IRAS objects. And in the AllWISE based catalog, there are 5301 OAGB-WISE and
3576 CAGB-WISE objects. Because IRAS has lower angular resolution and weaker
sensitivity (see Section 2.1) than WISE, the AGB-IRAS objects are generally
brighter than the AGB-IRAS objects.

In this work, we revise and update the catalog using the new sample of OGLE4
Mira variables presented by Iwanek et al. (2022). We use the color-selection
method (see Section 2.4) to classify the OGLE4 Mira variables into OAGB or
CAGB stars, except for the objects for which more reliable IR spectroscopic
data are available.

We present a new catalog of 19,196 OAGB-IRAS, 4118 CAGB-IRAS, 45,413 OAGB-
WISE, and 5366 CAGB-WISE objects (see Tables 3 and 4). For the new list of AGB
stars, we cross-identify the IRAS, AKARI, MSX, AllWISE, 2MASS, OGLE4, Gaia
DR3, and AAVSO counterparts and present the new catalog data, which contain
all of the available information.

Among the 65,981 OGLE4 Mira variables from Iwanek et al. (2022), 62,989
objects (15,248 OAGB-IRAS, 797 CAGB-IRAS, 44,977 OAGB-WISE, and 1967 CAGB-
WISE) are included in the new catalog of OAGB and CAGB stars and 2000 objects
are classified into NI-O4 or NW-O4 subgroups (see Tables 3 and 4). And 992
objects are not included in the catalog because 984 objects do not have the
AllWISE counterparts and can not be cross-identified as any known AGB stars
and another 8 objects can be cross-identified as known S type or silicate-
carbon stars.

Table 2: IRAS PSC sources

IRAS PSC1 | AKARI PSC2 | MSX2 | AKARI BSC2 | AllWISE3 | 2MASS3 | OGLE4M3,4 | Gaia3,5
---|---|---|---|---|---|---|---
245,889 | 159,031 | 82,001 | 61,705 | 245,756 | 245,049 | 16,379 | 98,522

1version 2.1. 2The number of counterparts within 60′′ using the IRAS position
(duplicate sources are exlcuded). 3The number of counterparts within 60′′
using the best position from IRAS, AKARI, and MSX (duplicate sources are
exlcuded). 4The number of positive counterparts within 60′′ (see Section 2.7).
5Gaia DR3 LPV object (Lebzelter et al. 2022).

Table 3: Sample of AGB stars based on the IRAS PSC (AGB-IRAS)

Class | Subgroup | Reference | Number | Selected | Duplicate | AE1 | Remaining
---|---|---|---|---|---|---|---
OAGB-IRAS | OAGB-IRAS(2021) | Suh (2021) | 5908 | 5908 | 0 | -12 | 5907
OAGB-IRAS | OI-O4S | This work | 16,3793 | 22734 | 6465 | 0 | 1627
OAGB-IRAS | OI-O4C | This work | 16,3793 | 12,9206 | 12585 | 0 | 11,662
OAGB-IRAS | OAGB-IRAS(2022) | Total | - | - | - | - | 19,196
CAGB-IRAS | CAGB-IRAS(2021) | Suh (2021) | 3596 | 3596 | 0 | -17 | 3595
CAGB-IRAS | CI-O4S | This work | 16,3793 | 2958 | 2439 | 0 | 52
CAGB-IRAS | CI-O4C | This work | 16,3793 | 4936 | 2210 | 0 | 471
CAGB-IRAS | CAGB-IRAS(2022) | Total | - | - | - | - | 4118
No class | NI-O4 | This work | 16,3793 | 3986 | 375 | - | 361

1added or excluded. 2IRAS 17056-3959 is a carbon star (now in CI-O4S). 3the
number (16379) of positive IRAS counterparts from all of the OGLE4 Mira
variables (see Section 2.7). 4objects with O-rich spectra from IRAS LRS and/or
SIMBAD. 5in OAGB-IRAS(2021) or CAGB-IRAS(2021) or in the list of S stars or
silicate-carbon stars in AGB-IRAS (Suh 2021). 6color-selected OAGB, CAGB, or
No class objects (see Section 2.4). 7IRAS 17514-3354 is an OAGB star (now in
OI-O4C). 8objects with C-rich spectra from SIMBAD. 9in CI-2021. 10in CI-2021
or OI-2021.

Table 4: Sample of AGB stars based on the AllWISE source catalog (AGB-WISE)

Class | Subgroup | Reference | Number | Selected | Duplicate | Remaining
---|---|---|---|---|---|---
OAGB-WISE | OAGB-WISE(2021) | Suh (2021) | 5301 | 5301 | 41+492 | 5248
OAGB-WISE | OW-O4C | This work | 64,9623 | 59,9284 | 14,9365+48276 | 40,165
OAGB-WISE | OAGB-WISE(2022) | Total | - | - | - | 45,413
CAGB-WISE | CAGB-WISE(2021) | Suh (2021) | 3576 | 3576 | 27+28 | 3572
CAGB-WISE | CW-O4C | This work | 64,9623 | 26744 | 7289+1526 | 1794
CAGB-WISE | CAGB-WISE(2022) | Total | - | - | - | 5366
No class | NW-O4 | This work | 64,9623 | 23604 | 71510+66 | 1639

1in OAGB-IRAS(2021). 2in OI-O4S or OI-O4C. 3the number of OGLE4 Mira variables
with AllWISE counterparts (see Table 1). 4color-selected OAGB, CAGB, or No
class objects (see Section 2.4). 5in OI-O4S, OI-O4C, or CI-O4S. 6in OAGB-
WISE(2021) or CAGB-WISE(2021). 7in CAGB-IRAS(2021) or OAGB-IRAS(2021). 8in
CI-O4S or CI-O4C. 9in CI-O4S, CI-O4C, or OI-O4S. 10in NI-O4.

###  3.1 AGB stars based on the IRAS PSC

Suh (2021) presented a catalog AGB stars in our Galaxy that is based on the
IRAS PSC: 5908 OAGB-IRAS and 3596 CAGB-IRAS objects. We revise and update the
catalog using the new sample of OGLE4 Mira variables presented by Iwanek et
al. (2022).

In identifying important dust features of AGB stars, IRAS Low Resolution
Spectrograph (LRS; $\lambda$ = 8$-$22 $\mu$m) data are useful (Kwok et al.
1997). For OAGB stars, the IRAS LRS class E objects show the 10 $\mu$m
silicate feature in emission and class A objects show the 10 $\mu$m silicate
feature in absorption. For CAGB stars, the LRS class C objects show the 11.3
$\mu$m SiC feature in emission.

For CAGB stars in the IRAS LRS class C, the accuracy ratio of the color-
selection method (see Section 2.4) is 82 %, while the accuracy ratio of the
Gaia RP spectral method (see Section 2.5) is 55 %. For OAGB stars in the IRAS
LRS class E or A, the inaccuracy ratios of the two methods are about the same.

Table 3 lists a new catalog of AGB-IRAS objects. We correct some misidentified
objects from Suh (2021). We add the OGLE4 Mira variables (presented by Iwanek
et al. 2022) with the positive IRAS counterparts (see Section 2.7), which are
classified into OAGB and CAGB stars using the color-selection method using the
WISE-2MASS 2CD (see Section 2.4).

In classifying into OAGB or CAGB, the more reliable method using the spectral
information (the spectral type from the IRAS LRS or SIMBAD) gets priority over
the color-selection method using IR 2CDs.

Finally, we present a new list of 23,314 AGB-IRAS objects, which are
classified into 19,196 OAGB-IRAS and 4118 CAGB-IRAS stars.

###  3.2 AGB stars based on the AllWISE source catalog

Suh (2021) presented a catalog AGB stars in our Galaxy that is based on the
AllWISE source catalog: 5301 OAGB-WISE and 3576 CAGB-WISE objects. We revise
and update the catalog using the new sample of OGLE4 Mira variables presented
by Iwanek et al. (2022).

Table 4 lists a new catalog of AGB-WISE objects. They are identified based on
the AllWISE source catalog and do not have positive IRAS counterparts. We make
corrections for some misidentified AGB-WISE objects from Suh (2021). We add
the OGLE4 Mira variables (presented by Iwanek et al. 2022) with the AllWISE
counterparts.

Now, we present a new list of 50,779 AGB-WISE objects, which are classified
into 45,413 OAGB-WISE and 5366 CAGB-WISE stars.

Figure 5: IRAS-2MASS 2CDs for AGB-IRAS stars in the previous catalog of Suh
(2021) and new AGB-IRAS stars identified from the OGLE4 Mira variables (see
Tables 3) compared with theoretical models (see Section 4.1). The number of
objects and the number of plotted objects (in parenthesis) with good-quality
observed data are shown for each class.

Figure 6: WISE-2MASS 2CDs for all AGB stars in the previous catalog of Suh
(2021) and new AGB stars identified from the OGLE4 Mira variables (see Tables
3 and 4) compared with theoretical models (see Section 4.1). The number of
objects and the number of plotted objects (in parenthesis) with good-quality
observed data are shown for each class.  Figure 7: Averaged observed IR colors
for different classes of AGB stars (see Tables 3 and 4; Figures 4 and 5). The
number of objects and the number of considered objects (in parenthesis) with
good-quality observed data for the horizontal axis are shown for each class.
Figure 8: Number density distributions of observed K[2.2]$-$W3[12] indices for
AGB-IRAS and AGB-WISE stars. The number of objects and the number of
considered objects (in parenthesis) with good-quality observed data are shown
for each class.  Figure 9: Number density distributions of the observed OGLE4
periods for AGB-IRAS and AGB-WISE stars.

##  4 IR 2CDs and Discussion

Figures 5 and 6 show various 2CDs for AGB stars in the new revised and
extended catalog (see Section 3; Tables 3 and 4). On all of the 2CDs in
Figures 1, 5, and 6, theoretical models for AGB stars are also plotted to be
compared with the observations.

###  4.1 Theoretical Dust Shell Models

In this paper, we use the same theoretical dust shell models as the ones used
in Suh (2020) and Suh (2021) (see Section 3 and Table 7 in Suh (2021) for the
summary). We use the radiative transfer code DUSTY (Ivezić & Elitzur 1997),
which assumes a spherically symmetric dust shell around a central star. We
assume that the central stars is a black body (T = 3000-2000 K).

For all models, we use a simple power law ($\rho\propto r^{-2}$) for the dust
density distribution and assume that the dust formation temperature ($T_{c}$)
is 1000-500 K. The outer radius of the dust shell is taken to be $10^{4}$
times the inner radius, the fiducial wavelength of the dust shell optical
depth ($\tau_{10}$) is 10 $\mu$m, and the uniform radius of spherical dust
grains is 0.1 $\mu$m.

For OAGB stars, we use dust optical constants of warm and cold silicates (Suh
1999), amorphous alumina (Suh 2016), and Fe0.9Mg0.1O (Henning et al. 1995).
For CAGB stars, we use dust optical constants of AMC (Suh 2000), SiC (Pégourié
1988), and Mg0.9Fe0.1S (Begemann et al. 1994).

On the 2CDs in Figures 1, 4, and 5, theoretical models for AGB stars are
plotted. For OAGB models (silicate $T_{c}$ = 1000 K), the model results for
$\tau_{10}$ = 0.001, 0.01, 0.05, 0.1, 0.5, 1, 3, 7, 15, 30, and 40 are plotted
from left to right on the 2CDs. For CAGB models (AMC $T_{c}$ = 1000 K), the
model results for $\tau_{10}$ = 0.001, 0.01, 0.1, 0.5, 1, 2, 3, and 5 are
plotted from left to right on the 2CDs.

For physical implications of the theoretical dust shell models, refer to Suh
(2020) and Suh (2021). Suh (2020) investigated IR properties of AGB Stars in
Our Galaxy and the Magellanic Clouds and discussed the physical implications
(e.g., mass-loss rates, masses, ages, and metallicities) of the theoretical
model parameters for different classes of AGB stars. Those model parameters
were obtained from the comparison of the theoretical dust shell models with
the observations on various IR 2CDs.

###  4.2 Comparison between Theory and Observations

On various IR 2CDs in Figures 5 and 6, we compare the observations with the
theoretical dust shell models for AGB stars. Generally, the theoretical dust
shell models for OAGB and CAGB stars can reproduce the observed points fairly
well.

The upper panel of Figure 5 plots an IRAS 2CD using [25]$-$[60] versus
[12]$-$[25]. A group of CAGB stars in the upper-left region are visual carbon
stars that show excessive flux at 60 $\mu$m due to the remnant of the earlier
phase when they were OAGB stars (e.g. Chan & Kwok 1990). Another group of CAGB
stars in the lower region, which extends to the right side, are IR carbon
stars. The IR carbon stars with redder [12]$-$[25] colors are believed to be
more evolved or massive stars with thicker dust envelope (see Section 1 in Suh
2020 for origin and evolution of carbon stars).

The lower panel of Figure 5 plots an 2MASS-IRAS 2CD using IR[12]$-$IR[60]
versus K[2.2]$-$IR[12]. On this 2CD, the distinguishment between OAGB and CAGB
stars is clearer. Visual carbon stars are located in the lower-left region on
this 2CD, because they show bluer K[2.2]$-$IR[12] and IR[12]$-$IR[60] colors
than IR carbon stars.

Figure 6 shows 2CDs using WISE and 2MASS colors. The upper panel of Figure 6
shows a WISE 2CDs using W3[12]$-$W4[22] versus W2[4.6]$-$W3[12]. The lower
panel of Figure 6 shows a WISE-2MASS 2CD using W3[12]$-$W4[22] versus
K[2.2]$-$W3[12]. On these 2CDs, the theoretical dust shell models for CAGB
stars with various C-rich dust grains can reproduce a wider range of observed
W3[12]$-$W4[22] colors.

In the lower regions of the 2CDs in Figure 6, we find that there are many
newly identified CAGB-IRAS or CAGB-WISE objects in this work, which are color-
selected from the OGLE4 Mira variables (see Section 2.4). We may need to
confirm the C-rich nature of these objects by using a more reliable IR
spectroscopic method. The objects in the upper-left regions on the 2CDs in
Figure 6 are likely to be early stage OAGB stars or visual carbon stars with
thin and detached dust shells ($T_{c}$ = 200-300 K; see Suh 2020).

###  4.3 Statistical properties of AGB stars

Figure 7 shows the error bar plots of the observed colors on the various 2CDs
used for all classes of AGB stars presented in Figures 5 and 6. Because IRAS
LRS class A objects (mostly OH/IR stars) and class C objects (IR carbon stars)
have thicker dust envelopes, they are believed to be more evolved or massive
AGB stars. We find that these more evolved OAGB or CAGB objects with thicker
dust envelopes are show redder colors for all of the IR color indices except
for the [25]$-$[60] color. This is because visual carbon stars with thin dust
envelopes show excessive flux at 60 $\mu$m. Generally, AGB-WISE objects show
bluer colors compared with AGB-IRAS objects. This would be because most AGB-
WISE objects are less evolved and have thinner dust shells than AGB-IRAS
objects.

Figure 8 shows number density distributions of K[2.2]$-$W3[12] indices for
AGB-IRAS and AGB-WISE stars. Compared with AGB-WISE objects, more AGB-IRAS
objects show redder colors because they have thicker dust envelopes. Compared
with CAGB-IRAS objects, OAGB-IRAS objects show redder colors.

Figure 9 shows number density distributions of the OGLE4 periods for AGB-IRAS
and AGB-WISE stars. We find that more AGB-IRAS objects show longer averaged
pulsation periods than AGB-WISE objects, which are generally less evolved or
less massive. We also find that CAGB-IRAS stars show generally longer periods
than OAGB-IRAS stars. This would be due to a selection effect of OGLE4
observations. Because the OGLE4 periods are obtained from the $I$ filter
photometry, most OAGB-IRAS stars with thick dust envelopes and longer
pulsation periods (OH/IR stars) are not observed by the OGLE4 observations.

##  5 Summary

We have investigated infrared properties of OGLE4 Mira variables in our
Galaxy. For each object, we have cross-identified the AllWISE, 2MASS, and IRAS
counterparts. We have presented various infrared two-color diagrams (2CDs) and
period-magnitude and period-color relations for the Mira variables. Generally,
the Mira variables with longer periods are brighter in the IR fluxes and
redder in IR colors.

We have revised and updated the catalog of AGB stars in our Galaxy from Suh
(2021) using the new sample of OGLE4 Mira variables. Now, we present a new
catalog of 74,093 (64,609 OAGB and 9484 CAGB) stars in our Galaxy. A group of
23,314 (19,196 OAGB-IRAS and 4118 CAGB-IRAS) stars are identified based on the
IRAS PSC and another group of 50,779 (45,413 OAGB-WISE and 5366 CAGB-WISE)
stars are identified based on the AllWISE source catalog.

For all of the AGB stars in the new catalog, we have cross-identified the
IRAS, AKARI, MSX, AllWISE, 2MASS, OGLE4, Gaia DR3, and AAVSO counterparts and
presented a new catalog data, which contains all of the available information.

We have presented various IR 2CDs for all AGB stars in the catalog. Comparing
the observations with the theory on the 2CDs, we have found that basic
theoretical dust shell models can account for the IR observations fairly well
for most of the AGB stars.

We have found that the AGB-IRAS objects show longer averaged pulsation periods
and redder averaged IR colors than AGB-WISE objects, which are generally less
evolved or less massive.

The new revised and updated data for the catalog of AGB stars in our Galaxy
(version 2022) will be accessible through the author’s webpage
(https://web.chungbuk.ac.kr/~kwsuh/agb.htm).

###### Acknowledgements.

This research was supported by Basic Science Research Program through the
National Research Foundation of Korea (NRF) funded by the Ministry of
Education (2022R1I1A3055131). This research has made use of the VizieR
catalogue access tool and the SIMBAD database, operated at CDS, Strasbourg,
France. This research has made use of the NASA/IPAC Infrared Science Archive,
which is operated by the Jet Propulsion Laboratory, California Institute of
Technology, under contract with the National Aeronautics and Space
Administration.

## References

  * Bailer-Jones et al. (2021) Bailer-Jones, C., Rybizki, J., Fouesneau, M., et al. 2021,Estimating Distances from Parallaxes. V. Geometric and Photogeometric Distances to 1.47 Billion Stars in Gaia Early Data Release 3, AJ, 161, 147
  * Begemann et al. (1994) Begemann, B., Dorschner, J., Henning, T., Mutschke, H., & Thamm, E. 1994, A Laboratory Approach to the Interstellar Sulfide Dust Problem, ApJ, 423, L71
  * Begemann et al. (1997) Begemann, B., Dorschner, J., Henning, T., et al. 1997, Aluminum Oxide and the Opacity of Oxygen-Rich Circumstellar Dust in the 12-17 Micron Range, ApJ, 476, 199
  * Beichman et al. (1988) Beichman, C. A., Neugebauer, G., Habing H., Clegg, P. E., & Chester, T. C. 1988, IRAS Catalogs and Atlases: Explanatory Supplement, NASA RP-1190 (Washington: NASA)
  * Chan & Kwok (1990) Chan, S. J., & Kwok, S., 1990, Evolution of Infrared Carbon Stars, A&A, 237, 354
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., Van Dyk, S., et al. 2003, VizieR Online Data Catalog: 2MASS All-Sky Catalog of Point Sources
  * Egan et al. (2003) Egan, M. P., Price, S. D., & Kraemer, K. E. 2003, The Midcourse Space Experiment Point Source Catalog v2.3 (2003 October) Air Force Research Lab. Tech. Rep. AFRL-VS-TR-2003-1589 (Springfield, VA: NTIS)
  * Henning et al. (1995) Henning, T., Begemann, B., Mutschke, H., & Dorschner, J. 1995, Optical Properties Of Oxide Dust Grains, A&AS, 112, 143
  * Höfner & Olofsson (2018) Höfner, S., & Olofsson, H., 2018, Mass loss of stars on the asymptotic giant branch: Mechanisms, models and measurements, A&A Rev., 26, 1
  * Ishihara et al. (2010) Ishihara, D., Acke, B., Verhoelst, T., et al. 2010, The AKARI/IRC Mid-Infrared All-Sky Survey, A&A, 514, A1
  * Iwanek et al. (2022) Iwanek, P., Soszyński, I.. Kozłowski, S. et al. 2022, The OGLE Collection of Variable Stars: Nearly 66,000 Mira Stars in the Milky Way, ApJS, 260, 46
  * Ivezić & Elitzur (1997) Ivezić, A., & Elitzur, M. 1997, Self-Similarity and Scaling Behaviour of Infrared Emission from Radiatively Heated Dust - I. Theory, MNRAS, 287, 799
  * Kwok et al. (1997) Kwok, S., Volk, K., & Bidelman, W. P. 1997, Classification and Identification of IRAS Sources with Low-Resolution Spectra, ApJS, 112, 557
  * Lebzelter et al. (2022) Lebzelter, T., Mowlavi, N., Lecoeur-Taibi, I., et al. 2022, Gaia Data Release 3: The second Gaia catalogue of Long-Period Variable candidates, A&A, in press
  * Murakami et al. (2007) Murakami, H., Baba, H., Barthel, P., et al. 2007, The Infrared Astronomical Mission AKARI, PASJ, 59, 369
  * Pégourié (1988) Pégourié, B. 1988, Optical Properties of Alpha Silicon Carbide, A&A, 194, 335
  * Soszyński et al. (2013) Soszyński, I., Wood, P. R., Udalski, A. 2013, ApJ, 779, 167
  * Suh (1999) Suh, K.-W. 1999, Optical Properties of the Silicate Dust Grains in the Envelopes around Asymptotic Giant Branch Stars, MNRAS, 304, 389
  * Suh (2000) Suh, K.-W. 2000, Optical Properties of the Carbon Dust Grains in the Envelopes around Asymptotic Giant Branch Stars, MNRAS, 315, 740
  * Suh (2015) Suh, K.-W. 2015, Infrared Two-Color Diagrams for AGB Stars, Post-AGB Stars, and Planetary Nebulae, ApJ, 808, 165
  * Suh (2016) Suh, K.-W. 2016, Optical Properties of Amorphous Alumina Dust in the Envelopes Around O-Rich AGB Stars, JKAS, 49, 127
  * Suh (2018) Suh, K.-W. 2018, Infrared Two-Color Diagrams of AGB stars and Planetary Nebulae using WISE data, JKAS, 51, 155
  * Suh (2020) Suh, K.-W. 2020, Infrared Properties of Asymptotic Giant Branch Stars in Our Galaxy and the Magellanic Clouds, ApJ, 891, 43
  * Suh (2021) Suh, K.-W. 2021, A New Catalog of Asymptotic Giant Branch Stars in Our Galaxy, ApJS, 256, 43
  * Suh & Hong (2017) Suh, K.-W., & Hong, J. 2017, A New Catalog of AGB Stars Based on Infrared Two-Color Diagrams, JKAS, 50, 131
  * van der Veen & Habing (1988) van der Veen, W. E. C. J., & Habing, H. J. 1988, The IRAS Two-Colour Diagram as a Tool for Studying Late Stages of Stellar Evolution, A&A, 194, 125
  * Watson et al. (2022) Watson, C., Henden, A. A., & Price, A. 2022, Vizier, B/vsx/vsx
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, The Wide-Field Infrared Survey Explorer (WISE): Mission Description And Initial On-Orbit Performance, AJ, 140, 1868
  * Yamamura et al. (2010) Yamamura, I., Makiuti, S., Ikeda, N., et al. 2010, AKARI/FIS All-Sky Survey Bright Source Catalogue Version 1.0 Release Note