# Neutral Utterances are Also Causes: Enhancing Conversational Causal Emotion
Entailment with Social Commonsense Knowledge

Jiangnan Li1,2, Joint work with Pattern Recognition Center, WeChat AI, Tencent
Inc, China.  Fandong Meng3 Zheng Lin1,2, Zheng Lin and Peng Fu are the
corresponding authors.  Rui Liu1,2 Peng Fu1,† Yanan Cao1,2 Weiping Wang1&Jie
Zhou3 1Institute of Information Engineering, Chinese Academy of Sciences,
Beijing, China
2School of Cyber Security, University of Chinese Academy of Sciences, Beijing,
China
3Pattern Recognition Center, WeChat AI, Tencent Inc, China
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Conversational Causal Emotion Entailment aims to detect causal utterances for
a non-neutral targeted utterance from a conversation. In this work, we build
conversations as graphs to overcome implicit contextual modelling of the
original entailment style. Following the previous work, we further introduce
the emotion information into graphs. Emotion information can markedly promote
the detection of causal utterances whose emotion is the same as the targeted
utterance. However, it is still hard to detect causal utterances with
different emotions, especially neutral ones. The reason is that models are
limited in reasoning causal clues and passing them between utterances. To
alleviate this problem, we introduce social commonsense knowledge (CSK) and
propose a Knowledge Enhanced Conversation graph (KEC). KEC propagates the CSK
between two utterances. As not all CSK is emotionally suitable for utterances,
we therefore propose a sentiment-realized knowledge selecting strategy to
filter CSK. To process KEC, we further construct the Knowledge Enhanced
Directed Acyclic Graph networks. Experimental results show that our method
outperforms baselines and infers more causes with different emotions from the
targeted utterance.111 The code and Appendix are available at
https://github.com/LeqsNaN/KEC

##  1 Introduction

Conversation is the most omnipresent medium for people to express daily
emotions. Capturing emotions is critical for machines to understand
conversations, which leads to researchers Majumder et al. (2019); Ghosal et
al. (2019); Shen et al. (2021b) devoting into emotion recognition in
conversations (ERC). However, further reasoning why an emotion is aroused in
conversations is still under explored Poria et al. (2021). Understanding
emotion causes is also a critical step for machines to understand
conversations, which can be further utilized in other high-level tasks like
empathetic generation Kim et al. (2021). Therefore, Poria et al. Poria et al.
(2021) propose a new task called Recognizing Emotion Cause in Conversations
(RECCON). RECCON involves an important subtask called Conversational Causal
Emotion Entailment (C2E2). C2E2 detects causal utterances for a targeted
utterance from conversations where the emotion is given.

Figure 1: A case that the baseline DAG-ERC fails while the selected
commonsense knowledge can help to make a right prediction.

Implied by the name, C2E2 is originally modeled as a “text pair” task by Poria
et al. Poria et al. (2021). In this way, a candidate utterance and the
targeted utterance are paired plus their conversational context and the
emotion. They are concatenated and then processed by a pretrained model.
However, this entailment style ignores the explicit interactions between
utterances, thus contextual information from conversations not fully captured.
To explicitly model the interactions between utterance representations, works
Ghosal et al. (2019); Shen et al. (2021b) in ERC adopt various types of graphs
for conversations, which can achieve decent performance.

As emotion plays a crucial role in C2E2, it is utilized in the entailment
Poria et al. (2021). When conversations are constructed as graphs, it is also
critical to consider emotion information. An intuitive way is to represent
emotion as emotion embedding Liang et al. (2021) and concatenate it with the
node representation. Experiments show that emotion information can remarkably
boost the performance. Specifically, we find that it can dramatically promote
the detection of causal utterances whose emotion is the same as the targeted
utterance. However, causal utterances with different emotions, especially
neutral ones (neutral causal utterances occupy 87$\%$ of this kind of causes),
is still hard to detect even with emotion information. We think the reason
lies in that models are limited in reasoning causal clues and passing them
between utterances. For example, in Fig. 1, the 2nd utterance $u_{2}$ is a
cause to the 3rd utterance $u_{3}$. Models cannot reason why a neutral
utterance about misdialing can evoke sadness to another utterance, leading to
a wrong prediction.

To boost models’ causal reasoning ability, we introduce social commonsense
knowledge (CSK) Hwang et al. (2021) and propose a graph-based structure to
properly utilize CSK. Therefore, causal utterances with different emotions
from the targeted utterance can be further detected. Specifically, social CSK
can infer a person’s social interactions like mental states and reactions.
Back to Fig. 1, equipped with the knowledge, models can infer that $u_{2}$
makes the speaker of $u_{3}$ feel upset and then give a right prediction. To
properly utilize knowledge, we propose a Knowledge Enhanced Conversation graph
(KEC). KEC passes knowledge through edges. As the knowledge should align with
speaker interactions and emotions, we therefore propose a speaker and
sentiment realized knowledge selecting strategy, which can pick up suitable
knowledge for utterances. Furthermore, to process the KEC, we extend the
Directed Acyclic Graph networks (DAG) for ERC Shen et al. (2021b) to build
Knowledge Enhanced DAG networks. To verify our method, we conduct experiments
on the RECCON dataset. Our method significantly outperforms the baseline of
RECCON and other highly related methods. In addition, our method promotes the
detection of causal utterances with different emotions from the targeted
utterance.

Figure 2: The structure of our model. It contains 3 modules: (1) Utterance
Encoder encodes every utterance; (2) KEC graph is constructed from a
conversation and knowledge attached in KEC graph is picked up by the knowledge
selecting strategy. Conversation Encoder then uses Knowledge Enhanced DAG
networks to process KEC graph; (3) Cause Predictor pairs every two utterances
to make predictions.

##  2 Related Work

Conversational Emotion Recognition (ERC) is a highly related task to C2E2. ERC
tends to predict the emotion of utterances in conversations. In ERC, a main
stream of works Majumder et al. (2019); Ghosal et al. (2019); Shen et al.
(2021a) focuses on studying the speaker modelling. Although ERC is a hot-spot
task for researchers to study conversations, further reasoning about emotions
and understanding the causes in conversations is still lack of exploration.

Emotion Cause Extraction (ECE) is to extract the causal clause for an emotion
clause from a document and Gui et al. Gui et al. (2016) construct the most
popular ECE dataset. Ding et al. Ding et al. (2019); Xia et al. Xia et al.
(2019) introduce position embeddings to model causal clauses appeared near the
emotion clause. As emotion and cause are highly related, Turcan et al. Turcan
et al. (2021) propose multi-task learning frameworks for emotion recognition
and cause extraction. To simultaneously extract emotion and cause clauses from
documents, Xia and Ding Xia and Ding (2019) propose a new task called Emotion
Cause Pair Extraction (ECPE) and handle it in a two-stage style. To overcome
the problems brought by two-stage processing, some works Ding et al. (2020a);
Wei et al. (2020); Ding et al. (2020b) propose their novel joint frameworks.
Both ECPE and ECE are only extract causes in articles, which does not extend
to the scene of conversations. Conversations are more complicated due to the
flowing emotion dynamics and speaker participation Ghosal et al. (2019). To
study the causal factor of emotion in conversations, Poria et al. Poria et al.
(2021) propose Recognizing Emotion Causes in Conversations (RECCON). RECCON
involves two subtasks: Conversational Causal Span Extraction (C2SE) and
Conversational Causal Emotion Entailment (C2E2). C2SE focuses on extracting
causal spans from conversations, which is modelled as a machine reading
comprehension task by Poria et al. Poria et al. (2021). In this work, we focus
on C2E2, the classification version of RECCON.

Commonsense Knowledge Utilization (CSK) can bring rational external knowledge
for models. In ERC, some works Zhong et al. (2019); Ghosal et al. (2020); Li
et al. (2021) study how to use CSK to help understand conversation context and
mental states of speakers. In ECE, Gui et al. Yan et al. (2021) propose
Knowledge-Aware Graph using CSK to alleviate position bias brought by position
embeddings. Turcan et al. Turcan et al. (2021) use CSK to enhance ECE and
Emotion Detection. In the task of Sentence Ordering, sequential commonsense
knowledge is structurally utilized by Ghosal et al., Ghosal et al. (2021).
However, for C2E2, it is still under explored how to properly use CSK.

##  3 Method

###  3.1 Task Definition

C2E2 is to find causal utterances for a targeted utterance from the historical
context. Given a conversation $C=[u_{1},...,u_{N}]$, with the corresponding
emotion sequence $E=[e_{1},...,e_{N}]$ and speaker sequence
$P=[p_{1},...,p_{N}]$, every utterance $u_{i}$ is paired with its contextual
utterances $u_{j}$ ($j\leq i$). If $u_{i}$ is a non-neutral utterance and
$u_{j}$ is its causal utterance, the pair $(u_{i},u_{j})$ is labeled with 1.
Otherwise, $(u_{i},u_{j})$ is labeled with 0. In our setting, neutral
utterances are involved as targeted utterances, leading to more negative non-
causal pairs and making the task more challenging.

###  3.2 Utterance Encoder

Similar as other conversation-related tasks (e.g., ERC), utterances in a
conversation are first encoded by an utterance encoder to form utterance
representations. We utilize the pretrained model RoBERTa Liu et al. (2019) to
encode an utterance $u_{i}=[w_{1},...,w_{L_{u}}]$ with $L_{u}$ words. We then
obtain the utterance representation by:

$\displaystyle
s_{i}=\mathrm{Linear}(\mathrm{Maxpooling}(\mathrm{RoBERTa}(u_{i}))),$ (1)

where $s_{i}\in\mathbb{R}^{d_{u}}$ and $d_{u}$ is the dimension of utterance
representation.

###  3.3 Conversation Encoder

As only reviewing conversation history for utterances, the conversation can be
modelled as a directed acyclic graph when not considering self-loop edges.
Based on the Directed Acylic Graph (DAG) networks Shen et al. (2021b), we
introduce our knowledge enhanced conversation graph, which is illustrated in
Fig. 2.

####  3.3.1 Knowledge Enhanced Conversation Graph

Knowledge Enhanced Conversation graph (KEC) is formulated as
$\mathcal{G}=(\mathbb{V},\mathrm{\mathbb{A}_{c}},\mathrm{\mathbb{A}_{k}})$.
$\mathbb{V}$ denotes the set of utterance nodes, $\mathrm{\mathbb{A}_{c}}$
denotes the adjacency matrix of conventional utterance interactions, and
$\mathrm{\mathbb{A}_{k}}$ denotes the adjacency matrix of knowledge passing
between utterances.

##### Utterance nodes.

$\mathbb{V}$ models all utterances as nodes $\mathrm{v}_{i}$ in a
conversation. $\mathrm{v}_{i}$ contains the attribute rep to store the
utterance representation (i.e., $\mathrm{v}_{i}.$rep$=s_{i}$).

##### Utterance interaction adjacency matrix.

Conversational context modeling is crucial to conversation-related tasks,
which can enrich targeted utterances with contextual messages passed from
other utterances. As only the dialogue history is considered, we construct
$\mathrm{\mathbb{A}_{c}}$ as a lower triangle matrix.
$\mathrm{\mathbb{A}_{c}}$ contains two attributes: item and rel, where item
stores 0 or 1 to denote the existence of an edge and rel stores the relation
type of an edge.

For nodes $\mathrm{v}_{i}$, $\mathrm{v}_{j}$,
$\mathrm{\mathbb{A}_{c}}[i,j].\texttt{item}=1$ is present if there is an edge
from $\mathrm{v}_{j}$ to $\mathrm{v}_{i}$, otherwise
$\mathrm{\mathbb{A}_{c}}[i,j].\texttt{item}=0$. For a targeted utterance, it
has been studied that local contextual information is more important than that
in the remote context Shen et al. (2021b). Therefore, a context window with
size of $w_{c}$ is sliding through conversations. In the window, every
contextual utterance sends an edge to the targeted one. Speaker information is
another crucial factor. Speaker interactions can be modelled as two relations:
self-dependency (SD) and inter-speaker dependency (ID) Ghosal et al. (2019).
For $\mathrm{v}_{i}$, $\mathrm{v}_{j}$,
$\mathrm{\mathbb{A}_{c}}[i,j].\texttt{rel}=\mathrm{SD}$ means that the speaker
identity $p_{i}$ is the same as $p_{j}$, otherwise
$\mathrm{\mathbb{A}_{c}}[i,j].\texttt{rel}=\mathrm{ID}$.

##### Knowledge passing adjacency matrix.

We pass the knowledge through edges to enhance the ability of inferring
causes. Similar as $\mathrm{\mathbb{A}_{c}}$, $\mathrm{\mathbb{A}_{k}}$ is
also a lower triangle matrix. $\mathrm{\mathbb{A}_{k}}$ contains two
attributes: item and klg, where item plays the same role as
$\mathrm{\mathbb{A}_{c}}.\texttt{item}$ and klg stores the knowledge attached
on an edge.

The commonsense knowledge (CSK) we utilize is the social CSK provided by the
CSK graph ATOMIC-2020 Hwang et al. (2021). Hwang et al., Hwang et al. (2021)
train a generating model called COMET Bosselut et al. (2019) on the CSK graph
so that texts out of the scope of the CSK graph can be freely inferred.
ATMOIC-2020 contains nine social-interaction relations while not all relations
can be fit into our condition. According to the direction of inference and
relevance, we pick four relations: xEffect, xReact, oEffect, and oReact.
x(o)Effect reasons the effect on the speaker-self (other speakers) after an
utterance, and x(o)React describes how the speaker-self (other speakers) feels
after an utterance. We exemplify these relations in Appendix A. Given an
utterance $u_{i}$ and one of the relations, COMET takes them as the input and
generates knowledge with the beam search of 5. We denote the generation as
$\texttt{CT}(u_{i},\texttt{xE})$, where CT and xE are the abbreviations of
COMET and xEffect respectively.

As COMET does not guarantee to generate emotion-related knowledge, it is
necessary to pick up the knowledge related to the emotion of the targeted
utterance. Furthermore, knowledge passing should align with speakers’
interactions. Therefore, we propose a knowledge selecting strategy, which is
sentiment realized and fits the condition of conversations.

Algorithm 1 Construction of Knowledge adjacency matrix

Input: Dialog $[u_{1},...,u_{N}]$, emotion $[e_{1},...,e_{N}]$, speaker
$[p_{1},...,p_{N}]$, COMET CT, window $w_{k}$, func $split$($\cdot$),
$e2s$($\cdot$).
Output: Knowledge passing adjacency matrix $\mathrm{\mathbb{A}_{k}}$.

1: $\forall i,j\in[1,...,N]$, initializing $\mathrm{\mathbb{A}_{k}}$ with
$\mathrm{\mathbb{A}_{k}}[i,j].\mathrm{item}=0$,
$\mathrm{\mathbb{A}_{k}}[i,j].\mathrm{klg}=$“none”. ($i$: target, $j$: source)

2: for $j$ in $[1,N]$ do

3: $\mathrm{posK_{x}}$, $\mathrm{negK_{x}}$, $\mathrm{neuK_{x}}$ =
$split$(CT($u_{j}$, xE/xR))

4: $\mathrm{posK_{o}}$, $\mathrm{negK_{o}}$, $\mathrm{neuK_{o}}$ =
$split$(CT($u_{j}$, oE/oR))

5: for $i$ in $[j,\mathrm{min}(j+w_{k},N)]$ do

6: if $p_{i}$ == $p_{j}$ then

7: $\mathrm{K}$ = {pos: $\mathrm{posK_{x}}$, neg: $\mathrm{negK_{x}}$, neu:
$\mathrm{neuK_{x}}$}

8: else

9: $\mathrm{K}$ = {pos: $\mathrm{posK_{o}}$, neg: $\mathrm{negK_{o}}$, neu:
$\mathrm{neuK_{o}}$}

10: end if

11: if $e2s$($e_{i}$) $\neq$ neu then

12: $\mathrm{\mathbb{A}_{k}}[i,j].\mathrm{item}=1$

13: if $e2s$($e_{j}$) == neu then

14:
$\mathrm{\mathbb{A}_{k}}[i,j].\mathrm{klg}=\mathrm{K}[\texttt{neu}]+[\mathrm{sep}]+\mathrm{K}[e2s(e_{i})]$

15: else

16: $\mathrm{\mathbb{A}_{k}}[i,j].\mathrm{klg}=\mathrm{K}[e2s(e_{i})]$

17: end if

18: end if

19: end for

20: end for

21: return $\mathrm{\mathbb{A}_{k}}$

The knowledge selection realizes two factors:

Speaker realization: If $p_{i}$ of the target $u_{i}$ equals to $p_{j}$ of the
source $u_{j}$, knowledge of xEffect and xReact is chosen, otherwise knowledge
of oEffect and oReact picked.

Sentiment realization: For the generated set of knowledge, we need to first
identify its emotion, and then pick up emotion-related pieces for $u_{i}$ and
$u_{j}$. As no standard emotion lexicon resource is provided, we replace the
emotion with the sentimental polarity. We utilize sentiwordnet to compute the
sentimental polarity of every piece of knowledge. As sentiwordnet rates every
word with scores of neutral, negative, and positive, we compute the
sentimental score of a piece of knowledge by:

$r_{s}=\left\\{\begin{matrix}r_{pos}-r_{neg},&\left|r_{pos}-r_{neg}\right|>r_{neu}\\\
0,&\mathrm{else}\\\ \end{matrix}\right.$ (2)

where $r_{pos}$ is the average positive score of all words in the knowledge,
and the knowledge is regarded as positiveness when $r_{s}>0$. According to the
$r_{s}$, we split knowledge into three sets and concatenate the knowledge in a
set with the token $[\mathrm{sep}]$ (denoted as a function $split(\cdot)$). If
a set is empty, the text “none” will be assigned to it. After splitting
knowledge, we further map the emotion of $u_{i}$ and $u_{j}$ into the
sentiment (denoted as a function $e2s(\cdot)$). The emotion happiness is
mapped to positive, neutral to neutral, and other emotions in the dataset to
negative. With the mapped sentiment, we can select the related knowledge for
utterances.

The knowledge selection is processed simultaneously along the construction of
$\mathrm{\mathbb{A}_{k}}$, which is illustrated in Algorithm 1. In line 14 of
Algorithm 1, we add the neutral knowledge, because we think that it can bring
some reasoning information to help understand causes.

####  3.3.2 Knowledge Enhanced DAG Networks

As the original DAG network considers no knowledge, we modify it and propose
the Knowledge Enhanced DAG network to better process the KEC graph.

We first encode the knowledge into representations by using the shared RoBERTa
of the utterance encoder and the same operation of encoding utterance. For an
item of knowledge $\mathrm{\mathbb{A}_{k}}[i,j].\texttt{klg}$, it is encoded
as $k_{i,j}\in\mathbb{R}^{d_{u}}$.

For a targeted utterance $u_{i}$, edge weights from neighboring nodes
$\mathrm{v}_{j}$ in the $l^{th}$ layer can be computed as:

$\alpha_{i,j}=\mathrm{Softmax}_{j\in\mathcal{N}_{i}}(\mathrm{W}_{\mathrm{w}}^{l}[h_{i}^{l-1}||(h_{j}^{l}+\mathrm{W}_{\mathrm{k}}^{l}k_{i,j})]),$
(3)

where $h_{i}^{l-1}$ is the node representation of the $(l-1)^{th}$ layer,
$h_{j}^{l}$ is the neighboring node representation that has been processed in
the $l^{th}$ layer, $\mathrm{W}_{\mathrm{w}}^{l}\in\mathbb{R}^{1\times
2d_{u}}$ and $\mathrm{W}_{\mathrm{k}}^{l}\in\mathbb{R}^{d_{u}\times d_{u}}$
are trainable parameters. For $h_{i}^{0}\in\mathbb{R}^{d_{u}}$, we initialize
it by $h_{i}^{0}=\mathrm{Linear}([s_{i}||eemb_{e_{i}}])$, where $eemb_{e_{i}}$
is the trainable embedding of the emotion $e_{i}$ and is initialized by the
weight vector in the MLM head of RoBERTa. Furthermore, the contextual message
and knowledge message are aggregated from neighboring nodes according to the
edge weight $\alpha_{i,j}$:

$\displaystyle
msg_{i}=\sum_{j\in\mathcal{N}_{i}}{\alpha_{i,j}\mathrm{W}_{\mathbb{A}_{\mathrm{c}[i,j].\texttt{rel}}}^{l}h_{j}^{l}},$
(4) $\displaystyle
nlg_{i}=\sum_{j\in\mathcal{N}_{i}}{\alpha_{i,j}\mathrm{W}_{\mathrm{k}}^{l}k_{i,j}},$
(5)

where
$\mathrm{W}_{\mathbb{A}_{\mathrm{c}[i,j].\texttt{rel}}}\in(\mathrm{W}_{\mathrm{SD}}^{l},\mathrm{W}_{\mathrm{ID}}^{l})$
with the size of $\mathbb{R}^{d_{u}\times d_{u}}$ is a trainable parameter.

|  | Train | Dev | Test
---|---|---|---|---
Num. of Causal Pair | Positive | 7027 | 328 | 1767
Negative | 45392 | 2842 | 14052
Num. of Dialgoue | 834 | 47 | 225
Num. of Utterance | 8206 | 493 | 2405
Avg. Len. of Utterance | 14 | 16 | 15

Table 1: Statistics of RECCON-DD. “Positive” means the true causal pair.

To capture the nodal and contextual information, DAG networks utilize two
GRUs: nodal unit ($\mathrm{GRU}_{n}$) and contextual unit
($\mathrm{GRU}_{c}$). The information is obtained by:

$\displaystyle nod_{i}=\mathrm{GRU}_{n}(h_{i}^{l-1},msg_{i}),$ (6)
$\displaystyle cxt_{i}=\mathrm{GRU}_{c}(msg_{i},h_{i}^{l-1}).$ (7)

To further consider the effect of knowledge from the context and the utterance
itself, we employ two additional GRUs: contextual knowledge unit
($\mathrm{GRU}_{k}$) and self-loop knowledge unit ($\mathrm{GRU}_{s}$).
Therefore, knowledge-related information can be obtained by:

$\displaystyle ckg_{i}=\mathrm{GRU}_{k}(nlg_{i},h_{i}^{l-1}),$ (8)
$\displaystyle skg_{i}=\mathrm{GRU}_{s}(k_{i,i},h_{i}^{l-1}).$ (9)

Finally, the node representation of utterance $u_{i}$ in the $l^{th}$ layer is
updated by summing the four types of information, i.e.
$h_{i}^{l}=nod_{i}+cxt_{i}+ckg_{i}+skg_{i}$.

###  3.4 Cause Predictor

After the encoding of two-level encoders, the final utterance representation
is computed by $h_{i}=||_{l=0}^{L}h_{i}^{l}$ Shen et al. (2021b). Whether
$u_{j}$ is the cause of $u_{i}$ is then computed by:

$p_{i,j}=sigmoid(\mathrm{MLP}([h_{i}||h_{j}])),$ (10)

where MLP maps the concatenation of $h_{i}$ and $h_{j}$ from the dimension of
${2d_{u}}$ to 1.

We use the cross entropy loss to train our model, and the loss on a
conversation is formulated as:

$\mathcal{L}=\sum_{i\leq N}\sum_{j\leq
i}(y_{i,j}\cdot\mathrm{log}p_{i,j}+(1-y_{i,j})\cdot(1-p_{i,j})),$ (11)

where $y_{i,j}\in\\{0,1\\}$ is the cause label of the pair $(u_{i},u_{j})$.

##  4 Experimental Setup

###  4.1 Dataset

We conduct experiments on the dataset call RECCON-DD Poria et al. (2021),
which is collected from the popular dyadic dialogue dataset DailyDialog (DD)
Li et al. (2017). Preprocessing the original RECCON-DD, we only consider
causes in the dialogue history for every utterance and remove the duplicated
causal pairs. Statistics of the processed RECCON-DD are presented in Tab. 1.

| Model | Neg. F1 | Pos. F1 | macro F1
---|---|---|---|---
1 | ECPE-2D△ | 94.96 | 55.50 | 75.23
ECPE-MLL△ | 94.68 | 48.48 | 71.59
RankCP△ | 97.30 | 33.00 | 65.15
2 | KAG | 94.49(0.22) | 55.52(2.39) | 75.02(1.13)
Adapted | 95.67(0.24) | 62.47(4.72) | 79.07(2.44)
SKAIG | 95.26(0.12) | 63.15(1.00) | 79.21(0.49)
3 | Entail | 94.83(0.47) | 58.59(3.78) | 76.66(1.66)
DAG-ERC | 95.33(0.25) | 63.56(2.10) | 79.44(1.16)
4 | KEC (ours) | 95.74(0.05) | 66.76${}^{*}_{(0.33)}$ | 81.25${}^{*}_{(0.17)}$

Table 2: Results of all models on RECCON-DD. ∗ denotes that our method is
significant against the best baseline DAG-ERC (p-value$<$0.05) with the paired
T-test. △ denotes the results referred from Poria et al. [2021].

###  4.2 Compared Methods

For the reason that C2E2 is a new task proposed by Poria et al. Poria et al.
(2021), there is only one baseline Poria et al. (2021). Therefore, we
additionally compare our model with other SOTA methods from conversation and
emotion cause related tasks.

Entail is the baseline of C2E2. It utilizes RoBERTa to process the
concatenation of utterances and their emotion and history context. The
concatenation is formed by “[cls] emotion [sep] $u_{i}$ [sep] $u_{j}$ [sep]
history context [sep]”. DAG-ERC Shen et al. (2021b) is the SOTA method of
Conversational Emotion Recognition (ERC). Our method is based on it. SKAIG Li
et al. (2021) is another SOTA method of ERC which utilize CSK to help the
prediction. It passes knowledge through edges using graph transformers Shi et
al. (2021), while no knowledge selecting strategy is applied. For the fair
comparison, we utilize the knowledge x/oEffect and x/oReact for SKAIG. KAG Yan
et al. (2021) is the SOTA method of Emotion Cause Extraction (ECE) using CSK
to overcome the position bias. It connects emotional utterances with a
contextual utterance if a path is matched in ConceptNet. It uses relational
GCNs to process their proposed graph. Knowledge used in KAG only participates
in the computation of edge weights. As other methods add no position
embeddings, we do not add position embeddings in KAG. Adapted Turcan et al.
(2021) is another SOTA method of ECE. It directly concatenates CSK (xReact and
oReact) after the utterance and utilizes RoBERTa to make prediction. For
baselines from Emotion Cause Pair Extraction (ECPE), following and referring
the results from Poria et al. Poria et al. (2021), we choose ECPE-2D Ding et
al. (2020a), ECPE-MLL Ding et al. (2020b), and RankCP Wei et al. (2020). To
consider emotion information, in the re-implementation, all graph-related
baselines utilize emotion embedding.

Method | Pos. F1 |  macro F1
---|---|---
KEC |  66.76(0.33) |  81.25(0.17)
$-$ contextual knowledge unit |  64.93(0.69) |  80.25(0.37)
$-$1.83(↓) |  $-$1.00(↓)
$-$ self-loop knowledge unit |  65.00(0.92) |  80.21(0.54)
$-$1.76(↓) |  $-$1.04(↓)
$-$ neutral knowledege |  65.80(0.96) |  80.67(0.52)
$-$0.96(↓) |  $-$0.58(↓)
$-$ emotion embs |  63.19(0.25) |  79.35(0.17)
$-$3.56(↓) |  $-$1.90(↓)
$-$ emoition embs & CSK |  46.66(0.76) |  69.90(0.48)
$-$20.1(↓) |  $-$11.4(↓)

Table 3: Ablation Study

###  4.3 Implementation

All compared methods and our model use base-size pretrained models (e.g.,
RoBERTa-base) from huggingface Transformers. We train our model using the
AdamW optimizer with 40 epochs, learning rate of 3e-5, L2 regularization of
1e-4. The batch is set to 4 and we step the gradient with the accumulating
step of 2. For Conversation Encoder, we set the dimension of utterance to 300,
the hidden size of GRUs to 300, the number of layers of DAG networks to 5. For
emotion embeddings, we utilize a linear unit mapping the dimension of 768 to
200. For Cause Predictor, dimensions of MLP is set to [600, 300, 300, 1] and
the dropout rate is set to 0.1. For the window size of knowledge and that of
context, we find that our model achieves better performance when they are the
same, and we set them to 2. Following Poria et al. Poria et al. (2021), we
report the F1 scores of both negative and positive causal pairs and the macro
F1 score of them. Our experiments are conducted using 5 random seeds.

##  5 Results and Discussions

Figure 3: “Same Emotion” reports the recall of causal pairs whose causal
utterances are with the same emotion as the targeted utterance. “Different
Emotion” refers to different emotions from the targeted utterance. The number
of a type of pairs is presented below the x-axis. Green bars denote DAG-ERC
without emotion embeddings, blue bars for DAG-ERC with emotion embeddings, and
orange bars for KEC.

###  5.1 Main Results

The results of our method and baselines are listed in Tab. 2. We split the
table into 4 rows: the 1st row for models of ECPE, the 2nd row for models
using CSK, the 3rd row for models with no CSK, and the 4th row for our method.

From the 1st row, models of ECPE are not competitive to models in other rows.
We think the reason is that the setting of ECPE is not suitable for C2E2.
Compared with models in the 2nd row, our method can achieve better results.
For KAG, we think the undesired performance can be attributed to that it is
designed to resist position bias for news articles. Conversations in RECCON
tends to be short and requires more sophisticated modeling than articles. In
addition, knowledge is only used to compute edge weight in KAG, which does not
propagate through edges like ours. For “Adapted”, knowledge is directly
concatenated after the utterance, and is not explicitly passed between
utterances. For SKAIG, knowledge is also not selected, which leads to
sentiment-nonaligned knowledge passed through the graph. On the contrary, our
KEC selects proper knowledge to align with emotions, and explicitly models and
propagates the knowledge in graphs. The comparison indicates that our method
is more suitable for utilizing knowledge in C2E2.

As the baseline of RECCON, “Entail” is not comparable to the graph-based SKAIG
and DAG-ERC, which demonstrates that explicit modeling of utterance
interactions is critical for C2E2. DAG-ERC is the SOTA baseline for ERC with
no knowledge involved. By introducing CSK and the knowledge selection, our
method can achieve significant improvements against DAG-ERC. Therefore, it is
critical to utilize proper knowledge for cause detection. To verify our
knowledge enhanced DAG networks, we construct a variant that knowledge
$\mathrm{W}_{\mathrm{k}}^{l}k_{j,i}$ is directly added on $h_{j}^{l}$ in eq. 4
and no knowledge-related GRUs are involved. The macro F1 and Pos. F1 achieved
by the variant are 80.38 and 65.24. This demonstrates that, instead of being
simply added in contextual messages, knowledge prefers to be exclusively
modelled by gated units in our knowledge enhanced DAG networks.

Figure 4: The effect of the window size. The y-axis denotes the macro F1 on
the development set. To better fit the trend, we draw the trend line of fifth-
order polynomial as the red dotted line.

###  5.2 Ablation Study

To study the effect of different modules equipped in our model, we remove
contextual knowledge unit $\mathrm{GRU}_{k}$, self-loop knowledge unit
$\mathrm{GRU}_{s}$, neutral knowledge added when source utterances are neutral
in line 14 of Algorithm 1, and emotion embedding respectively. The results are
in Tab. 3.

The decreased performance by removing $\mathrm{GRU}_{k}$ shows the necessity
of controlling knowledge from the context. As for removing $\mathrm{GRU}_{s}$
affects the performance, it indicates that self-knowledge enriching for
targeted utterances is critical because 37$\%$ causal utterances are the
targeted utterances themselves in RECCON. By taking off neutral knowledge, the
performance shows decline. This manifests that neutral knowledge can help the
model to better detect causes as it can provide some reasoning information.

When removing emotion embedding, the performance is still competitive to DAG-
ERC who is equipped with emotion embedding. This can be attributed to that our
knowledge selecting strategy can pick up sentiment related knowledge for the
model. As emotion embedding and knowledge are all removed, the model becomes
DAG-ERC without emotion information. The performance drops dramatically, which
verify the claim that emotion information are important to C2E2. As for
“Entail” without emotion, the performance is presented in Appendix C. In the
next subsection, we will further show what causes are preferred by the emotion
information.

###  5.3 Performance on Cause Pairs

There is a question about whether our model can promote the detection of
causal utterances with emotions different from that of the targeted utterance.
We therefore present the recalls of causal pairs whose utterances are with the
same emotion (SE pairs) and different emotions (DE pairs) in Fig. 3 for every
targeted emotion. From the figure, we can see that DAG-ERC with emotion
embedding can dramatically boost the detection of SE pairs while limited
improvement is achieved on DE pairs. Specifically, emotion information
increases the recall on all SE pairs by about 23$\%$ while only about 3$\%$ on
DE pairs. By introducing CSK, our model, compared with DAG-ERC with emotion
embedding, further increases the recall by 3$\%$ on SE pairs and by 12$\%$ on
DE pairs (especially boosting the detection of neutral causal utterances by
12$\%$). This manifests that our model can promote the detection of DE pairs,
especially on pairs whose targeted utterances are with negative emotions like
Anger, Sadness in Fig. 3.

###  5.4 Effect of Window Size

To show the effect of the window size, we illustrate the trend of macro F1
scores with the window size increasing in Fig. 4. The trend line shows that
our model achieves the best performance with a small window. As the window
size keeps increasing, the performance will drop and fluctuate. The reason for
this can be attributed to that over 80$\%$ causal utterances are located near
the targeted utterance in no more than 2 steps. It can be future work to solve
the problem that models bias to neighboring utterances in C2E2.

##  6 Conclusion

In this paper, we study Conversational Causal Emotion Entailment. We find that
adding emotion information into graphs cannot effectively detect causal
utterances with different emotions from the targeted utterance due to models’
limited causal reasoning ability. To alleviate this, we introduce social
commonsense knowledge and propose a graph-based structure to properly utilize
knowledge. Specifically, we propose a Knowledge Enhanced Conversation graph
(KEC) with a designed knowledge selecting strategy. To process KEC, we
construct Knowledge Enhanced Directed Acyclic Graph networks. Our method
outperforms baselines and detects more causes with different emotions from the
targeted utterance.

## Acknowledgments

This work was supported by National Natural Science Foundation of China (No.
61976207, No. 61906187).

## References

  * Bosselut et al. [2019] Antoine Bosselut, Hannah Rashkin, Maarten Sap, Chaitanya Malaviya, Asli Celikyilmaz, and Yejin Choi.  COMET: commonsense transformers for automatic knowledge graph construction.  In ACL, 2019.
  * Ding et al. [2019] Zixiang Ding, Huihui He, Mengran Zhang, and Rui Xia.  From independent prediction to reordered prediction: Integrating relative position and global label information to emotion cause identification.  In AAAI, pages 6343–6350, 2019.
  * Ding et al. [2020a] Zixiang Ding, Rui Xia, and Jianfei Yu.  ECPE-2D: emotion-cause pair extraction based on joint two-dimensional representation, interaction and prediction.  In ACL, pages 3161–3170, 2020.
  * Ding et al. [2020b] Zixiang Ding, Rui Xia, and Jianfei Yu.  End-to-end emotion-cause pair extraction based on sliding window multi-label learning.  In EMNLP, 2020.
  * Ghosal et al. [2019] Deepanway Ghosal, Navonil Majumder, Soujanya Poria, Niyati Chhaya, and Alexander F. Gelbukh.  Dialoguegcn: A graph convolutional neural network for emotion recognition in conversation.  In EMNLP-IJCNLP, pages 154–164, 2019.
  * Ghosal et al. [2020] Deepanway Ghosal, Navonil Majumder, Alexander F. Gelbukh, Rada Mihalcea, and Soujanya Poria.  COSMIC: commonsense knowledge for emotion identification in conversations.  In Findings of EMNLP, pages 2470–2481, 2020.
  * Ghosal et al. [2021] Deepanway Ghosal, Navonil Majumder, Rada Mihalcea, and Soujanya Poria.  Stack: Sentence ordering with temporal commonsense knowledge.  In EMNLP, 2021.
  * Gui et al. [2016] Lin Gui, Ruifeng Xu, Qin Lu, Dongyin Wu, and Yu Zhou.  Emotion cause extraction, A challenging task with corpus construction.  In SMP, 2016.
  * Hwang et al. [2021] Jena D. Hwang, Chandra Bhagavatula, Ronan Le Bras, Jeff Da, Keisuke Sakaguchi, Antoine Bosselut, and Yejin Choi.  (comet-) atomic 2020: On symbolic and neural commonsense knowledge graphs.  In AAAI, pages 6384–6392, 2021.
  * Kim et al. [2021] Hyunwoo Kim, Byeongchang Kim, and Gunhee Kim.  Perspective-taking and pragmatics for generating empathetic responses focused on emotion causes.  In EMNLP, pages 2227–2240, 2021.
  * Li et al. [2017] Yanran Li, Hui Su, Xiaoyu Shen, Wenjie Li, Ziqiang Cao, and Shuzi Niu.  Dailydialog: A manually labelled multi-turn dialogue dataset.  In IJCNLP, 2017.
  * Li et al. [2021] Jiangnan Li, Zheng Lin, Peng Fu, and Weiping Wang.  Past, present, and future: Conversational emotion recognition through structural modeling of psychological knowledge.  In Findings of EMNLP, 2021.
  * Liang et al. [2021] Yunlong Liang, Fandong Meng, Ying Zhang, Yufeng Chen, Jinan Xu, and Jie Zhou.  Infusing multi-source knowledge with heterogeneous graph neural network for emotional conversation generation.  In AAAI, 2021.
  * Liu et al. [2019] Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov.  Roberta: A robustly optimized BERT pretraining approach.  CoRR, abs/1907.11692, 2019.
  * Majumder et al. [2019] Navonil Majumder, Soujanya Poria, Devamanyu Hazarika, Rada Mihalcea, Alexander F. Gelbukh, and Erik Cambria.  Dialoguernn: An attentive RNN for emotion detection in conversations.  In AAAI, 2019.
  * Poria et al. [2021] Soujanya Poria, Navonil Majumder, Devamanyu Hazarika, Deepanway Ghosal, Rishabh Bhardwaj, Samson Yu Bai Jian, Pengfei Hong, Romila Ghosh, Abhinaba Roy, Niyati Chhaya, Alexander F. Gelbukh, and Rada Mihalcea.  Recognizing emotion cause in conversations.  Cogn. Comput., 13(5):1317–1332, 2021.
  * Shen et al. [2021a] Weizhou Shen, Junqing Chen, Xiaojun Quan, and Zhixian Xie.  Dialogxl: All-in-one xlnet for multi-party conversation emotion recognition.  In AAAI, 2021.
  * Shen et al. [2021b] Weizhou Shen, Siyue Wu, Yunyi Yang, and Xiaojun Quan.  Directed acyclic graph network for conversational emotion recognition.  In ACL/IJCNLP, pages 1551–1560, 2021.
  * Shi et al. [2021] Yunsheng Shi, Zhengjie Huang, Shikun Feng, Hui Zhong, Wenjing Wang, and Yu Sun.  Masked label prediction: Unified message passing model for semi-supervised classification.  In IJCAI, 2021.
  * Turcan et al. [2021] Elsbeth Turcan, Shuai Wang, Rishita Anubhai, Kasturi Bhattacharjee, Yaser Al-Onaizan, and Smaranda Muresan.  Multi-task learning and adapted knowledge models for emotion-cause extraction.  In ACL/IJCNLP, pages 3975–3989, 2021.
  * Wei et al. [2020] Penghui Wei, Jiahao Zhao, and Wenji Mao.  Effective inter-clause modeling for end-to-end emotion-cause pair extraction.  In ACL, 2020.
  * Xia and Ding [2019] Rui Xia and Zixiang Ding.  Emotion-cause pair extraction: A new task to emotion analysis in texts.  In ACL, pages 1003–1012, 2019.
  * Xia et al. [2019] Rui Xia, Mengran Zhang, and Zixiang Ding.  RTHN: A rnn-transformer hierarchical network for emotion cause extraction.  In IJCAI, 2019.
  * Yan et al. [2021] Hanqi Yan, Lin Gui, Gabriele Pergola, and Yulan He.  Position bias mitigation: A knowledge-aware graph model for emotion cause extraction.  In ACL/IJCNLP, pages 3364–3375, 2021.
  * Zhong et al. [2019] Peixiang Zhong, Di Wang, and Chunyan Miao.  Knowledge-enriched transformer for emotion detection in textual conversations.  In EMNLP-IJCNLP, pages 165–176, 2019.


# Mid-Infrared Period–Luminosity Relations for Miras in the Large Magellanic
Cloud

Patryk Iwanek Astronomical Observatory, University of Warsaw, Al. Ujazdowskie
4, 00-478 Warsaw, Poland  Igor Soszyński Astronomical Observatory, University
of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland  Szymon Kozłowski
Astronomical Observatory, University of Warsaw, Al. Ujazdowskie 4, 00-478
Warsaw, Poland

(Accepted June 25, 2021)

###### Abstract

We present the mid-infrared (mid-IR) period–luminosity relations (PLRs) using
over 1000 Mira variables in the Large Magellanic Cloud (LMC), for the four
Wide-Field Infrared Survey Explorer (WISE) and the four Spitzer bands. These
PLRs cover a mid-IR wavelength range from $3.4$ $\mu$m to $22$ $\mu$m and are
presented separately for the oxygen-rich (O-rich) and carbon-rich (C-rich)
Miras. These relations can be used to measure distances to individual O-rich
and/or C-rich Mira stars with the accuracy of $5\%$ and $12\%$, respectively.
They are the most accurate Mira PLRs in the mid-IR to date.

stars: AGB and post-AGB – stars: carbon – stars: variables: general –
Magellanic Clouds

##  1 Introduction

Mira variables are fundamental-mode pulsators with periods ranging from
$\sim$100 days to over $1000$ days. They belong to the class of late-type,
low- or intermediate-mass pulsating red giants also known as Long Period
Variables (LPVs). Miras are Asymptotic Giant Branch (AGB) tracers of old- and
intermediate-age populations (Iben & Renzini, 1983; Whitelock, 2012). Like
other AGB stars, they can be divided into oxygen-rich (O-rich) and carbon-rich
(C-rich) giants, depending on their surface composition (Riebel et al., 2010).
Miras, thanks to their large brightness variations, can be easily found in the
Milky Way and other Local Group galaxies, and they are an extensively studied
subclass of LPVs.

Since 1912, radially pulsating stars have been widely used as Galactic and
extragalactic distance indicators thanks to the discovery of period–luminosity
relation (PLR; Leavitt & Pickering 1912) for classical Cepheids. The history
of research on PLRs for Miras dates back to the 20s of the last century and
began when Gerasimovič (1928) noticed that Miras with longer periods are on
average fainter at optical wavelengths. This result has been confirmed by e.g.
Wilson & Merrill (1942).

Glass & Lloyd Evans (1981) were the first to report the existence of PLR for
Miras at the near-infrared (NIR) wavelengths. Their research was based on 11
Miras from the Large Magellanic Cloud (LMC), and it appeared that PLR for
Miras in NIR was much tighter than in the optical passbands. This NIR PLR was
redefined by Feast et al. (1989) and Hughes & Wood (1990).

The availability of long-term photometric data for millions of stars from
large-scale sky surveys such as MACHO (Alcock et al., 2000), EROS (Expérience
de Recherche d’Objets Sombres; Ansari, 1996), OGLE (Optical Gravitational
Lensing Experiment; Udalski et al., 2015) allowed an in-depth study of the
PLRs of LPVs. Based on the MACHO data, Wood et al. (1999) and then Wood (2000)
showed five distinct, parallel sequences in the period–luminosity (PL) plane.
It turned out that one of these sequences was composed of Miras (sequence C in
Fig. 1 in Wood 2000).

The research on pulsating red giants conducted by the OGLE team members has
significantly expanded our knowledge about PLRs for such stars. For instance,
Soszyński et al. (2004) showed that OGLE Small Amplitude Red Giants (OSARGs)
form a completely different set of PL sequences than Miras and Semiregular
Variables (SRVs). In turn, Soszyński et al. (2005) increased the complexity of
the PL distribution of LPVs and proposed a photometric method to distinguish
between oxygen-rich (O-rich) and carbon-rich (C-rich) AGB stellar populations.

Many modern studies on the Mira PLRs (e.g. Whitelock et al. 2006, 2008; Riebel
et al. 2010; Ita & Matsunaga 2011; Riebel et al. 2015; Whitelock et al. 2017;
Yuan et al. 2017a, b, 2018; Bhardwaj et al. 2019) showed that these stars can
be used as an excellent extension of the Cosmic Distance Ladder. Due to their
well-defined PLRs and high luminosity in the infrared passbands, which are
less affected by interstellar extinction, Miras can be used as an attractive
distance indicator (e.g. Matsunaga et al. 2009; Qin et al. 2018; Huang et al.
2018; Molina et al. 2019; Huang et al. 2020; Urago et al. 2020).

Recently, Iwanek et al. (2021) analyzed Mira variability in a wide range of
wavelengths, covering $0.1$ to $40$ microns. The authors used densely-covered,
accurate OGLE $I-$band light curves to create templates and fit them to NIR
and mid-infrared (mid-IR) data from multiple sky surveys, extending the study
to longer wavelengths (up to $40$ microns) by fitting spectral energy
distributions (SEDs). As a result, they determined the variability amplitude
ratios and phase lags for a range of wavelengths. Additionally, the authors
presented synthetic PLRs in near- and mid-IR for the existing and future sky
surveys.

The aim of this paper is to derive the mid-IR PLRs for LMC Miras discovered in
the OGLE-III data (Soszyński et al., 2009), using the Wide Field Infrared
Survey Explorer (WISE, Wright et al., 2010) and Spitzer (Werner et al., 2004)
data, the most accurate distance to the LMC (Pietrzyński et al., 2019), and
the reddening map of the LMC (Skowron et al., 2021).

Despite many previous works on the Mira PLRs, this subject has not been
comprehensively studied in the mid-IR. Using the reddening map of the LMC and
the distance to the LMC accurate to one percent, we provide the most accurate
Mira PLRs to date. Additionally, for the first time, we present PLRs
separately for the O-rich and C-rich Miras, in the WISE W1, W2, W3, and W4
bands.

##  2 Data

###  2.1 Sample of Miras

OGLE is one of the largest and the longest lasting time-domain variability sky
surveys worldwide. Since 1992, the Galactic bulge (BLG), and later also the
Magellanic Clouds (MCs) and the Galactic disk (GD) have been regularly
monitored to search for variability. These observations have enabled many
scientific discoveries, catalogs of variable stars, and two-band, multi-
decades-long time-series photometry.

One of the products of the third phase of the OGLE project (OGLE-III, Udalski
2003) is the catalog of LPVs in the LMC (Soszyński et al., 2009). The authors
discovered almost $100\;000$ LPVs, including OSARGs, SRVs and Mira-type stars.
The two latter subclasses were separated using the $I$-band pulsation
amplitude (for Miras $\Delta I$ $>$ $0.8$ mag). Using color-color ($(V-I)$ vs.
$(J-K)$) and Wesenheit ($W_{I}$ vs. $W_{JK}$) diagrams, Soszyński et al.
(2009) divided Miras into O-rich and C-rich stars. This division was evaluated
and confirmed by Ita & Matsunaga (2011). In the original catalog, Soszyński et
al. (2009) published two-band ($I$ and $V$ in the Johnson-Cousins photometric
system), 13-year-long light curves covering the time span of 1996–2009. The
full catalog can be found in the on-line databases through the OGLE
webpage111http://ogle.astrouw.edu.pl

###  2.2 Optical data

In this paper, we used 1663 published LMC Miras (Soszyński et al., 2009) with
provided coordinates, and surface chemistry classifications (1194 C-rich and
469 O-rich). The OGLE-III phase ended in 2009, and at the beginning of 2010
the OGLE project restarted its observations as the fourth phase (OGLE-IV). A
detailed description of the OGLE-IV phase, instrumentation, calibration, data
reduction, observations’ cadence, or the sky coverage can be found in Udalski
et al. (2015).

We supplemented publicly available Mira light curves with the non-public OGLE-
IV data collected over an additional 10 years since 2010. We revised their
pulsation periods using two-decades-long light curves. To find periods, we
used TATRY code based on the multiharmonic analysis of the variance algorithm
(ANOVA, Schwarzenberg-Czerny, 1996). We visually inspected each light curve,
phase-folded it with its new pulsation period, and manually removed obvious
outlying points. The updated OGLE light curves have on average $\sim$1300 and
$\sim$100 data points in the I-band and V-band, respectively, obtained between
December 29th, 1996 and March 15th, 2020.

###  2.3 Mid-IR data

The role of infrared observations in studying the PLRs and measuring cosmic
distances is invaluable because the influence of interstellar dust on the
stellar light decreases with the increasing wavelength. Moreover, Mira
variables are known to undergo the mass-loss phenomenon (e.g. Perrin et al.
2020), may be the cause of the formation of circumstellar dust shells around
these stars. The impact of circumstellar shells on the observed stellar
radiation, as for the interstellar dust, is smaller at longer wavelengths.

We cross-matched our catalog of Mira stars with two space mission databases
containing mid-IR observations: WISE (Wright et al., 2010) and Spitzer (Werner
et al., 2004). WISE is a 40-cm diameter infrared space telescope that has
monitored the entire sky in four bands $W1$ ($3.4$ $\mu$m), $W2$ ($4.6$
$\mu$m), $W3$ ($12$ $\mu$m) and $W4$ ($22$ $\mu$m). In the standard mode, the
WISE telescope has observed each sky location every six months collecting
several independent exposures during one flyby. The situation was slightly
different in the case of the LMC area, which was observed much more frequently
due to the polar trajectory of the telescope. As a result, the light curves of
stars located in the LMC were covered with much denser sampling. The main
mission of the WISE telescope ended in 2010, however, observations in $W1$ and
$W2$ have been continued since 2011 as part of the Near Earth Object WISE
Reactivation Mission (NEOWISE-R Mainzer et al., 2011, 2014). For this
analysis, we retrieved the data from AllWISE Multiepoch Photometry
Table222https://wise2.ipac.caltech.edu/docs/release/allwise/, and NEOWISE-R
Single Exposure (L1b) Source Table using NASA/IPAC Infrared Science
Archive333https://irsa.ipac.caltech.edu/applications/Gator/.

The Spitzer Space Telescope is an 85-cm diameter telescope with three infrared
instruments onboard: Infrared Array Camera (IRAC), Infrared Spectrograph
(IRS), and Multiband Imaging Photometry for Spitzer (MIPS). The LMC was
observed by Spitzer in the Surveying the Agents of Galaxy’s Evolution (SAGE,
Meixner et al. 2006) survey using IRAC in four bands: $3.6$ $\mu$m, $4.5$
$\mu$m, $5.8$ $\mu$m, and $8.0$ $\mu$m, also called $[3.6]$, $[4.5]$, $[5.8]$,
and $[8.0]$, or $I1$, $I2$, $I3$, and $I4$, respectively. We downloaded all
available measurements from the SAGE IRAC Epoch 1 and Epoch 2 Catalog using
NASA/IPAC Infrared Science Archive. Light curves from these databases contain
mostly two epochs.

The above-mentioned databases were searched for counterparts within 1 arcsec
around Miras’ coordinates. For further analysis, we used Spitzer light curves
containing two epochs in each IRAC band only. This left us with $1470$ light
curves (out of 1663) collected from October to November 2005. The WISE data
had to be cleared of significant outliers. After two-step cleaning procedure,
we removed light curves with less than $100$ data points in the $W1$\- and
$W2$\- and less than $3$ epochs in the $W3$\- and $W4$-bands from further
analyses. Finally, this left us with 1311 Miras observed by the WISE
telescope. The median number of data points per light curve were 645, 703, 63,
22 in $W1$-, $W2$-, $W3$-, and $W4$-bands, respectively. The WISE data were
collected between 8th February 2010 and 17th June 2010 (AllWISE data), and
between 13th December 2010 and 1st December 2019 (NEOWISE-R data). A detailed
description of data extraction and cleaning procedure can be found in Iwanek
et al. (2021). In Figure 1, we present transmission curves of the WISE and
Spitzer bands.

The analysis of spectral energy distributions (SEDs) of the LMC Miras
performed by Iwanek et al. (2021) showed, that photometric measurements in the
$W4-$band ($22$ $\mu$m) are typically significantly biased. While the OGLE,
VMC (The VISTA survey of the Magellanic Clouds system), Spitzer, WISE ($<20$
$\mu$m) and MIPS measurements (one of the Spitzer instrument that observed the
sky at $24$ $\mu$m), were well decribed by either a single (in the case of
O-rich Miras) or a double (in the case of C-rich Miras) Planck function, the
$W4$ measurements clearly deviated from the fit. We conducted the analysis of
the PLR in the $W4-$band for completeness purposes only, but they should not
be used in future analyses.

Figure 1: Transmission curves for WISE and Spitzer bands.

##  3 Methods

###  3.1 Corrections for the interstellar extinction

Interstellar extinction has a strong influence on stellar light and should be
taken into account when measuring accurate PLRs. Although this effect is
smaller at longer wavelengths and the amount of dust toward the LMC is much
smaller than, for example, toward the Galactic bulge, it should not be
neglected.

Recently, Skowron et al. (2021) published the most accurate optical reddening
maps of the Magellanic Clouds using Red Clump (RC) stars. The reddening
$E(V-I)$ and the upper $\sigma_{\mathrm{+,E(V-I)}}$ and lower
$\sigma_{\mathrm{-,E(V-I)}}$ uncertainties can be obtained for a given Right
Ascension (RA) and Declination (DEC) via an on-line form on the OGLE
webpage444http://ogle.astrouw.edu.pl/cgi-ogle/get_ms_ext.py. Following the
authors’ remarks, the extinction $A_{I}$ could be calculated as:

$A_{\mathrm{I}}\simeq 1.5\times E(V-I),$ (1)

where the coefficient can vary between 1.1 and 1.7, depending on the dust
characteristics. Therefore, we used the coefficient equal to 1.5 with the
uncertainty equal to 0.2, what is consistent with the coefficient $1.505$ for
$A_{I}$ derived for $R_{V}=A_{V}/E(B-V)=3.1$ by Schlafly & Finkbeiner (2011).
Knowing that:

$E(V-I)=A_{\mathrm{V}}-A_{\mathrm{I}},$ (2)

the extinction in V-band is represented by the relation:

$A_{\mathrm{V}}=(2.5\pm 0.2)\times E(V-I),$ (3)

with the uncertainty:

$\sigma_{\mathrm{A_{V}}}=2.5\times
E(V-I)\times\sqrt{\left(\frac{0.2}{2.5}\right)^{2}+\left(\frac{\sigma_{\mathrm{+,E(V-I)}}+\sigma_{\mathrm{-,E(V-I)}}}{2\times
E(V-I)}\right)^{2}}.$ (4)

The V-band extinction could be transformed to extinctions in the WISE and
Spitzer bands using extinction laws derived by Chen et al. (2018) and Wang &
Chen (2019). We summarize these transformations in Table 1 for WISE W1, W2,
W3, and Spitzer $[3.6]$ $\mu$m, $[4.5]$ $\mu$m, $[5.8]$ $\mu$m, $[8.0]$ $\mu$m
bands. Considering the uncertainty $\sigma_{\mathrm{A_{V}}}$ (see Equation
(4)) and the transformation uncertainties
$\sigma_{\mathrm{A_{\lambda}}/A_{\mathrm{V}}}$ presented in Table
LABEL:tab:extinction_coefficients, the extinction uncertainties
$\sigma_{\mathrm{A_{\lambda}}}$ in a given WISE or Spitzer band can be
calculated as:

$\sigma_{\mathrm{A_{\lambda}}}=A_{\mathrm{\lambda}}\times
A_{\mathrm{V}}\times\sqrt{\left(\frac{\sigma_{\mathrm{A_{\lambda}}/A_{\mathrm{V}}}}{A_{\mathrm{\lambda}}/A_{\mathrm{V}}}\right)^{2}+\left(\frac{\sigma_{\mathrm{A_{V}}}}{A_{\mathrm{V}}}\right)^{2}}.$
(5)

As there is a lack of information in the literature on the extinction
transformation from the $V-$band to $W4-$band, we simply assumed that
$A_{\mathrm{W4}}$ is equal $0$. Iwanek et al. (2021) showed that the median
$A_{\mathrm{[8.0]}}$ is equal to $0.007$, and is one order of magnitude
smaller than extinction in $J-$band ($1.25$ $\mu$m). Knowing that the
influence of interstellar dust on stellar light decreases with increasing
wavelength (see e.g. Cardelli et al. 1989), the true extinction in $W4-$band
toward the LMC must be close to $0$.

Moreover, Iwanek et al. (2021) compared the extinction obtained with the
method described above with the classical reddening law derived by Cardelli et
al. (1989). The authors concluded that the difference between these two
methods is less than $10\%$, while extinction in the mid-IR is, in general,
comparable to a typical photometric uncertainty for a survey.

The mean value of extinction uncertainties in WISE (with an exception of the
$W4-$band) and Spitzer bands are presented in Table
LABEL:tab:extinction_calculated. Throughout the paper, all stellar brightness
and colors are extinction corrected.

Table 1: Extinction transformation coefficients between the OGLE V-band and WISE or Spitzer bands with the uncertainties. The extinction $A_{\mathrm{\lambda}}$ can be simply obtained by multiplying $A_{\mathrm{V}}$ (see Equation (3)) and the coefficient $A_{\mathrm{\lambda}}/A_{\mathrm{V}}$. $\lambda$ | $A_{\mathrm{\lambda}}/A_{\mathrm{V}}$ | $\sigma_{\mathrm{A_{\lambda}}/A_{\mathrm{V}}}$
---|---|---
WISE W1 | 0.039 | 0.004
WISE W2 | 0.026 | 0.003
WISE W3 | 0.040 | 0.009
Spitzer $[3.6]$ | 0.037 | 0.003
Spitzer $[4.5]$ | 0.026 | 0.003
Spitzer $[5.8]$ | 0.019 | 0.003
Spitzer $[8.0]$ | 0.025 | 0.003
Table 2: Mean values of calculated extinction with the uncertainties in WISE (with the exception of $W4-$band) and Spitzer bands. $\lambda$ | $A_{\lambda,\mathrm{mean}}$ | $\sigma_{\mathrm{A_{\lambda,\mathrm{mean}}}}$
---|---|---
WISE W1 | $0.012$ | $0.003$
WISE W2 | $0.008$ | $0.002$
WISE W3 | $0.012$ | $0.003$
Spitzer $[3.6]$ | $0.011$ | $0.002$
Spitzer $[4.5]$ | $0.008$ | $0.002$
Spitzer $[5.8]$ | $0.006$ | $0.001$
Spitzer $[8.0]$ | $0.008$ | $0.002$
Table 3: Basic parameters for 1663 LMC Miras from the Soszyński et al. (2009) catalog. We provide the coordinates, surface chemistry classification, pulsation periods, and mean magnitudes with uncertainties for all WISE and Spitzer bands. Mean magnitudes presented in this table are corrected for the interstellar extinction. Number | RA $(^{h}:^{m}:^{s})$ | Decl. $(^{\circ}:^{m}:^{s})$ | type |  $P$ (d) |  $m_{W1}$ (mag) |  $\sigma_{m,W1}$ (mag) | … |  $m_{[8.0]}$ (mag) |  $\sigma_{m,[8.0]}$ (mag)
---|---|---|---|---|---|---|---|---|---
00055 | 04:29:49.87 | -70:19:00.7 | C | 289.78 | 9.315 | 0.018 | … | -9.999 | -9.999
00082 | 04:30:44.96 | -69:50:41.0 | O | 162.33 | -9.999 | -9.999 | … | 10.959 | 0.061
00094 | 04:30:54.95 | -69:28:35.5 | C | 335.87 | 10.048 | 0.017 | … | 9.585 | 0.019
00096 | 04:30:59.53 | -69:57:16.1 | C | 385.50 | 9.495 | 0.022 | … | 7.884 | 0.035
00098 | 04:31:03.28 | -69:34:15.3 | C | 323.23 | 10.064 | 0.009 | … | 8.606 | 0.023
00115 | 04:31:27.40 | -70:40:41.4 | C | 176.04 | -9.999 | -9.999 | … | -9.999 | -9.999
00144 | 04:31:54.49 | -68:42:26.3 | C | 369.67 | 9.879 | 0.014 | … | 7.974 | 0.144
00225 | 04:33:17.24 | -68:09:28.8 | C | 487.44 | 9.539 | 0.063 | … | 6.908 | 0.061
00265 | 04:33:43.68 | -70:09:50.5 | C | 444.04 | 9.885 | 0.041 | … | 7.198 | 0.044
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | $\ddots$ | ⋮ | ⋮
91928 | 06:16:49.54 | -70:43:03.7 | C | 370.72 | 9.892 | 0.018 | … | -9.999 | -9.999

###  3.2 Template light curves and mean magnitudes

To measure mean magnitudes from sparsely sampled variable mid-IR light curves,
we first modeled the high-cadence OGLE $I$-band data with the semi-parametric
Gaussian process regression (GPR) model (He et al., 2016). In short, the model
consists of a number of independent parts that include the mean magnitude, a
low-frequency trend across cycles, a periodic term that models the main
variability, and finally the stochastic term that absorbs any short-term
deviations. We then fitted the GPR templates by modifying their amplitude,
shifting in time and magnitude to match the mid-IR light curves. Having the
best parameters for the variability amplitude ratio, phase-shift, and
magnitude shift, we transformed the light curve template to the mid-IR
datasets of interest. A detailed description of making templates and fitting
them to mid-IR data can be found in Iwanek et al. (2021).

We used the transformed templates to measure the mean magnitude in the WISE
and Spitzer bands. Each template light curve was transformed to the flux
scale, fitted with a third-order truncated Fourier series, integrated to
determine the mean brightness, and finally transformed back to the magnitude
scale. The uncertainties of the mean magnitude were determined from the
$\chi^{2}$ surface. In Table 3, we present basic parameters for each Mira star
from our sample (ID, coordinates, type, pulsation period) and measured mean
brightness with uncertainty in each WISE and Spitzer band.

##  4 mid-IR Period–Luminosity and Period–Luminosity–Color Relations for Miras

PLRs for Miras are tighter and better defined at infrared wavelengths than in
optical bands, as Mira PLRs are generally flat in optical (see e.g., Iwanek et
al. 2021 for an explanation of this effect), but also blurred by stronger
effects of the interstellar and circumstellar extinction, and hence harder to
work with. The use of the Wesenheit index improves the PLRs in optical bands,
however, they are still quite broad. Adding the color-term to the relations in
optical bands (thus fitting period-luminosity-color relations; PLC) makes
these relations tighter (see e.g. Bhardwaj et al., 2019). In the NIR (e.g.
$JHK$-bands), the PLRs for O-rich Miras are tight, while PLRs for C-rich Miras
have a large scatter. For Miras of both types, PLRs are present in the mid-IR
bands, but the scatter of C-rich Mira PRLs is significantly larger than that
of O-rich Miras (for comparison see e.g. Figure 2 in Ita & Matsunaga 2011 or
Figure 8 in Yuan et al. 2017b).

In recent years, various research groups have presented different approaches
to PLRs for Miras. For instance, Ita & Matsunaga (2011) compared PLRs for
O-rich and C-rich Miras and concluded that the short-period Miras (with
$\log{P}\lesssim 2.3$) of both types, indeed obey the same relation. On the
other hand, this sample represents only a small fraction of the entire range
of pulsation periods occupied by Miras.

Due to the large scatter of C-rich Mira PLRs, Ita & Matsunaga (2011) decided
to analyze the O-rich Miras only. They introduced PLRs in the Spitzer mid-IR
bands as a combination of two linear fits, with a kink at $\log{P}=2.6$.
Bhardwaj et al. (2019) found that this kink shifts to longer periods with the
increasing wavelength. The authors showed that in the optical bands this kink
is at $\log{P}=2.48$, while in the NIR at $\log{P}=2.54$. These results are
consistent with Ita & Matsunaga (2011) confirming a kink at even longer
periods in the mid-IR bands.

Figure 2: Color-magnitude diagrams (CMDs; upper row) and period–color diagrams
(lower panel) for the LMC Miras. CMDs are colored by the pulsation period,
while period–color diagrams present division into C-rich (red) and O-rich
(blue) Miras. Left plots present NIR color $(J-K)$, while middle and right
plots present mid-IR colors, $(W1-W2)$ and $([3.6]-[5.8])$, respectively.

Another approach to O-rich Miras assumed construction of PLRs by using a
quadratic fit. For example, this approach was used by Yuan et al. (2017b) to
determine PLRs for less than 180 O-rich Miras from the central region of the
LMC, using time-series observations in $JHK_{s}$ bands. The same method was
used to construct PLRs for Miras discovered in the M33 galaxy (Yuan et al.,
2017a, 2018), or to construct LMC Miras PLRs based on the Gaia DR2 data (Gaia
Collaboration et al., 2018; Bhardwaj et al., 2019).

Other authors proposed a linear fit to both O- and C-rich Miras separately
(e.g. Riebel et al. 2010, 2015), or even a single linear fit to both groups of
Miras in the LMC (e.g. Urago et al. 2020).

###  4.1 Fitting procedure

We used the observed mean magnitudes and pulsating periods presented in Table
3 to create PLRs. Our goal was to fit C-rich and O-rich Miras separately,
therefore we adopted the same division as provided by Soszyński et al. (2009,
which we also present in Table 3). The PLRs for the C-rich Miras were fitted
with the linear function:

$m_{\mathrm{\lambda,fit}}=a_{\mathrm{0,\lambda}}+a_{\mathrm{1,\lambda}}\times(\log{P}-2.3),$
(6)

where $\lambda$ is any of the WISE and Spitzer bands.

Since the PLRs for O-rich Miras appear to be non-linear, we made both linear
and quadratic fits. The quadratic fit was made for the entire range of
pulsation periods and had a form:

$m_{\mathrm{\lambda,fit}}=a_{{0,\mathrm{\lambda}}}+a_{{\mathrm{1,\lambda}}}\times(\log{P}-2.3)+a_{{\mathrm{2,\lambda}}}\times(\log{P}-2.3)^{2}.$
(7)

The O-rich Miras with pulsation periods shorter than $\sim 400$ days (i.e.
$\log P=2.6$) follow the linear PLRs (e.g. Ita & Matsunaga 2011). Therefore,
we also performed fits in the form as presented by Equation (6) for Miras with
such pulsation periods.

We fitted models to all O-rich and C-rich Miras using the weighted least
squares method, using sigma-clipping procedures, with weights that are
represented by the covariance matrix $\mathcal{C}_{kl}$, and its inverse
$\mathcal{B}_{kl}\equiv\mathcal{C}_{kl}^{-1}$. Following Gould (2003)
notation:

$\mathcal{C}_{kl}=\delta_{kl}\times(\sigma^{2}_{\mathrm{m,\lambda,k}}+\sigma_{\lambda}^{2}),$
(8)

$\chi^{2}=\sum_{k=1}^{N}\sum_{l=1}^{N}(m_{\mathrm{\lambda},k}-m_{\mathrm{\lambda,fit},k})\mathcal{B}_{kl}(m_{\mathrm{\lambda},l}-m_{\mathrm{\lambda,fit},l}),$
(9)

where $\delta_{kl}$ is the Kronecker delta ($\delta_{kl}=1$ when $k=l$, or
$\delta_{kl}=0$, when $k\neq l$). The diagonal components of the
$\mathcal{C}_{kl}$ matrix (Equation (8)) contain the sum of the squares of the
uncertainty of the mean magnitude $\sigma_{\mathrm{m,\lambda}}$ and the
intrinsic PLR scatter $\sigma_{\lambda}$, which is updated iteratively during
a sigma-clipping procedure. The off-diagonal elements of the matrix are equal
to $0$. For each fit, we also calculated $\chi^{2}$ in the form as presented
by Equation (9).

PLRs need to be calibrated using objects with known and well-measured
distances. The distance to the LMC was measured with $1\%$ accuracy by
Pietrzyński et al. (2019) and it is $d=49.59\pm 0.09$ (statistical) $\pm$ 0.54
(systematic) kpc (distance modulus $\mu=18.477$). We calculated the absolute
zero-point $a_{\mathrm{0,\lambda,abs}}$ as the difference
$a_{\mathrm{0,\lambda}}-\mu$, while the uncertainty of the absolute zero point
is the sum of the zero-point $a_{\mathrm{0,\lambda}}$ uncertainty resulting
from the fit, and both statistical and systematic uncertainties of the
distance, all combined in quadrature.

###  4.2 PLC Relations for C-rich Miras

A strong correlation between colors and pulsation periods of Mira variables
exists in the optical and NIR bands. Therefore, the PLC usually has much lower
scatter than the PLRs, especially for C-rich Miras (see e.g. Bhardwaj et al.,
2019). The PLC can be described in the form of the following equation:

$m_{\mathrm{\lambda,fit}}=a_{\mathrm{0,\lambda}}+a_{\mathrm{1,\lambda}}\times(\log{P}-2.3)+c\times(m_{\mathrm{\lambda_{1}}}-m_{\mathrm{\lambda_{2}}}),$
(10)

where $c\times(m_{\mathrm{\lambda_{1}}}-m_{\mathrm{\lambda_{2}}})$ is the
color-term. In the case of optical and NIR PLC relations, the coefficient $c$
is statistically significant with relatively small uncertainty (see Table 2 in
Bhardwaj et al., 2019).

We fitted PLC relations (Equation (10)) to the C-rich Miras in the WISE bands
using the $(W1-W2)$ color index, while in the Spitzer bands we tried fits with
two color indices, i.e. $([3.6]-[4.5])$ and $([3.6]-[5.8])$. In each case, the
fit including the color-term resulted in a statistically insignificant $c$
coefficient – $c$ was close to zero with a comparable uncertainty. This led us
to the conclusion that in the mid-IR the dependence between the pulsation
period and color does not exist, or is very small.

To investigate this topic further, we used magnitudes in NIR bands ($J$ and
$K$) extracted from the 2MASS All-Sky Catalog of Point Sources (Cutri et al.,
2003). We searched for objects within 1 arcsec radius around LMC Miras’
coordinates. We found 1656 (out of 1663) Miras’ counterparts. We corrected $J$
and $K$ magnitudes for the interstellar extinction using the same method as
described in Section 3.1, with transformations between bands derived by Chen
et al. (2018) and Wang & Chen (2019).

In Figure 2, we present color-magnitude diagrams (CMDs, upper row) for all
Miras plotted with NIR and mid-IR data mentioned above. The color map in each
case presents the pulsation period of Miras. In the lower row of Figure 2, we
present period–color diagrams using the same color-term as in CMDs. In the
period–color diagrams, we plotted C-rich (red) and O-rich (blue) Miras,
separately.

From the CMDs it is visible why the color-term must have an influence on the
PLRs scatter in optical and NIR bands and why it is negligible in the mid-IR.
For the $K$-band, we can find very bright and blue Miras with the exact same
periods as much fainter and red Miras. This is not the case for the mid-IR
data, where the period seems to depend on the brightness only.

In the left plot of the lower panel in Figure 2, we present period-color
diagram $(J-K)$ vs. $\log P$, which is commonly used to divide AGB stars into
C- and O-rich (Soszyński et al., 2009; Ita & Matsunaga, 2011). It is visible
that both groups of Miras form two distinct groups. A similar division is also
seen for mid-IR colors (middle and right plot), but it seems that the overlap
of both groups is slightly larger than that in NIR. However, the mid-IR
period–color diagram could also be used for the division of AGB stars into
C-rich and O-rich groups in other galaxies.

Table 4: Mid-IR Period-Luminosity Relations for C-rich and O-rich Miras in
WISE and Spitzer bands.

C-rich Mirasa
---
$\lambda$ | $a_{\mathrm{0,\lambda}}$ | $a_{\mathrm{0,\lambda,abs}}^{\dagger}$ | $a_{\mathrm{1,\lambda}}$ | $a_{\mathrm{2,\lambda}}$ | $\sigma_{\mathrm{\lambda}}^{c}$ | $\chi^{2}/\mathrm{dof}^{d}$ | $N_{\mathrm{in}}$ | $N_{\mathrm{fin}}$ | $N_{\mathrm{out}}$
WISE W1 (Fig. 3a) | $11.006\pm 0.023$ | $-7.471\pm 0.033$ | $-4.198\pm 0.074$ | $0.000$ | $0.258$ | $0.97$ | $999$ | $982$ | $17$
WISE W2 (Fig. 3b) | $10.916\pm 0.030$ | $-7.561\pm 0.038$ | $-6.445\pm 0.094$ | $0.000$ | $0.326$ | $0.97$ | $998$ | $986$ | $12$
WISE W3 (Fig. 3c) | $10.545\pm 0.052$ | $-7.932\pm 0.058$ | $-9.040\pm 0.164$ | $0.000$ | $0.584$ | $0.99$ | $1018$ | $1012$ | $6$
Spitzer $[3.6]$ (Fig. 6a) | $10.734\pm 0.028$ | $-7.743\pm 0.037$ | $-4.147\pm 0.091$ | $0.000$ | $0.344$ | $0.89$ | $1072$ | $1064$ | $8$
Spitzer $[4.5]$ (Fig. 6b) | $10.775\pm 0.033$ | $-7.702\pm 0.041$ | $-5.780\pm 0.106$ | $0.000$ | $0.406$ | $0.93$ | $1072$ | $1066$ | $6$
Spitzer $[5.8]$ (Fig. 6c) | $10.789\pm 0.040$ | $-7.689\pm 0.047$ | $-7.183\pm 0.130$ | $0.000$ | $0.494$ | $0.76$ | $1072$ | $1066$ | $6$
Spitzer $[8.0]$ (Fig. 6d) | $10.691\pm 0.044$ | $-7.786\pm 0.050$ | $-8.329\pm 0.143$ | $0.000$ | $0.541$ | $0.96$ | $1072$ | $1065$ | $7$
O-rich Miras
linear fita, $\log P\leq 2.6$
WISE W1 (Fig. 4a) | $11.228\pm 0.009$ | $-7.249\pm 0.025$ | $-3.807\pm 0.066$ | $0.000$ | $0.116$ | $0.55$ | $196$ | $194$ | $2$
WISE W2 (Fig. 4b) | $11.030\pm 0.013$ | $-7.447\pm 0.027$ | $-3.794\pm 0.096$ | $0.000$ | $0.170$ | $0.80$ | $196$ | $193$ | $2$
WISE W3 (Fig. 4c) | $10.251\pm 0.028$ | $-8.226\pm 0.036$ | $-3.900\pm 0.210$ | $0.000$ | $0.365$ | $0.94$ | $199$ | $194$ | $5$
Spitzer $[3.6]$ (Fig. 7a) | $11.073\pm 0.013$ | $-7.404\pm 0.027$ | $-3.522\pm 0.103$ | $0.000$ | $0.236$ | $0.73$ | $342$ | $327$ | $15$
Spitzer $[4.5]$ (Fig. 7b) | $10.970\pm 0.014$ | $-7.507\pm 0.028$ | $-3.290\pm 0.111$ | $0.000$ | $0.256$ | $0.78$ | $342$ | $328$ | $14$
Spitzer $[5.8]$ (Fig. 7c) | $10.792\pm 0.014$ | $-7.685\pm 0.028$ | $-3.313\pm 0.111$ | $0.000$ | $0.254$ | $0.78$ | $342$ | $323$ | $19$
Spitzer $[8.0]$ (Fig. 7d) | $10.595\pm 0.017$ | $-7.882\pm 0.029$ | $-3.473\pm 0.131$ | $0.000$ | $0.303$ | $0.82$ | $342$ | $330$ | $12$
quadratic fitb
WISE W1 (Fig. 5a) | $11.252\pm 0.010$ | $-7.224\pm 0.026$ | $-3.749\pm 0.075$ | $-1.810\pm 0.154$ | $0.128$ | $0.92$ | $238$ | $227$ | $11$
WISE W2 (Fig. 5b) | $11.069\pm 0.013$ | $-7.408\pm 0.027$ | $-3.679\pm 0.102$ | $-2.911\pm 0.207$ | $0.173$ | $0.95$ | $238$ | $231$ | $7$
WISE W3 (Fig. 5c) | $10.335\pm 0.029$ | $-8.142\pm 0.037$ | $-3.459\pm 0.224$ | $-6.107\pm 0.422$ | $0.383$ | $0.91$ | $240$ | $232$ | $8$
Spitzer $[3.6]$ (Fig. 8a) | $11.094\pm 0.015$ | $-7.383\pm 0.028$ | $-3.648\pm 0.115$ | $-2.201\pm 0.239$ | $0.267$ | $0.74$ | $398$ | $381$ | $17$
Spitzer $[4.5]$ (Fig. 8b) | $11.000\pm 0.016$ | $-7.477\pm 0.029$ | $-3.430\pm 0.123$ | $-2.895\pm 0.257$ | $0.287$ | $0.84$ | $398$ | $385$ | $13$
Spitzer $[5.8]$ (Fig. 8c) | $10.813\pm 0.017$ | $-7.664\pm 0.029$ | $-3.568\pm 0.128$ | $-2.736\pm 0.256$ | $0.300$ | $0.81$ | $398$ | $385$ | $13$
Spitzer $[8.0]$ (Fig. 8d) | $10.640\pm 0.018$ | $-7.837\pm 0.030$ | $-4.037\pm 0.270$ | $-4.363\pm 0.201$ | $0.315$ | $0.82$ | $398$ | $385$ | $13$

a linear model in a form:
$M_{\mathrm{\lambda,fit}}=a_{\mathrm{0,\lambda}}+a_{\mathrm{1,\lambda}}\times(\log{P}-2.3)$
b quadratic model in a form:
$M_{\mathrm{\lambda,fit}}=a_{{0,\mathrm{\lambda}}}+a_{{\mathrm{1,\lambda}}}\times(\log{P}-2.3)+a_{{\mathrm{2,\lambda}}}\times(\log{P}-2.3)^{2}$
c scatter in the period–luminosity relation
d $\chi^{2}$ calculated using Equation 9 divided by the degrees of freedom
† absolute zero-point with uncertainty calculated using distance modulus
$\mu=18.477$ mag (Pietrzyński et al., 2019), and systematic and statistic
uncertainties of the distance to the LMC

##  5 Results

Using the PLR models and fitting procedures, described in Section 4.1, we
performed fits to the datasets using the weighted least squares method with a
$3\sigma$-clipping procedure. We iteratively fitted the linear model in a form
of Equation (6) to the O-rich and C-rich Miras, and the quadratic model in a
form of Equation (7) to the O-rich Miras only. In each iteration, we rejected
points deviating by more than $3\sigma$ from the model, until no outliers were
left. The final parameters of the fitted models (mid-IR PLRs), the absolute
zero-point $a_{\mathrm{0,\lambda,abs}}$ with its uncertainty, are presented in
Table LABEL:tab:PL_results, along with dispersions
$\sigma_{\mathrm{\lambda}}$, $\chi^{2}/\mathrm{dof}$, the number of initial
data points $N_{\mathrm{in}}$, the final number of data points after rejecting
outliers $N_{\mathrm{fin}}$, and the number of outliers rejected during
$3\sigma$-clipping procedure $N_{\mathrm{out}}$. PLRs for each band are
presented in Figures: 3 (WISE, C-rich Miras), 4 (WISE, O-rich Miras, linear
fit), 5 (WISE, O-rich Miras, quadratic fit), 6 (Spitzer, C-rich Miras), 7
(Spitzer, O-rich Miras, linear fit), 8 (Spitzer, O-rich Miras, quadratic fit).

As mentioned in Section 2.3, the SEDs fitting showed that photometric
measurements in the $W4$-band are significantly biased. Therefore, in Figures
3, 4, and 5, we present the PLRs in $W4-$band for completeness purposes only.
We do not provide the fit parameters in Table LABEL:tab:PL_results and in
Figures 3, 4, and 5, because PLRs based on biased photometric magnitudes lead
to biased PLRs, and should not be used in future studies.

PLRs for the C-rich Miras, both in the WISE and Spitzer bands, have
significantly larger dispersions $\sigma_{\mathrm{\lambda}}$ than the ones
obtained for the O-rich Mira PLRs. One of the reasons for that is the presence
of gas and dust ejecta from these stars that can significantly alter their
luminosity (by introducing long-term trends or variations), which in turn can
lead to larger dispersions. As a rule of thumb, changes of brightness (the
mean and amplitude) of O-rich Miras are reasonably stable over time, so that
scatter along PLRs is smaller than for the C-rich Miras.

###  5.1 Comparison of mean magnitude measuring methods

One of the basic parameters that characterizes a star is their mean magnitude.
The mean magnitude should not depend on the measurement method, however,
various methods come with their own advantages and disadvantages, the latter
leading to biases. Since the mean magnitudes are extremely important for the
PLRs analyses, we investigated the influence of several methods on the PLRs
zero-point $a_{\mathrm{0,\lambda}}$, slope $a_{\mathrm{1,\lambda}}$, and
dispersion $\sigma_{\lambda}$.

The PLRs in this paper (Figures 3–8 and Table LABEL:tab:PL_results) are based
on the mean magnitudes measured from the fitted and scaled to mid-IR OGLE
templates, as the integrated third-order Fourier series (see Section 3.2 in
this paper, and Iwanek et al. 2021 for more details). In principle, in some
cases this method could fail, as the Mira brightness in addition to typical
periodic variability may exhibit aperiodic variations of the mean magnitude
and amplitude. We compared these PLRs (in each Spitzer and WISE band) with
PLRs constructed using mean magnitudes measured in two different ways. The
PLRs fitting method is the same and is described in Section 4.1.

In contrast to more sophisticated methods (i.e., the integrated Fourier series
fit), we also measured the mean magnitudes using two simple methods: (1) the
regular flux-based mean of the scaled OGLE templates to the mid-IR data, and
(2) the sum of maximum and minimum fluxes divided by two and converted to
magnitude. The magnitudes have been corrected for the interstellar extinction
as described in Section 3.1. We then constructed PLRs separately for the
C-rich Miras (in the full range of pulsation periods) and for the O-rich Miras
(with linear fit to periods $\log P\leq 2.6$).

Let us compare the differences between the PLR parameters for the mean
magnitudes obtained with different methods to the corresponding uncertainties
of our base models (presented in Table LABEL:tab:PL_results). In both cases
(base—method (1) and base—method (2)), the smallest differences from our base
models were obtained for the mean magnitudes calculated as a regular mean
(method (1)). For the C-rich Miras, the differences in zero-points are at the
level of $0.5\sigma$, while for the O-rich Miras these differences are less
than $0.5\sigma$, for both Spitzer and WISE data. The slope of PLRs varies
around $3\sigma$ for our base models of C-rich Miras, while this differences
for the O-rich Miras are smaller than $0.5\sigma$. Dispersions are the
smallest, almost in all cases, for the mean magnitudes obtained from the
integration of the Fourier series and are used in PLRs presented in Table
LABEL:tab:PL_results and in Figures 3–8. Method (2) significantly
underestimates the mean magnitudes and causes much worse results than the base
model.

##  6 Conclusions

In this paper, we used the catalog of Mira-type variable stars in the LMC from
the OGLE survey (Soszyński et al., 2009) to construct the most accurate mid-IR
PLRs to date. We cross-matched the known LMC Miras with the WISE (Wright et
al., 2010) and Spitzer databases (Werner et al., 2004) to obtain mid-IR light
curves. We used two-decades-long OGLE light curves to create templates using
the GPR model (He et al., 2016), and we fitted these templates to the mid-IR
data (Iwanek et al., 2021). We measured the mean magnitudes in mid-IR WISE and
Spitzer bands by fitting a third-order truncated Fourier series, and
integrating them. The mean magnitudes were corrected for the interstellar
extinction using the reddening map of the LMC (Skowron et al., 2021), and
extinction curves in the mid-IR bands (Chen et al., 2018; Wang & Chen, 2019).
We calibrated the Mira PLRs using the most accurate distance measurement to
the LMC (Pietrzyński et al., 2019).

The PLRs presented in Table LABEL:tab:PL_results can be used to measure
distances to the Mira-type stars in the Milky Way and nearby galaxies. The
distance to an individual Mira star can be measured with the accuracy at the
level of $5\%$ and $12\%$ for the O-rich and C-rich Mira, respectively.

In future work, we are going to measure the mid-IR distances to the tens of
thousands of Mira-type stars in the Milky Way discovered in the OGLE data. It
will allow us to study the three-dimensional structure of the Milky Way.

We are deeply grateful to the anonymous referee for the detailed report. The
referee raised several points that turned out to be of significant importance.
The suggestions given by the referee let us to greatly improved this
manuscript, but more importantly, these comments motivated us to explored the
subject of Miras variability in-depth. As a result, the referee’s remarks led
to the creation of a second paper entitled ”Multiwavelength Properties of
Miras” (Iwanek et al., 2021). We thank Drs. Jan Skowron, Przemek Mróz, Radek
Poleski, and Mariusz Gromadzki, and Marcin Wrona for comments that helped to
improve this manuscript. This work has been supported by the National Science
Centre, Poland, grant MAESTRO no. 2016/22/A/ST9/00009 to IS. PI is partially
supported by the Kartezjusz programme no. POWR.03.02.00-00-I001/16-00 founded
by the National Centre for Research and Development, Poland. SK acknowledges
the financial support of the Polish National Science Center through grant no.
2018/31/B/ST9/00334. This publication makes use of data products from the
Wide-field Infrared Survey Explorer (WISE), which is a joint project of the
University of California, Los Angeles, and the Jet Propulsion
Laboratory/California Institute of Technology, funded by the National
Aeronautics and Space Administration (NASA). This work is based in part on
archival data obtained with the Spitzer Space Telescope, which was operated by
the Jet Propulsion Laboratory, California Institute of Technology under a
contract with NASA.

Figure 3: Period–Luminosity Relations for C-rich Miras from the LMC in mid-IR
WISE bands: (a) W1, (b) W2, (c) W3, (d) W4. In each case, the left (right)
y-axis of the top panel presents the observed (absolute) mean magnitude in a
given band $m_{\mathrm{{\lambda}}}$ as a function of the logarithm of the
pulsation period $P$ (in days), while the bottom panel presents residuals
calculated as the difference between the measured and predicted absolute
magnitudes from the fit. In the top panel, we present the best fit with the
solid line, while dashed lines present the range of $\pm 3\sigma$, where
$\sigma$ means the dispersion along the fit. The model is given in each plot
(with the exception of panel (d) for the $W4$-band due to the photometric
issues in this band) in the form:
$m_{\mathrm{\lambda,fit}}=a_{\mathrm{0,\lambda}}+a_{\mathrm{1,\lambda}}\times(\log{P}-2.3)$.
Final samples of C-rich Miras are plotted with red points. Errorbars of
$\log{P}$ are smaller than the point size. Rejected stars during
$3\sigma$-clipping procedure are plotted in gray. The final $3\sigma$ range,
above which there are no outliers, is plotted as dashed line in the bottom
plot. Panel (d) for the $W4-$band is shown here for completeness.

Figure 4: Period–Luminosity Relations for O-rich Miras from the LMC in mid-IR
WISE bands: (a) W1, (b) W2, (c) W3, (d) W4. In each case, the left (right)
y-axis of the top panel presents the observed (absolute) mean magnitude in a
given band $m_{\mathrm{{\lambda}}}$ as a function of the logarithm of the
pulsation period $P$ (in days), while the bottom panel presents residuals
calculated as the difference between the measured and predicted absolute
magnitudes from the fit. In the top panel, we present the best fit with the
solid line, while dashed lines present the range of $\pm 3\sigma$, where
$\sigma$ means the dispersion along the fit. The model is given in each plot
(with the exception of panel (d) for $W4$-band due to the photometric issues
in this band) in the form:
$m_{\mathrm{\lambda,fit}}=a_{{0,\mathrm{\lambda}}}+a_{{\mathrm{1,\lambda}}}\times(\log{P}-2.3)$.
Final samples of O-rich Miras are plotted with blue points. Errorbars of
$\log{P}$ are smaller than the point size. Rejected stars during
$3\sigma$-clipping procedure are plotted in gray. The final $3\sigma$ range,
above which there are no outliers, is plotted as dashed line in the bottom
plot. Panel (d) for the $W4-$band is shown here for completeness.

Figure 5: Period–Luminosity Relations for O-rich Miras from the Large
Magellanic Cloud in mid-IR WISE bands: (a) W1, (b) W2, (c) W3, (d) W4. In each
case, the left (right) y-axis of the top panel presents the observed
(absolute) mean magnitude in a given band $m_{\mathrm{{\lambda}}}$ as a
function of the logarithm of the pulsation period $P$ (in days), while the
bottom panel presents residuals calculated as the difference between the
measured and predicted absolute magnitudes from the fit. In the top panel, we
present the best fit with the solid line, while dashed lines present the range
of $\pm 3\sigma$, where $\sigma$ means the dispersion along the fit. The model
is given in each plot (with the exception of panel (d) for $W4$-band due to
the photometric issues in this band) in the form:
$m_{\mathrm{\lambda,fit}}=a_{{0,\mathrm{\lambda}}}+a_{{\mathrm{1,\lambda}}}\times(\log{P}-2.3)+a_{{\mathrm{2,\lambda}}}\times(\log{P}-2.3)^{2}$.
Final samples of O-rich Miras are plotted with blue points. Errorbars of
$\log{P}$ are smaller than the point size. Rejected stars during
$3\sigma$-clipping procedure are plotted in gray. The final $3\sigma$ range,
above which there are no outliers, is plotted as dashed line in the bottom
plot. Panel (d) for the $W4-$band is shown here for completeness.

Figure 6: Same as in Figure 3, but for the Spitzer bands: (a) $[3.6]\;\mu m$,
(b) $[4.5]\;\mu m$, (c) $[5.8]\;\mu m$, (d) $[8.0]\;\mu m$. Figure 7: Same as
in Figure 4, but for the Spitzer bands: (a) $[3.6]\;\mu m$, (b) $[4.5]\;\mu
m$, (c) $[5.8]\;\mu m$, (d) $[8.0]\;\mu m$. Figure 8: Same as in Figure 5, but
for the Spitzer bands: (a) $[3.6]\;\mu m$, (b) $[4.5]\;\mu m$, (c) $[5.8]\;\mu
m$, (d) $[8.0]\;\mu m$.

## References

  * Alcock et al. (2000) Alcock, C., Allsman, R. A., Alves, D. R., et al. 2000, The Astrophysical Journal, 542, 281, doi: 10.1086/309512
  * Ansari (1996) Ansari, R. 1996, Vistas in Astronomy, 40, 519, doi: 10.1016/S0083-6656(96)00037-2
  * Bhardwaj et al. (2019) Bhardwaj, A., Kanbur, S., He, S., et al. 2019, ApJ, 884, 20, doi: 10.3847/1538-4357/ab38c2
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Chen et al. (2018) Chen, X., Wang, S., Deng, L., & de Grijs, R. 2018, ApJ, 859, 137, doi: 10.3847/1538-4357/aabfbc
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, 2MASS All Sky Catalog of point sources.
  * Feast et al. (1989) Feast, M. W., Glass, I. S., Whitelock, P. A., & Catchpole, R. M. 1989, MNRAS, 241, 375, doi: 10.1093/mnras/241.3.375
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gerasimovič (1928) Gerasimovič, B. P. 1928, Proceedings of the National Academy of Science, 14, 963, doi: 10.1073/pnas.14.12.963
  * Glass & Lloyd Evans (1981) Glass, I. S., & Lloyd Evans, T. 1981, Nature, 291, 303, doi: 10.1038/291303a0
  * Gould (2003) Gould, A. 2003, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0310577
  * He et al. (2016) He, S., Yuan, W., Huang, J. Z., Long, J., & Macri, L. M. 2016, AJ, 152, 164, doi: 10.3847/0004-6256/152/6/164
  * Huang et al. (2018) Huang, C. D., Riess, A. G., Hoffmann, S. L., et al. 2018, ApJ, 857, 67, doi: 10.3847/1538-4357/aab6b3
  * Huang et al. (2020) Huang, C. D., Riess, A. G., Yuan, W., et al. 2020, ApJ, 889, 5, doi: 10.3847/1538-4357/ab5dbd
  * Hughes & Wood (1990) Hughes, S. M. G., & Wood, P. R. 1990, AJ, 99, 784, doi: 10.1086/115374
  * Iben & Renzini (1983) Iben, I., J., & Renzini, A. 1983, ARA&A, 21, 271, doi: 10.1146/annurev.aa.21.090183.001415
  * Ita & Matsunaga (2011) Ita, Y., & Matsunaga, N. 2011, MNRAS, 412, 2345, doi: 10.1111/j.1365-2966.2010.18056.x
  * Iwanek et al. (2021) Iwanek, P., Kozłowski, S., Gromadzki, M., et al. 2021, arXiv e-prints, arXiv:2107.03397.  https://arxiv.org/abs/2107.03397
  * Leavitt & Pickering (1912) Leavitt, H. S., & Pickering, E. C. 1912, Harvard College Observatory Circular, 173, 1
  * Mainzer et al. (2011) Mainzer, A., Bauer, J., Grav, T., et al. 2011, ApJ, 731, 53, doi: 10.1088/0004-637X/731/1/53
  * Mainzer et al. (2014) Mainzer, A., Bauer, J., Cutri, R. M., et al. 2014, ApJ, 792, 30, doi: 10.1088/0004-637X/792/1/30
  * Matsunaga et al. (2009) Matsunaga, N., Kawadu, T., Nishiyama, S., et al. 2009, MNRAS, 399, 1709, doi: 10.1111/j.1365-2966.2009.15393.x
  * Meixner et al. (2006) Meixner, M., Gordon, K. D., Indebetouw, R., et al. 2006, AJ, 132, 2268, doi: 10.1086/508185
  * Molina et al. (2019) Molina, C. N., Borissova, J., Catelan, M., et al. 2019, MNRAS, 482, 5567, doi: 10.1093/mnras/sty3041
  * Perrin et al. (2020) Perrin, G., Ridgway, S. T., Lacour, S., et al. 2020, A&A, 642, A82, doi: 10.1051/0004-6361/202037443
  * Pietrzyński et al. (2019) Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Nature, 567, 200, doi: 10.1038/s41586-019-0999-4
  * Qin et al. (2018) Qin, W., Nataf, D. M., Zakamska, N., Wood, P. R., & Casagrande, L. 2018, ApJ, 865, 47, doi: 10.3847/1538-4357/aad7fb
  * Riebel et al. (2010) Riebel, D., Meixner, M., Fraser, O., et al. 2010, ApJ, 723, 1195, doi: 10.1088/0004-637X/723/2/1195
  * Riebel et al. (2015) Riebel, D., Boyer, M. L., Srinivasan, S., et al. 2015, ApJ, 807, 1, doi: 10.1088/0004-637X/807/1/1
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schwarzenberg-Czerny (1996) Schwarzenberg-Czerny, A. 1996, ApJ, 460, L107, doi: 10.1086/309985
  * Skowron et al. (2021) Skowron, D. M., Skowron, J., Udalski, A., et al. 2021, ApJS, 252, 23, doi: 10.3847/1538-4365/abcb81
  * Soszyński et al. (2004) Soszyński, I., Udalski, A., Kubiak, M., et al. 2004, Acta Astron., 54, 129.  https://arxiv.org/abs/astro-ph/0407057
  * Soszyński et al. (2005) —. 2005, Acta Astron., 55, 331.  https://arxiv.org/abs/astro-ph/0512578
  * Soszyński et al. (2009) Soszyński, I., Udalski, A., Szymański, M. K., et al. 2009, Acta Astron., 59, 239.  https://arxiv.org/abs/0910.1354
  * Udalski (2003) Udalski, A. 2003, Acta Astron., 53, 291.  https://arxiv.org/abs/astro-ph/0401123
  * Udalski et al. (2015) Udalski, A., Szymański, M. K., & Szymański, G. 2015, Acta Astron., 65, 1\.  https://arxiv.org/abs/1504.05966
  * Urago et al. (2020) Urago, R., Omodaka, T., Nagayama, T., et al. 2020, ApJ, 891, 50, doi: 10.3847/1538-4357/ab70b1
  * Wang & Chen (2019) Wang, S., & Chen, X. 2019, ApJ, 877, 116, doi: 10.3847/1538-4357/ab1c61
  * Werner et al. (2004) Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1, doi: 10.1086/422992
  * Whitelock (2012) Whitelock, P. A. 2012, Ap&SS, 341, 123, doi: 10.1007/s10509-012-0993-x
  * Whitelock et al. (2006) Whitelock, P. A., Feast, M. W., Marang, F., & Groenewegen, M. A. T. 2006, MNRAS, 369, 751, doi: 10.1111/j.1365-2966.2006.10322.x
  * Whitelock et al. (2008) Whitelock, P. A., Feast, M. W., & Van Leeuwen, F. 2008, MNRAS, 386, 313, doi: 10.1111/j.1365-2966.2008.13032.x
  * Whitelock et al. (2017) Whitelock, P. A., Kasliwal, M., & Boyer, M. 2017, in European Physical Journal Web of Conferences, Vol. 152, European Physical Journal Web of Conferences, 01009, doi: 10.1051/epjconf/201715201009
  * Wilson & Merrill (1942) Wilson, R. E., & Merrill, P. W. 1942, ApJ, 95, 248, doi: 10.1086/144391
  * Wood (2000) Wood, P. R. 2000, PASA, 17, 18, doi: 10.1071/AS00018
  * Wood et al. (1999) Wood, P. R., Alcock, C., Allsman, R. A., et al. 1999, in Asymptotic Giant Branch Stars, ed. T. Le Bertre, A. Lebre, & C. Waelkens, Vol. 191, 151
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Yuan et al. (2017a) Yuan, W., He, S., Macri, L. M., Long, J., & Huang, J. Z. 2017a, AJ, 153, 170, doi: 10.3847/1538-3881/aa63f1
  * Yuan et al. (2017b) Yuan, W., Macri, L. M., He, S., et al. 2017b, AJ, 154, 149, doi: 10.3847/1538-3881/aa86f1
  * Yuan et al. (2018) Yuan, W., Macri, L. M., Javadi, A., Lin, Z., & Huang, J. Z. 2018, AJ, 156, 112, doi: 10.3847/1538-3881/aad330