# The period–luminosity relation for Mira variables in the Milky Way using
Gaia DR3: a further distance anchor for $H_{0}$

Jason L. Sanders
Department of Physics and Astronomy, University College London, London WC1E
6BT, UK
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Gaia DR3 parallaxes are used to calibrate preliminary period–luminosity
relations of O-rich Mira variables in the 2MASS $J$, $H$ and $K_{s}$ bands
using a probabilistic model accounting for variations in the parallax
zeropoint and underestimation of the parallax uncertainties. The derived
relations are compared to those measured for the Large and Small Magellanic
Clouds, the Sagittarius dwarf spheroidal galaxy, globular cluster members and
the subset of Milky Way Mira variables with VLBI parallaxes. The Milky Way
linear $JHK_{s}$ relations are slightly steeper and thus fainter at short
period than the corresponding LMC relations suggesting population effects in
the near-infrared are perhaps larger than previous observational works have
claimed. Models of the Gaia astrometry for the Mira variables suggest that,
despite the intrinsic photocentre wobble and use of mean photometry in the
astrometric solution of the current data reduction, the recovered parallaxes
should be on average unbiased but with underestimated uncertainties for the
nearest stars. The recommended Gaia EDR3 parallax zeropoint corrections
evaluated at $\nu_{\mathrm{eff}}=1.25\,\mu\mathrm{m}^{-1}$ require minimal
($\lesssim 5\,\mu\mathrm{as}$) corrections for redder five-parameter sources,
but over-correct the parallaxes for redder six-parameter sources, and the
parallax uncertainties are underestimated, at most by a factor $\sim 1.6$ at
$G\approx 12.5\,\mathrm{mag}$. The derived period–luminosity relations are
used as anchors for the Mira variables in the Type Ia host galaxy NGC 1559 to
find $H_{0}=(73.7\pm 4.4)\,\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$.

###### keywords:

stars: AGB – stars: variables: general – stars: distances – cosmological
parameters

††pubyear: 2021††pagerange: The period–luminosity relation for Mira variables
in the Milky Way using Gaia DR3: a further distance anchor for $H_{0}$–C.2

##  1 Introduction

Mira variables are thermally pulsating asymptotic giant branch (AGB) stars
with characteristic periods of between $100$ and $1000$ days, and high
amplitudes ($\gtrsim 2.5$ in $V$ and between $\sim 0.3$ and $\sim 1$ in
$K_{s}$, Matsunaga et al., 2009; Catelan & Smith, 2015). Primarily from their
study in the Large Magellanic Cloud (LMC, Glass & Evans, 1981; Wood et al.,
1999; Soszyński et al., 2009), they are known to follow period–luminosity
relations (with a typical scatter of $\sim 0.2\,\mathrm{mag}$ from single-
epoch $K_{s}$ data and $\sim 0.1\,\mathrm{mag}$ for mean $K_{s}$ measurements,
Yuan et al., 2017b). As AGB stars, Mira variables have chemistry dominated by
either carbon-rich or oxygen-rich species as determined by the strength of
dredge-up episodes, largely a reflection of their initial mass and composition
(Höfner & Olofsson, 2018). Both C-rich and O-rich Mira variables satisfy
period–luminosity relations (e.g. the recent calibrations from Iwanek et al.,
2021a) although the O-rich relations are typically tighter than the C-rich
relations in the near-infrared due to the presence of significant
circumstellar dust in the C-rich Mira variables (Ita & Matsunaga, 2011). This
makes O-rich Mira variables powerful distance tracers for both Galactic and
cosmological studies.

The need for reliable well-calibrated distance indicators has received
significant recent interest in light of the ‘Hubble tension’. The current
expansion rate of the Universe, the Hubble constant $H_{0}$, can be measured
using Type Ia supernovae in nearby galaxies or alternatively extrapolated from
the early Universe using the best-fitting $\Lambda$CDM model of the cosmic
microwave background radiation (Planck Collaboration et al., 2014). An
absolute calibration, or anchor, of the Hubble diagram is required to utilise
the Type Ia supernovae, and traditionally the most precise and well-studied
calibrators have been the classical Cepheids (Freedman et al., 2001; Riess et
al., 2011; Riess et al., 2021; Riess et al., 2022a). The problem then becomes
anchoring the Cepheid scale, which can be done with local Cepheids using Gaia
parallax measurements (Gaia Collaboration et al., 2021) of individual Cepheids
(Riess et al., 2021) or those of their host cluster (Riess et al., 2022b),
eclipsing binaries in the Magellanic Clouds (Pietrzyński et al., 2019; Graczyk
et al., 2020) or the water maser in NGC 4258 (Reid et al., 2019). The latest
estimates of the Hubble constant from Riess et al. (2022a, b) using classical
Cepheids with a combination of all three anchors are in tension at the $\sim
5\sigma$ level with the early Universe extrapolation from Planck Collaboration
et al. (2014) possibly pointing towards new physics beyond the standard
cosmological model (Di Valentino et al., 2021). However, the discrepancy could
also arise from systematics in the use of Cepheids (e.g. Efstathiou, 2020).
There have been many proposed and applied alternatives to classical Cepheids
such as the tip of the giant branch (e.g. Freedman, 2021) which produces a
more intermediate result between that of Planck Collaboration et al. (2014)
and Riess et al. (2021), the J-AGB method (Madore & Freedman, 2020),
gravitational lensing (Wong et al., 2020) and masers (Pesce et al., 2020).
Mira variables offer another interesting alternative to the usual classical
Cepheid variables as (i) they are less biased to young populations so are
present in a broad range of galaxies, in particular the full range of Type Ia
supernovae hosting galaxies, (ii) as intermediate age tracers they are likely
not in crowded or dust-obscured regions of their host galaxies so the
photometric systematics are weaker, and (iii) they can be brighter than
Cepheid variables in the infrared so can be utilised in more distant galaxies,
especially in the era of the James Webb Space Telescope. Recently Huang et al.
(2020) have used a sample of Mira variables in NGC 1559 anchored to Mira
variables in the LMC and/or NGC 4258 to estimate the distance to SN 2005df and
measure $H_{0}=(73.3\pm 4.0)\,\mathrm{km\,s}^{-1}\mathrm{kpc}^{-1}$ in good
agreement with other local measurements (as well as the early Universe
extrapolated value from Planck Collaboration et al., 2014, at the $\sim
1.5\sigma$ level).

Mira variables have also found significant use as a tracer of Galactic and
Local Group structure. Thanks to their brightness in the infrared and their
representation across a range of intermediate age populations, they are useful
probes of structure across the Galactic disc (Feast & Whitelock, 2000b; Grady
et al., 2019, 2020), the Galactic bulge (Catchpole et al., 2016), the heavily-
extincted nuclear stellar region (Glass et al., 2001; Matsunaga et al., 2009;
Sanders et al., 2022) and the Magellanic Clouds (e.g. Deason et al., 2017).
Furthermore, their periods are linked to their age (and possibly metallicity),
as confirmed empirically by variations of velocity dispersion with period
(Feast & Whitelock, 2000b) and demonstrated theoretically in non-linear
pulsation calculations (Trabucchi & Mowlavi, 2022). Recently, Grady et al.
(2020) have used the empirical period–age relation for O-rich Mira variables
to map the age structure of the Milky Way’s bar-bulge and disc.

Typically, the period–luminosity relation of Mira variables has been
calibrated using Mira variables in the LMC (Glass & Evans, 1981; Feast et al.,
1989; Ita et al., 2004; Groenewegen, 2004; Fraser et al., 2008; Riebel et al.,
2010; Ita & Matsunaga, 2011; Yuan et al., 2017a, b; Bhardwaj et al., 2019;
Iwanek et al., 2021b). However, population effects (e.g. metallicity and age
variations) can alter the period–luminosity relation (e.g. Qin et al., 2018).
For both extragalactic and Galactic studies, a calibration based on the
perhaps more representative Milky Way Mira variables could be preferable.
Whitelock et al. (2008) used a sample of $184$ O-rich Mira variables observed
by the Hipparcos satellite in combination with Mira variables in globular
clusters and those observed with VLBI to derive a near-infrared
period–luminosity relation of $M_{K}=(-7.25\pm 0.07)+3.50(\log_{10}P-2.38)$,
within $\sim 0.02$ of their derived LMC relation (correcting for the updated
LMC distance modulus from Pietrzyński et al., 2019). This already suggests the
population effects on the ($K_{s}$-band) period–luminosity relation are small.

The arrival of data from the Gaia satellite (Gaia Collaboration et al., 2016,
2018, 2021) has opened up the possibility of an updated fully geometric
calibration of the Mira period–luminosity relation, particularly as Gaia’s
multi-epoch observations have enabled all-sky catalogues of Mira variables to
be extracted from the data (Mowlavi et al., 2018; Lebzelter et al., 2022).
However, significant care must be taken when using astrometric data. Large
parallax uncertainties can introduce a Lutz-Kelker bias when converting
parallax measurements to distances which must be avoided with more careful
probabilistic inversions (e.g. Luri et al., 2018; Bailer-Jones et al., 2018).
Furthermore, systematic variations in the Gaia parallax zeropoint are present
at the $\sim 10\,\mu\mathrm{as}$ level and vary with magnitude, colour, on-sky
location and other more subtle variables (Lindegren et al., 2021b), and the
formal parallax uncertainties from Gaia EDR3 are believed to be underestimated
particularly at the bright end by a few $10$s of percent (El-Badry et al.,
2021; Maíz Apellániz, 2022). However, the existence of period–luminosity
relations for certain stellar types opens up the possibility of measuring
these systematic effects (e.g. Ren et al., 2021) and indeed a fully
probabilistic model can simultaneously calibrate the properties of standard
candles _and_ measure systematic issues with the data (e.g. Sesar et al.,
2017; Chan & Bovy, 2020).

In this paper, new period–luminosity relations for O-rich Mira variables in
the Milky Way are provided using data from Gaia Data Release 3. The relations
are derived using a probabilistic model incorporating distance priors and a
model for Gaia parallax systematics. The new relations are then used to
estimate $H_{0}$ using Mira variables in the Type Ia supernova host galaxy,
NGC 1559. Section 2 describes the dataset employed in this work to measure the
period–luminosity relations before the methodology is described in Section 4.
The new O-rich Mira variable period–luminosity calibrations are presented and
discussed in Section 5 before they are utilised for the estimation of the
Hubble constant in Section 6, taking into account the non-negligible C-rich
contamination. The conclusions are presented in Section 7. In three
appendices, the approximate completeness of the Gaia DR3 Mira variable
catalogue is discussed (Appendix A), the expected Gaia performance for
pulsating AGB stars is presented (Appendix B) and the period–luminosity
relations for the LMC, SMC and the Sgr dwarf spheroidal galaxy are estimated
(Appendix C).

##  2 O-rich Mira variables in Gaia DR3

Figure 1:  Period–magnitude diagram for the OGLE long-period variable sample
in the LMC (Soszyński et al., 2009). The left panel shows the logarithmically
coloured density of the full sample and the central panel highlights those
stars with $\Delta G>\Delta G_{\mathrm{thresh}}=0.865\,\mathrm{mag}$. The
different sequences from Wood et al. (1999), Wood (2000) and Ita et al. (2004)
are marked as dashed lines and labelled (note ‘C’ sequence should not be
confused with C-rich). In the central panel, a fraction
$f_{\mathrm{contam.}}=0.043$ of the selected sources fall off the C sequence
(defined by the grey shaded area). Restricting further to those with $\Delta
G_{\mathrm{Fourier}}>0.865$ produces a contamination fraction of
$f^{\mathrm{amp}}_{\mathrm{contam.}}=0.013$. The right panel shows the (median
and $\pm 1\sigma$) ratio of measured to true $\Delta G$ for a set of simulated
sinusoidal light curves with periods assigned from the dataset used in this
paper and randomly sampled phases sampled using the EDR3 photometric scanning
law. Results for both the DR2 and DR3 sampling period are shown. $\Delta G$ is
predominantly biased low, particularly at the aliasing period of
$190\,\mathrm{day}$. This is illustrated in the inset for the simulated
sampled light curve of a $183\,\mathrm{day}$ source (with the solid dots the
DR3 measurements and circles the DR2 measurements).  Figure 2: Properties of
the O-rich Mira sample: the top left panel shows the distribution of the
Wesenheit index difference from Lebzelter et al. (2018) for all stars
classified here as Mira variables split by O-rich or C-rich classification
(the thin dotted histograms show the distribution of those sources cross-
matched to the Suh & Hong (2017) catalogue using their classifications). The
lower left panel shows the distribution of the Wesenheit index difference vs.
period with the blue shading showing the logarithmic density of the O-rich
Mira variables and the orange contour containing $90\,\mathrm{per\,cent}$ of
C-rich Mira variables. The top left panel shows the period distribution of the
O-rich Mira variables (the small depletion around $\sim 190$ days is due to
Gaia’s scanning strategy). The lower right panel shows the top-down
Galactocentric view (from the North Galactic Pole) of the O-rich Mira
variables using the LMC period–luminosity relation. In the analysis, stars in
the midplane ($|b|<3\,\mathrm{deg}$) and those in the bulge region (shown in
projection by the grey wedge) are also removed.

The primary data source is the long period variable (LPV) candidate catalogue
(Lebzelter et al., 2022) from Gaia DR3 (Gaia Collaboration et al., 2016,
2022). Gaia DR3 includes $34$ months of data with a mean number of
observations per source of $43$. The Gaia variability processing consists of
two stages: an initial classification of all likely variable sources (Holl et
al., 2018; Rimoldini et al., 2019, 2022) and then a series of specific object
studies (SOS) that further process each variability class. The initial
classification was performed on all sources with at least $5$ Gaia field-of-
view transits in their processed and cleaned photometric time series, and that
were classified as likely variable when comparing to the variability of
literature variable objects and the $75\%$ least variable Gaia sources at each
magnitude. Classification into separate variability classes was then performed
using features including time-series summary statistics, Lomb-Scargle periods,
colours and parallax, and a training set composed of literature
classifications. Gaia DR3 published all stars classified as LPV with $G$
$5$th-$95$th percentile greater than $0.1\,\mathrm{mag}$,
$G_{\mathrm{BP}}-G_{\mathrm{RP}}>0.5$, at least $10$ visibility_periods_used,
a reported renormalized unit weight error (RUWE), more than $9$ $G$
observations and a ratio of number of $G_{\mathrm{RP}}$ measurements in the
cleaned time series to number of $G$ measurements $>0.5$. Of these a stricter
subset with more than $12$ $G$ observations and number of $G_{\mathrm{RP}}$
measurements to number of $G$ measurements ratio of $>0.8$ were considered in
the SOS (along with $522$ sources that satisfy all SOS LPV requirements but
were mostly classified as symbiotic stars). Periods were found using a
generalised Lomb-Scargle method and were published if the period was
$>35\,\mathrm{day}$ and shorter than the time series duration, the $G$ band
signal-to-noise was greater than $15$ and no strong correlation was detected
between the photometric time-series and the image parameter determination time
series. This resulted in $392\,240$ LPV candidates with published periods from
$2\,326\,297$ sources in Gaia DR3 classified as LPV. The completeness of the
full LPV candidates catalogue and the subset with published periods is briefly
assessed in Appendix A. In conclusion, the completeness of the Milky Way
sample with periods is $\gtrsim 90\,\mathrm{per\,cent}$ for
$|b|>3\,\mathrm{deg}$ and $\Delta G>\Delta G_{\mathrm{thresh}}$ with respect
to the full Gaia DR3 source catalogue.

Due to Gaia’s scanning strategy, periods around $\sim 190$ days and below
$120$ days are susceptible to aliasing. Cross-matching those LPVs later
defined as Mira variables with the AAVSO International Variable Star Index
(VSX, Watson et al., 2006, downloaded 30th April 2022), ASAS-SN (Jayasinghe et
al., 2018, 2019a, 2019b), and the OGLE LPV sample (Soszyński et al., 2009;
Iwanek et al., 2022), the fraction of likely aliases (periods disagreeing by
more than $25\,\mathrm{per\,cent}$ – the approximate width of the one-to-one
relations upon cross-matching) is $3.4\,\mathrm{per\,cent}$,
$3.8\,\mathrm{per\,cent}$ and $0.8\,\mathrm{per\,cent}$ respectively,
indicating aliasing is a minor issue and largely the Gaia periods are accurate
(Lebzelter et al., 2022). Mowlavi et al. (2018) report that the Gaia DR2 LPV
catalogue is contaminated at the few percent level by young stellar objects
(YSO) and the same is expected for Gaia DR3. With some parallax information,
these can be identified as intrinsically fainter than the Mira variables. A
conservative cut is employed by removing a handful of sources with
$G-5\log_{10}(100\,\mathrm{mas}/(\varpi-3\sigma_{\varpi}))>1.75(G_{\mathrm{BP}}-G_{\mathrm{RP}})-3$.
After this cut, the vast majority of the sample cross-matched to VSX and ASAS-
SN are classified by these collections as LPVs (Mira variables, semi-regular
variables or otherwise) with the largest contaminant being YSOs but only at
the $\lesssim 0.2\,\mathrm{per\,cent}$ level.

The Gaia DR3 catalogue of candidate LPVs is complemented with variables from
VSX (Watson et al., 2006, downloaded 30th April 2022). VSX is a compilation of
variables initially built from the General Catalogue of Variable Stars (Samus’
et al., 2017). All sources labelled as type ‘M’ (Mira), ‘M:’ (ambiguous Mira),
‘SR’ (semi-regular), ‘SRA’ (semi-regular variables similar to Mira but with
small amplitude) and ‘LPV’ (long-period variable) are selected and cross-
matched to Gaia DR3, removing variables already in the Gaia DR3 LPV catalogue.
This adds $2587$ stars to our Mira variable sample and $632$ to our more
restricted sample used for fitting defined later.

Only LPVs with 2MASS photometry are used (cross-matched within
$1\,\mathrm{arcsec}$ using proper motions to account for the epoch
difference). 2MASS observations are single-epoch so will produce additional
scatter about any fitted period–luminosity relation. However, the relations
should be unbiased representations of the arithmetic mean magnitude
period–luminosity relations (as required in the later $H_{0}$ analysis). Mean
$J$, $H$ and $K_{s}$ magnitudes could be estimated using the Gaia light
curves. However, the epoch difference ($\sim 17$ year) is large enough that
the typical Gaia frequency uncertainties ($\Delta\nu\approx
0.05\,\mathrm{year}^{-1}$, $\Delta\nu/\nu\approx 4\,\mathrm{per\,cent}$)
produce $\Delta\phi/\phi\approx\Delta\nu\times(17\,\mathrm{year})\approx
85\,\mathrm{per\,cent}$ uncertainties in the phase at the 2MASS epoch. This
simple consideration does not account for uncertainties in the light curve
fits at fixed period, the uncertainty in the amplitude ratios between the
$JHK_{s}$ and $G$ bands, or any stochastic cycle-to-cycle variation that can
be observed in Mira variables (Ou & Ngeow, 2022; Iwanek et al., 2022).
Therefore, it appears with the current data any attempt to find the mean
magnitudes from the single-epoch data will only add noise. For this reason,
only the single-epoch measurements are used here.

###  2.1 Selecting Mira variables

To isolate a sample of Mira variables from the combined Gaia DR3 and VSX LPV
candidates catalogue, two amplitude measures are combined: $\Delta
G_{\mathrm{Fourier}}$, the amplitude derived from a Fourier fit provided in
the Gaia DR3 LPV candidate catalogue (the amplitude column gives the $G$-band
semi-amplitude from the Fourier fit i.e. half the required value) and $\Delta
G$, the $G$-band amplitude measure computed from the reported Gaia
uncertainties (Belokurov et al., 2017). This latter quantity is defined as

$\Delta G=\frac{5\sqrt{2}}{\ln
10}\frac{\sqrt{\texttt{phot\\_g\\_n\\_obs}}}{\texttt{phot\\_g\\_mean\\_flux\\_over\\_error}}.$
(1)

For light curves that are near sinusoidal and sampled fairly over period, this
measure will be equal to the Fourier amplitude. Although the two measures
correlate strongly with each other (Sanders & Matsunaga, 2023), both
quantities are used for defining Mira variables as they behave differently for
poorly sampled lightcurves. If the lightcurve is undersampled, $\Delta
G_{\mathrm{Fourier}}$ overestimates the amplitude as only a limited range of
phases are used in the fit. However, as $\Delta G$ is a measure of the data
scatter, it will be underestimated in this regime. In the right panel of Fig.
1 the ratio of the measured to true $G$-band amplitudes is shown for a set of
simulated sinusoidal light curves with periods assigned from the dataset used
in this work and randomly drawn phases. The light curves are sampled according
to the DR2 and DR3 photometric scanning laws (using the Gaia DR2 scanning law
from Boubert et al. 2021 and the DR3 nominal scanning law both as part of the
scanninglaw package, Green, 2018; Boubert et al., 2020; Everall et al., 2021,
with the data-taking gaps from Riello et al. 2021). Although $\Delta G$ can be
underestimated, particularly around the troublesome $190$ day period or for
stars that were part of the ecliptic pole scanning law which had approximately
two thirds of their observations taken within a month, it is rarely
significantly overestimated so when selecting using $\Delta G$ very low
contamination from lower amplitude non-Mira variables is expected.
Furthermore, when combined with $\Delta G_{\mathrm{Fourier}}$ it is quite
certain that the LPVs are high amplitude.

Following Grady et al. (2019), the cut $\Delta G>\Delta G_{\mathrm{thresh}}$
and $\Delta G_{\mathrm{Fourier}}>\Delta G_{\mathrm{thresh}}$ where $\Delta
G_{\mathrm{thresh}}=(5\sqrt{2}/\ln 10)10^{-0.55}\approx 0.865\,\mathrm{mag}$
is employed to isolate Mira variable stars. For the small set of stars from
VSX without counterparts in the Gaia DR3 LPV catalogue, $\Delta
G_{\mathrm{Fourier}}$ is not measured so only the $\Delta G>\Delta
G_{\mathrm{thresh}}$ cut is employed. In Fig. 1, the sample of OGLE LMC LPV
stars from Soszyński et al. (2009) is shown along with those that satisfy
$\Delta G>\Delta G_{\mathrm{thresh}}$. These selected stars predominantly lie
along the ‘C’ sequence associated with fundamental mode pulsation (Wood et
al., 1999; Wood, 2000; Ita et al., 2004) with only a fraction
$f_{\mathrm{contam.}}=0.04$ consistent with membership of a different
sequence. $21\,\mathrm{per\,cent}$ of the $\Delta G$–selected OGLE LPVs are
classified as semi-regular variables by Soszyński et al. (2013) on the basis
of their $I$ amplitudes but as acknowledged by these authors and Trabucchi et
al. (2021b) the traditional definitions of Mira variables are possibly not
appropriate as lower amplitude variables or irregular Mira variables follow
the same period–luminosity relation (as evident from Fig. 1) and are probably
governed by the same physics. If the set of OGLE LPVs with Gaia DR3 Fourier
amplitudes is considered, cutting on both $\Delta G>G_{\mathrm{thresh}}$ and
$\Delta G_{\mathrm{Fourier}}>G_{\mathrm{thresh}}$ reduces the contamination
from non-C-sequence stars to $f_{\mathrm{contam.}}=0.01$.

###  2.2 Separation of O-rich and C-rich Mira variables

LPVs exhibit oxygen-rich and carbon-rich chemistry depending on the initial
mass and metallicity of the star (Höfner & Olofsson, 2018). Of these two
populations, the O-rich subset are more useful as they follow a tighter
period–luminosity relation (Ita & Matsunaga, 2011). Although significant
within the LMC, C-rich Mira variables are rarer within the Galactic disc
(Blanco et al., 1984) and tend to be confined to the outer disc. C-rich Mira
variables are typically redder and dustier than their O-rich counterparts.
Lebzelter et al. (2018) showed that O-rich and C-rich Mira variables within
the LMC can be separated in the plane of $W_{\mathrm{RP,BP-RP}}-W_{Ks,J-Ks}$
vs. $K_{s}$. Here the two Wesenheit indices are $W_{\mathrm{RP,BP-
RP}}=G_{\mathrm{RP}}-1.3(G_{\mathrm{BP}}-G_{\mathrm{RP}})$ and
$W_{Ks,J-Ks}=K_{s}-0.686(J-K_{s})$. The boundary employed by Lebzelter et al.
(2018) is slightly curved in ‘colour’-magnitude space but the curvature is
weak and a pure $W_{\mathrm{RP,BP-RP}}-W_{Ks,J-Ks}$ cut performs similarly.

Lebzelter et al. (2022) have discussed how O-rich and C-rich LPVs can be
distinguished using the Gaia DR3 BP/RP spectra due to the distinct separation
of a set of bandheads arising from TiO for O-rich stars and CN for C-rich
stars. As acknowledged by Lebzelter et al. (2022), the bandhead separation
diagnostic performs poorly for very red sources leading to the
misclassification of many O-rich sources as C-rich. Sanders & Matsunaga (2023)
utilised an unsupervised classification approach using the BP/RP spectra that
uses the UMAP (Uniform Manifold Approximation and Projection, McInnes et al.,
2018) algorithm on the normalized coefficients. This approach performs better
than the published Gaia DR3 classifications for highly-extincted stars. For
those stars without BP/RP spectra, Sanders & Matsunaga (2023) used a
supervised classification algorithm (Chen & Guestrin, 2016, XGBoost) trained
on Gaia and 2MASS photometric data, periods and amplitudes for the stars with
unsupervised classifications. This produces a $95\,\mathrm{per\,cent}$ purity
C-rich sample and $99.5\,\mathrm{per\,cent}$ purity O-rich sample (due to the
dominance of O-rich sources in the sample). Here the BP/RP unsupervised
classifications are used when available, falling back to the supervised
photometric classifications when no BP/RP spectra is provided in Gaia DR3.
Fig. 2 shows the distribution of the Mira variable sample in the Wesenheit
‘colour’ vs. period where the separation of the O-rich and C-rich populations
is clear. A simple cut of $W_{\mathrm{RP,BP-RP}}-W_{Ks,J-Ks}<1$ would remove
most C-rich sources but would also remove some longer period O-rich sources.
The stars in the final sample that are also in the catalogue of Suh & Hong
(2017) are shown by the dotted histogram separated using these authors’
classification. The classifications are a combination of low-resolution
spectroscopic, maser and photometric classifications. Using our
classifications to isolate O-rich stars results in only $5$ of the $783$
matches ($0.6\,\mathrm{per\,cent}$) with the Suh & Hong (2017) catalogue being
classified by them as C-rich (using the updated IRAS PSC catalogue of Suh
(2021) results in $12$ of $867$ matches classified as C-rich,
$1.4\,\mathrm{per\,cent}$).

###  2.3 Summary of selections

In summary, the Gaia DR3 LPV candidates with reported periods have been
combined with additional LPVs from VSX. Mira variables have been isolated by
cutting on the $G$-band Fourier amplitude, $\Delta G_{\mathrm{Fourier}}$, and
$G$-band scatter, $\Delta G$, and potential YSO contaminants have been removed
with a parallax cut. O-rich and C-rich separation has been performed using the
BP/RP spectra where available, and otherwise using broadband Gaia and 2MASS
photometry combined with periods and amplitudes. Considering the issues of
period aliasing, YSO contamination, non-Mira LPV contamination and C-rich
contamination altogether, it seems the cuts defined here produce a O-rich Mira
variable catalogue with a reliability upwards of $95\,\mathrm{per\,cent}$.

For fitting the period–luminosity relations, only stars with $G<17$,
$G_{\mathrm{BP}}-G_{\mathrm{RP}}>1.9$, distances $<25\,\mathrm{kpc}$ (as
estimated a priori using the LMC period–luminosity relations in Appendix C),
periods more than $100$ days and less than $1000$ days, period uncertainties
$<50\,\mathrm{per\,cent}$ (the median period uncertainty is
$5\,\,\mathrm{per\,cent}$ and the 95th percentile is $11\,\mathrm{per\,cent}$)
and Gaia EDR3 RUWE $<1.4$ (see next section) are retained. Furthermore, stars
in the bulge region ($|\ell|<15,|b|<10$), those in the Galactic mid-plane
($|b|<3\,\mathrm{deg}$), those within $15$ deg of the LMC, those within $10$
deg of the SMC and those at distances greater than $18\,\mathrm{kpc}$ within
$15$ deg of the Sgr dSph are removed. These on-sky selections are visualized
in Fig. 12. With this set of cuts, there remain $15159$ O-rich and $875$
C-rich Mira variables (from an initial catalogue of $86477$ stars with the
$|b|$ and bulge cuts most severely reducing the sample). The lower right panel
of Fig. 2 shows the view of the sample from the Galactic North Pole using the
LMC period–luminosity relation.

##  3 Astrometric data quality

###  3.1 Initial considerations

Figure 3: Column-normalized distributions of the reduced unit weight error
(RUWE, top) and the astrometric excess noise (bottom) against various
quantities for our Milky Way Mira variable sample. The two horizontal lines in
the top panels show the RUWE cuts ($1.4$ and $2$) employed in this work. The
blue points in the second lower panel shows the binned distribution (median
$\pm 1\sigma$) of an M dwarf sample defined in the text. The dashed line in
the second lower panel is the median trend subtracted in each of the other
lower panels to produce the black line.

To confidently use the Gaia EDR3 astrometric data for period–luminosity
calibration, their quality must be assessed (note Gaia DR3 did not update the
astrometry so EDR3 and DR3 astrometry refer to the same thing). This is a
particular concern for Mira variables as they are some of the reddest sources
observed by Gaia. Additionally, their variability (in both colour and
magnitude) makes the astrometry challenging, and as discussed in Mowlavi et
al. (2018) in the current Gaia data releases epoch photometry is not utilised
in the astrometric solution (Lindegren et al., 2021a) which could lead to
errors for variable sources (Pourbaix et al., 2003, see Appendix B) . There
are a number of recommended quality cuts for handling Gaia data (Fabricius et
al., 2021) but the only quality criterion used here is the renormalized unit-
weight error (RUWE) from Gaia EDR3 by ensuring all stars have RUWE $<1.4$ (a
test with $<2$ is also run). Although nearly all of the sample has significant
($>3$) astrometric excess noise, Lindegren et al. (2021a) caution against
using astrometric excess noise for very red sources
($G_{\mathrm{BP}}-G_{\mathrm{RP}}>3$) as it likely reflects shortcomings of
the instrument and attitude modelling. Also, a large fraction of the sample
have ipd_gof_harmonic_amplitude $>0.2$ ($62\,\mathrm{per\,cent}$) and
ipd_frac_multi_peak $>2$ ($33\,\mathrm{per\,cent}$) which is indicative of
poor LSF/PSF fits due to possible binarity (Lindegren et al., 2021a). However,
the LSF/PSF calibrations (Rowell et al., 2021) have only been performed down
to $\nu_{\mathrm{eff}}=1.24\,\mu\mathrm{m}^{-1}$ so it is anticipated that
redder sources will not have well fitting LSF/PSFs. Furthermore, for six-
parameter solutions a default LSF/PSF at
$\nu_{\mathrm{eff}}=1.43\,\mu\mathrm{m}^{-1}$ is utilised which perhaps makes
the IPD statistics unreliable for the significantly redder sources. Finally,
these sources typically fall outside the advised adjusted BP/RP excess factor
range as a function of magnitude but this is probably due to their variability
(as already highlighted by figure 21 of Riello et al., 2021).

In Fig. 3 the column-normalized distributions of RUWE and astrometric excess
noise are shown against various other quantities for the Mira variable sample.
The RUWE distributions are largely flat with all plotted quantities except for
an enhancement in the Galactic midplane and a small uptick at bluer
$(G_{\mathrm{BP}}-G_{\mathrm{RP}})$. There is some slight evidence of an
increase in RUWE for nearby, brighter sources. The astrometric excess noise
shows strong trends, particularly with colour. However, the astrometric excess
noise vs. $(G_{\mathrm{BP}}-G_{\mathrm{RP}})$ is shown for a sample of M dwarf
stars from Gaia DR3 with $(G_{\mathrm{BP}}-G_{\mathrm{RP}})>3.5$, RUWE $<1.4$
and $\varpi>10\,\mathrm{mas}$. This traces the trend in the Mira variables in
the overlapping region. It therefore appears the astrometric excess noise
trend arises from poor characterisation of the instrument performance, rather
than anything intrinsic. In other panels the trends can be related to the
fundamental trend in colour i.e. redder stars are fainter, typically higher
amplitude and found more often in the midplane. This is corroborated by the
black lines which depict the median trends after subtracting the median colour
dependence (shown as a black dashed line in the second lower panel).

The use of mean photometry in the astrometric solution leads to two effects:
(i) the centroids have a residual uncorrected offset due to the colour
variation of the sources and (ii) an average astrometric error instead of the
epoch astrometric errors is used. Using the pseudocolour uncertainties for the
six-parameter solutions, the typical centroid shift with effective wavenumber
is estimated as $\sim 2\,\mathrm{mas}\,\mu\mathrm{m}$ (see also de Bruijne et
al., 2006; Lindegren et al., 2021a) which using the typical amplitudes and
colours of the Mira variable sample is $\sim 6\,\mathrm{per\,cent}$ of the
reported uncertainties in the median. The variation of the epoch uncertainties
due to the typical $G$ and $(G_{\mathrm{BP}}-G_{\mathrm{RP}})$ amplitudes of
the sample is $\sim 20\,\mathrm{per\,cent}$. The full analysis presented in
Appendix B demonstrates that in combination these effects lead to a modest
underestimate of the astrometric uncertainties of at most
$10\,\mathrm{per\,cent}$ with the largest underestimates arising from the
highest amplitude stars.

###  3.2 Intrinsic photocentre wobble

A further concern is that AGB stars have large radii, $\sim 1\,\mathrm{AU}$,
and complex surface dynamics and, as highlighted recently by Chiavassa et al.
(2018) (see also van Belle et al., 2002), the motion of the atmosphere can
lead to shifts of the photocentre typically of order
$5-10\,\mathrm{per\,cent}$ of the radius. Appendix B investigates this issue
in considerable detail and here only simple arguments as to its impact on the
Gaia EDR3 astrometry are presented. The previous comparison with the M-dwarf
sample suggests the quality of the astrometry for the Mira variables arises
from Gaia’s performance rather than any intrinsic noise, but this is validated
further here.

Chiavassa et al. (2011) presented a simulation of the red supergiant
Betelgeuse finding a $G$-band photocentre wobble of $0.065\,\mathrm{AU}$,
about $2\,\mathrm{per\,cent}$ of its radius, whilst Chiavassa et al. (2018)
presented $8$ simulations of AGB stars with typical $G$-band photocentre
wobbles of $5-10\,\mathrm{per\,cent}$ with longer period stars (or more
precisely longer pressure scaleheight) having a larger wobble. In the optical,
the photocentre wobble is composed of long variations on the order of years
due to large convective cells covering of order $\sim 1/3$ the stellar radius
(more evident in infrared observations) with shorter variations on the order
of months due to smaller convective cells in the upper atmospheres of size
$10\,\mathrm{per\,cent}$ the stellar radius. For nearby AGB stars, this
photocentre wobble can be a significant observable effect presenting a
fundamental error floor for the astrometry. However, due to the stochasticity
of the AGB photocentre wobble and the lack of preferred direction relative to
the parallax ellipse and proper motion vector, it is expected that over long
enough timespans (or averaged over many stars) the wobble should manifest as
an additional random uncertainty and the astrometric parameters will be
unbiased but possibly with poorly estimated uncertainties (Chiavassa et al.,
2011).

Figure 4:  Approximate single-epoch astrometric uncertainty from Gaia for the
O-rich Mira variable sample (in AU). The along-scan astrometric uncertainty,
$\sigma_{\mathrm{AL}}$, is approximated from the Gaia parallax uncertainty,
$\sigma_{\varpi}$, as $0.53\sqrt{N}\sigma_{\varpi}$. There are $18$
observations per transit ($9$ CCD observations in each field-of-view) such
that, ignoring systematics, the single-epoch astrometric uncertainty is
$\sigma_{\mathrm{AL}}/\sqrt{18}$. This is transformed into AU using the
parallax computed from the LMC period–luminosity relation. The expected
photocentre wobble for ‘normal’ and dusty Mira variables is shown by the
orange solid and pink dashed lines assuming the period–radius relations from
van Belle et al. (2002, vB+2002) and a $3\,\mathrm{per\,cent}$ radial
variation (the errorbars give the amplitude of the uncertainty in the
relations and the thicker parts of the lines are the regions over which
vB+2002 had data). The AGB models from Chiavassa et al. (2018) are shown as
red points. The measured photocentre wobble for o Cet (Mira, large hexagon)
and $\alpha$ Her and $\alpha$ Ori (small faint hexagons) are shown (note the
latter two stars are red supergiants and their periods have been used for
convenience to place them in the plot). Essentially all of the sample lies
above the models suggesting the astrometric uncertainties are not dominated by
photocentre wobble and the parallaxes are reliable.

Photocentre wobble is only detectable when it is similar to or greater than
the Gaia single-epoch astrometric uncertainty. Lindegren et al. (2021a)
provides the median along-scan astrometric uncertainty in Gaia EDR3,
$\sigma_{\mathrm{AL}}$, as a function of $G$ but with no information on the
uncertainty as a function of colour. Belokurov et al. (2020) have demonstrated
that $\sigma_{\mathrm{AL}}$ is approximately related to the reported parallax
uncertainty, $\sigma_{\varpi}$, as $\sigma_{\mathrm{AL}}\approx
0.53\sqrt{N}\sigma_{\varpi}$ where $N$ is the number of observations
(astrometric_n_good_obs_al) allowing for the estimation of
$\sigma_{\mathrm{AL}}$ for a range of different magnitudes, colours, on-sky
positions etc. Gaia typically makes $18$ astrometric observations in a short
timespan ($9$ CCDs for each of the two fields of view) such that in the
absence of systematic uncertainties, the single-epoch along-scan astrometric
uncertainty is $\sim\sigma_{\mathrm{AL}}/\sqrt{18}$. It is this uncertainty
that must be compared with the expected AGB photocentric wobble. Using the
parallax to transform this astrometric uncertainty into a physical scale gives
$\sim\sqrt{N/350}(\sigma_{\varpi}/0.1\mathrm{mas})(\mathrm{mas}/\varpi)(0.24\mathrm{AU})$
where typical values for $N$ and $\sigma_{\varpi}$ for the sample are used.
For a star at $1\,\mathrm{kpc}$, the single-epoch astrometric uncertainty is
larger than the expected $5-10\,\mathrm{per\,cent}$ of the radius wobble
(assuming the radius is $1\,\mathrm{AU}$). In Fig. 4 the single-epoch
astrometric uncertainty, $(1/\sqrt{18})0.53\sqrt{N}\sigma_{\varpi}/\varpi$ in
AU, is displayed for the full O-rich Mira variable sample with RUWE$<1.4$ and
$G<17$ using $\varpi$ estimated from the LMC period–luminosity relation
(Appendix C). This can be compared to the AGB models from Chiavassa et al.
(2018), the measured photocentre wobble from o Cet (Mira) and the two
supergiants, $\alpha$ Her and $\alpha$ Ori, (Chiavassa et al., 2011) and a
simple model of the photocentre wobble using the period-radius models from van
Belle et al. (2002) and a $3\,\mathrm{per\,cent}$ radius wobble that fits the
Mira observation well.

We see that because the bulk of the sample has $\varpi\ll 1\,\mathrm{mas}$,
the physical scale Gaia is capable of probing for the sample is significantly
greater than $0.1\,\mathrm{AU}$ and the uncertainty budget is dominated by
Gaia’s limitations. If the astrometric excess noise is instead used as the
measure of along-scan astrometric uncertainty a similar result is found. This
gives confidence that for the majority of the considered sample the astrometry
should be free from any effects arising from intrinsic photocentre wobble and
that Gaia is capable of providing precision measurements for this type of
star. However, this may be more of a concern with future data releases with
improved astrometric uncertainties for red stars. For example, Chiavassa et
al. (2011) estimated that the photocentre wobble should be a measurable effect
from the Gaia uncertainties for stars within $4.4\,\mathrm{kpc}$ assuming the
predicted wobble from models of Betelgeuse. Their analysis assumed final Gaia
parallax uncertainties of $7.8\,\mu\mathrm{as}$ whilst the typical uncertainty
for the present sample is an order of magnitude larger around
$0.1\,\mathrm{mas}$. However, it should be stressed that improved measurements
over longer baselines will likely not produce on average biased astrometric
results, but more affect the reported uncertainties.

One caveat here is that the AGB model expectation might be very wrong and
$\sigma_{\mathrm{AL}}$ in fact does reflect the photocentre wobble rather than
any limitation on Gaia’s performance. This is unlikely considering for the
bulk of stars the single-epoch astrometric uncertainty is of order the radius
of the star and also that the measured photocentre wobble of Mira suggests if
anything the AGB models of Chiavassa et al. (2018) produce too large a
photocentre wobble. Furthermore, a comparison with M dwarf stars at similar
colours and magnitudes shows similar $\sigma_{\mathrm{AL}}$ and astrometric
excess noise to the sample used here (see Fig. 3 and Fig. 13). No photocentre
wobble is expected for these sources suggesting in the majority of cases
$\sigma_{\mathrm{AL}}$ is governed by Gaia’s limitations.

A much fuller analysis of the expected Gaia performance for AGB stars is
presented in Appendix B and reaches the same conclusions as the simpler
considerations presented here. The analysis demonstrates that on average the
astrometric parameters for the sample of stars used in this work are unbiased
but the uncertainties are underestimated for $G\lesssim 11$ and
$\varpi>0.5\,\mathrm{mas}$ (a small fraction of the total sample).

##  4 Period–luminosity relation for O-rich Mira variables

A probabilistic model is introduced to measure the period–luminosity relation
for the sample of O-rich Mira variable stars from Gaia DR3. This allows the
inclusion of uncertainties in the data and a prior when transforming from the
uncertain parallax measurements to absolute magnitudes (Bailer-Jones et al.,
2018; Luri et al., 2018). Furthermore, the fact the Mira variables appear to
follow a period–luminosity relation can be used to simultaneously calibrate
this relation and measure the parallax zeropoint and parallax uncertainties of
the sample (e.g. Sesar et al., 2017; Chan & Bovy, 2020). As highlighted
previously, variations in the accuracy of the astrometry with both colour and
magnitude are anticipated. This is particularly important for the Mira
variables as they are some of the reddest sources observed by Gaia and many
fall outside the effective wavenumber range covered by previously-published
zeropoint corrections (Lindegren et al., 2021b).

As a preliminary illustration of the sample and an indication of the ability
to measure the period–luminosity relation accurately, the $K_{s}$ absolute
magnitude computed from the Gaia EDR3 parallax vs. period is shown in Fig. 5.
The $K_{s}$ magnitudes have been corrected for extinction as described later
in Section 4.2 and the Gaia EDR3 parallaxes have been zeropoint corrected
using the Lindegren et al. (2021b) corrections evaluated at
$\nu_{\mathrm{eff}}=1.25\,\mu\mathrm{m}^{-1}$ as described later in Section
4.4. For comparison, the period–luminosity relation for the LMC as derived in
Appendix C is shown. The subset of stars with parallax uncertainties better
than $10\,\mathrm{per\,cent}$ align nicely with the LMC relation, possibly
falling slightly under in the mean, although this effect is partly due to the
selection on parallax errors biasing the measurements towards higher
parallaxes and hence higher absolute magnitudes. The fuller sample shows a
significant scatter about the expected period–luminosity relation due to the
parallax uncertainties. In the following sections, the model for the data is
introduced, before the handling of the parallax zeropoint modelling is
described in more detail.

Figure 5: Period–luminosity distribution for the Milky Way O-rich Mira sample.
The background shows a log-scaled histogram of the subsample with parallax
signal-to-noise greater than $1$ (the large scatter is due primarily to the
parallax uncertainties). Absolute magnitudes have been computed using the Gaia
EDR3 parallaxes corrected by the Lindegren et al. (2021b) zeropoint
corrections evaluated at $\nu_{\mathrm{eff}}=1.25\,\mu\mathrm{m}^{-1}$ as
described in Section 4.4 and using the extinction corrections described in
Section 4.2. The orange points are the subset with parallax errors smaller
than $10\,\mathrm{per\,cent}$ (note this selection has the effect of biasing
the measurements towards higher parallaxes and absolute magnitudes). The
errorbar shows the minimum formal uncertainty from the parallax measurements
alone. The pink line shows the LMC O-rich Mira relation as derived in Appendix
C along with its scatter (a combination of that due to single epoch
observations and any intrinsic scatter due to population variations). The
black line shows similar for the MW relation derived in this work.

###  4.1 Probabilistic model

The joint single-star likelihood of the Gaia EDR3 parallax $\varpi$ and a
magnitude $m$ given the $G$ magnitude, effective wavenumber
$\nu_{\mathrm{eff}}$ (pseudo-colour for six-parameter astrometric solutions),
period $P$ and on-sky location $(\ell,b)$ (and corresponding uncertainties) is
expressed as

$\begin{split}p(\varpi,m|G,\nu_{\mathrm{eff}},P,\ell,b)=\int\mathrm{d}s\,p(\varpi|&s,G,\nu_{\mathrm{eff}},\ell,b)\\\
\times&p(m|s,P)p(s|\ell,b),\end{split}$ (2)

where

$\begin{split}p(\varpi|s,G,c,\ell,b)&=\mathcal{N}(\varpi|1/s+\varpi_{0},f^{2}_{\varpi}\sigma_{\varpi}^{2}+\sigma^{2}_{\varpi,0}),\\\
p(m|s,P)&=\sum_{j=1}^{j=2}\vartheta_{j}\mathcal{N}(m|m_{\mathrm{abs}}(P)+\mu,\sigma_{m,j}^{2}(P)),\\\
p(s|\ell,b)&=\frac{s^{2}}{2L^{3}(\ell,b)}\mathrm{e}^{-s/L(\ell,b)}.\end{split}$
(3)

$\mathcal{N}(x|\mu,\sigma^{2})$ is a normal distribution with mean $\mu$ and
standard deviation $\sigma$. Here $s$ is the true distance (with corresponding
distance modulus $\mu$), $\varpi_{0}(G,\nu_{\mathrm{eff}},\ell,b)$ is a
colour-, magnitude- and spatially-dependent parallax zeropoint offset
(described in a later subsection), $\sigma_{\varpi}$ are the reported parallax
uncertainties with $f_{\varpi}(G,\nu_{\mathrm{eff}})$ a colour- and magnitude-
dependent scaling (again specified later) and $\sigma_{\varpi,0}$ an
additional systematic error floor. A two-component Gaussian mixture model is
employed for the magnitude distribution about the predicted magnitude
$m_{\mathrm{abs}}(P)+\mu$ where $m_{\mathrm{abs}}(P)$ is the period–luminosity
relation. This mixture model accounts for possible outliers using a mixing
simplex $\vartheta_{j}$ ($\vartheta_{1}+\vartheta_{2}=1$). In Section 2 and
Fig. 2 the contamination was estimated to be at the few per cent level.
$p(s|\ell,b)$ is the distance prior. The various modelling choices are
discussed in the following subsections.

####  4.1.1 Period–magnitude relation

The adopted period-magnitude relation $m_{\mathrm{abs}}(P)$ for magnitude $m$
is given by

$m_{\mathrm{abs}}(P)=a_{m}+\begin{cases}b_{m}(\log_{10}P-2.3),&\mathrm{if}\log_{10}P\leq
2.6,\\\ 0.3b_{m}+c_{m}(\log_{10}P-2.6),&\mathrm{otherwise}.\end{cases}$ (4)

Period–luminosity relations for O-rich Mira variables have been computed using
those stars in the LMC (e.g. Ita & Matsunaga, 2011; Yuan et al., 2017a, b).
Typically a linear relation is appropriate for $P<400\,\mathrm{days}$ beyond
which a break occurs and the period–luminosity relation is steeper (Ita &
Matsunaga, 2011; Bhardwaj et al., 2019). This is often attributed to
additional luminosity arising from the onset of hot-bottom burning for stars
with $P>400\,\mathrm{days}$ (Whitelock et al., 2003). Following Ita &
Matsunaga (2011), a break is placed at $\log_{10}P=2.6$, which is validated by
fits to the LMC (see Appendix C) although Bhardwaj et al. (2019) advocate for
a slightly lower break at $300$ days. Often the entire period range is
modelled with a quadratic relation (Yuan et al., 2017b). Quadratic relations
are weakly disfavoured over broken linear relations for the LMC data (see
Appendix C) and also have the tendency to bias the relation for short periods
when attempting to fit the curvature at long periods. Furthermore, the broken
linear relation are made continuous (c.f. Ita & Matsunaga, 2011) as this form
is perhaps more physically motivated and reduces the number of parameters by
one.

In Appendix C period–luminosity relations for O-rich and C-rich Mira variables
in the LMC are provided using the form of the period-magnitude relation in
equation (4). The resulting $(b_{m},c_{m})$ posterior distributions are used
as priors for the Milky Way O-rich sample.

####  4.1.2 Period–amplitude relation

The scatter about the period-magnitude relation for each component is given by

$\sigma_{m,j}^{2}(P)=\sigma^{2}_{\mu}(P)+\mathrm{Var}(m_{\mathrm{abs,LMC}}(P))+\sigma_{m,\mathrm{obs}}^{2}+\sigma^{2}_{\mu,0,j}.$
(5)

The scatter consists of four terms: the first term $\sigma^{2}_{\mu}(P)$ gives
the intrinsic scatter about the period–amplitude relation. The bulk of the
spread arises from using single-epoch observations. Longer period variables
have larger amplitudes so a model of the form

$\sigma_{\mu}(P)=\sigma_{2.3}+\begin{cases}m_{\sigma-}(\log_{10}P-2.3),&\mathrm{if}\log_{10}P\leq
2.6,\\\
0.3m_{\sigma-}+m_{\sigma+}(\log_{10}P-2.6),&\mathrm{otherwise},\end{cases}$
(6)

is employed. The choice here mirrors the period–magnitude relation of equation
(4) as the break in the period–luminosity relation potentially due to hot-
bottom burning is accompanied by a break in the period–amplitude relation
(e.g. Matsunaga et al., 2009). Even if multi-epoch data from which accurate
mean magnitudes could be estimated were available, some intrinsic scatter
might be expected due to other hidden dependencies (e.g. age and metallicity,
Qin et al., 2018) so $\sigma_{\mu}(P)$ is considered as a quadrature sum of
the single-epoch scatter and intrinsic scatter. Again the posterior
distributions for fits to the (single-epoch) LMC data (Appendix C) are used as
priors for $F\equiv(\sigma_{2.3},m_{\sigma-},m_{\sigma+})$.

The second term in the scatter is $\mathrm{Var}(m_{\mathrm{abs,LMC}}(P))$
which gives the variance arising from the uncertainty in the period. For
simplicity, the additional spread in the magnitude
$\mathrm{Var}(m_{\mathrm{abs,LMC}}(P))$ is then computed using the fitted
period-magnitude relations for the LMC (see Appendix C). For large period
uncertainties, the prior understanding of the width of the period distribution
is also important. The Gaussian in $\log_{10}P$ with mean $q$ and width
$\sigma_{q}$ fitted to the LMC data in Appendix C is used. The uncertainty in
the period is then computed by combining with the prior distribution. The
third term in equation (5), $\sigma_{m,\mathrm{obs}}^{2}$, is the variance
arising from the photometric uncertainties, uncertainties in the extinction
and uncertainties in the extinction coefficients. The final term
$\sigma^{2}_{\mu,0,j}$ is an additional residual only employed for the outlier
component such that $\sigma^{2}_{\mu,0,1}=0$.

###  4.2 Extinction corrections

The magnitudes $m$ must be corrected for the effects of extinction. When
available, the Green et al. (2019, Bayestar 2019) extinction estimates and
their uncertainties are evaluated at the distance of each Mira variable using
the LMC Wesenheit period–luminosity relations of Appendix C. The reported
extinctions are assumed to be exactly equal to $E(B-V)$ on the Schlegel et al.
(1998) scale (validated as $E(B-V)=1.02E(g_{\mathrm{PS}}-r_{\mathrm{PS}})$
from Wang & Chen 2019) so must be adjusted to account for the
$14\,\,\mathrm{per\,cent}$ reduction reported by Schlafly & Finkbeiner (2011)
and to convert to ‘true’ $E(B-V)$. The extinction estimates are flagged as
possibly unreliable if stars are beyond the faintest main-sequence star in
Pan-STARRS at a given on-sky location. Green et al. (2019) also use giant
stars in their extinction estimates such that the extinctions beyond the
faintest main-sequence star can be constrained. However, the giant models are
less certain than the main sequence models. To account for this, the
extinction uncertainties are arbitrarily inflated by a factor two for the
estimates flagged as unreliable. Outside the Pan-STARRS footprint, the
extinction map from Schlegel et al. (1998) is used accounting for the
recalibration from Schlafly & Finkbeiner (2011) and a $16\,\mathrm{per\,cent}$
uncertainty is employed.

For computing the extinction in a general band, the extinction coefficients
from Wang & Chen (2019) are used and their provided uncertainties in the
extinction coefficients are propagated. An alternative to explicitly
correcting for extinction is to use the Wesenheit magnitudes given by
$m=W_{x,y-x}\equiv x-e(y-x)$ where the extinction coefficient $e\equiv
A(x)/E(y-x)$ from Wang & Chen (2019) (or from Yuan et al., 2013, as a model
variant) and $x$ and $y$ are the observed extincted magnitudes. Variations in
$e$ are somewhat degenerate with changes to the period–luminosity relation so
it is preferable to keep $e$ fixed although it is allowed to vary in one model
variant.

###  4.3 Distance prior

$p(s|\ell,b)$ in equation (3) is the prior on distance. Luri et al. (2018)
emphasised the importance of using an appropriate prior when working with
parallax data. Whilst it is tempting to simultaneously constrain a Galactic
density model prior alongside calibrating the parallax data and
period–luminosity relation, this is non-trivial as the sample is subject to
complex selection effects (see Appendix A). For example, the effects of the
Gaia scanning law are visible on small scales. More severe, however, is the
incompleteness in the plane due to extinction. Astraatmadja & Bailer-Jones
(2016) explored using a fixed Galactic prior for finding distances from
parallaxes with and without photometric information, but find that the simple
exponentially decreasing space density prior
$p(\boldsymbol{x})\propto\exp(-s/L)$ produces similar (but in the case of the
Galactic centre regions significantly less biased) distance estimates and is
significantly simpler to work with. Bailer-Jones et al. (2018) used the
exponentially decreasing space density prior to estimate distances from the
full Gaia DR2 dataset, adopting a scalelength $L(\ell,b)$ that varies with on-
sky position. The adopted functional dependence is determined in on-sky bins
from a Gaia mock catalogue and fitted with a spherical harmonic series.
Bailer-Jones et al. (2021) updated this procedure for Gaia EDR3 by introducing
an additional parameter into the prior ($p(\boldsymbol{x})\propto
s^{\beta-2}\exp(-(s/L)^{\alpha})$ with $\alpha$, $\beta$ and $L$ all functions
of on-sky position. The simpler single-parameter exponentially decreasing
prior is chosen adopting a spherical harmonic series in $\ln L$ given by

$\ln
L/L_{0}=\sum_{n=1}^{n_{\mathrm{max}}}\sum_{m=0}^{n}\Big{[}s_{nm}P_{n}^{m}(\sin
b)\sin m\ell+c_{nm}P_{n}^{m}(\sin b)\cos m\ell\Big{]}.$ (7)

Here $P_{n}^{m}(x)$ are associated Legendre polynomials and $s_{n0}=0$. A QR
re-parametrization for this series is used which significantly improves
sampling111Stan Development Team. 2018. Stan Modeling Language Users Guide and
Reference Manual, Version 2.18.0. http://mc-stan.org.. The $P_{n}^{m}(\sin
b)\sin m\ell$ and $P_{n}^{m}(\sin b)\cos m\ell$ terms are combined into a
single matrix M of dimensions $(N_{\mathrm{data}},N_{\mathrm{series}})$ where
$N_{\mathrm{series}}=n_{\mathrm{max}}(n_{\mathrm{max}}+2)$, and the
coefficients $s_{nm}$ and $c_{nm}$ into a vector $\boldsymbol{S}$ of length
$N_{\mathrm{series}}$. $\textbf{{M}}=\textbf{{QR}}$ is decomposed into the
thin QR decomposition and then samples are taken in the transformed vector
$\boldsymbol{\tilde{S}}=\textbf{{R}}\boldsymbol{S}$. A shrinkage prior is
placed on $\boldsymbol{\tilde{S}}\sim\mathcal{N}(0,\tau)$ where $\tau$ follows
a unit half-Cauchy prior. The prior scalelength for datum $i$ is $\ln
L/L_{0}=(\textbf{{Q}}\boldsymbol{\tilde{S}})_{i}$. $n_{\mathrm{max}}$ is set
to $10$. The bar–bulge region is not used in the modelling to avoid biases
introduced by an inappropriate prior for this region.

###  4.4 Parallax zeropoint model

As reported initially by Lindegren et al. (2018) for Gaia DR2 and by Lindegren
et al. (2021b) and Fabricius et al. (2021) for Gaia EDR3, the reported Gaia
parallaxes and proper motions have zeropoint offsets and typically
underestimated uncertainties due to limitations in the instrument and attitude
modelling. Lindegren et al. (2021b) reported an approximation for the
zeropoint offset of the Gaia EDR3 parallaxes using samples of quasars,
binaries and stars in the LMC. The sources with five- and six-parameter
astrometric solutions were treated separately. The zeropoint correction was
approximated as a function of $G$ magnitude, ecliptic latitude and colour
(using $\nu_{\mathrm{eff}}$ for the five-parameter solutions and the pseudo-
colour for the six-parameter solutions). The implementation is available at
https://gitlab.com/icc-ub/public/gaiadr3_zeropoint. Several works (Riess et
al., 2021; Zinn, 2021; Huang et al., 2021) have validated the Lindegren et al.
(2021b) corrections, typically with some adjustment needed for bright stars
($G\lesssim 11$). Groenewegen (2021) presented an independent analysis of the
Gaia EDR3 parallax zeropoint using a sample of quasars and wide binaries. This
analysis differed from that presented by Lindegren et al. (2021b) by not
separating five- and six-parameter solutions, and using on-sky bins rather
than polynomials to capture the spatial dependence of the zeropoint. Maíz
Apellániz (2022) carried out a similar investigation of the Gaia EDR3
zeropoint to Lindegren et al. (2021b) using a sample of open clusters,
globular clusters and Magellanic Cloud data finding agreement with Lindegren
et al. (2021b) for faint objects ($G>13$) but some discrepancy for the
brighter objects.

In summary, these previous analyses have shown that the Gaia EDR3 parallax
zeropoint, $\varpi_{0}$, is observed to vary at the $\sim 30\,\mu\mathrm{as}$
level as a function of colour, magnitude, on-sky position and the type of
astrometric solution (Lindegren et al., 2021b). Ideally, all possible
variations would be included in the modelling here and the parallax zeropoint
behaviour simultaneously constrained. However, initial tests demonstrated that
magnitude dependence of $\varpi_{0}$ cannot be simultaneously constrained
alongside the period–luminosity relation. A similar phenomenon was reported by
Chan & Bovy (2020). In a similar vein, the variation of the parallax zeropoint
with on-sky position is degenerate with the on-sky distance prior variation
$p(s|\ell,b)$ (again see Chan & Bovy, 2020). Without additional information
(e.g. other tracer populations) the magnitude or on-sky dependence of the
zeropoint are not constrained so instead previously determined zeropoint
models are used with some additional colour dependence i.e.
$\varpi_{0,i}=\varpi^{f}_{0,i}(G,\ell,b)+\varpi^{e}_{0,i}(\nu_{\mathrm{eff}})$
where $i\in\\{5,6\\}$ denotes whether five- or six-parameter astrometric
solutions are considered. For the base model $\varpi^{f}_{0,i}$, three options
are used:

  1. 1.

the zeropoint corrections of Lindegren et al. (2021b) evaluated at
$\nu_{\mathrm{eff}}=1.25\,\mu\mathrm{m}^{-1}$ (this wavenumber is within the
interpolation grid for both five- and six-parameter solutions) accounting for
the $15\,\mu\mathrm{as}$ overestimate reported by Riess et al. (2021) and Zinn
(2021) for $G<10.8$,

  2. 2.

the colour-independent Healpix level 1 corrections from Groenewegen (2021)
(also incorporating the inflation of uncertainties he suggested) and

  3. 3.

the zeropoint model from Maíz Apellániz (2022) evaluated at
$\nu_{\mathrm{eff}}=1.25\,\mu\mathrm{m}^{-1}$.

For the additional modelled colour-dependent zeropoint,
$\varpi^{e}_{0,i}(\nu_{\mathrm{eff}})$, a quadratic is used with different
parameters for the five- and six-parameter solutions such that in summary the
model is

$\varpi_{0,i}(G,\nu_{\mathrm{eff}},\ell,b)=\varpi^{f}_{0,i}(G,\ell,b)+\sum_{j=0}^{j=2}q_{i,j}(\nu_{\mathrm{eff}}-1.1\,\mu\mathrm{m}^{-1})^{j}.$
(8)

There are then three free parameters $q_{i,j}$ for each of the five- and six-
parameter solutions.

####  4.4.1 Parallax uncertainty underestimate model

For the scaling factor of the parallax uncertainties, $f_{\varpi}$, two
quadratics in $G$ and $\nu_{\mathrm{eff}}$ for the five- and six-parameter
solutions are used:

$\ln
f_{\varpi,i}(G,\nu_{\mathrm{eff}})=\sum_{k,l\in\\{0,1,2\\}}r_{i,k,l}(G-14)^{k}(\nu_{\mathrm{eff}}-1.1\,\mu\mathrm{m}^{-1})^{l},$
(9)

where $i\in\\{5,6\\}$. This choice is motivated by the Gaia astrometric
performance being sensitive to colour and magnitude. As highlighted in Section
2, the parallax uncertainties may also be underestimated due to AGB
photocentre wobble. In Appendix B an additional parallax-dependent term is
included in $f_{\varpi}$ which does not affect the overall period–luminosity
relation fits.

###  4.5 Implementation

The models are implemented in Stan (Carpenter et al., 2017) using the python
interface PyStan 222Stan Development Team. 2018. PyStan: the Python interface
to Stan, Version 2.17.1.0. http://mc-stan.org.. The following priors are
adopted:

  1. 1.

$a_{m}\sim\mathcal{N}(a_{m,\mathrm{LMC}},0.5)$,

  2. 2.

$(b_{m},c_{m})\sim\mathcal{N}((b_{m,\mathrm{LMC}},c_{m,\mathrm{LMC}}),25\Sigma_{bc,\mathrm{LMC}})$,

  3. 3.

$F\sim\mathcal{N}(F_{\mathrm{LMC}},((3,1,1)\otimes(3,1,1))\Sigma_{F,\mathrm{LMC}})$
where $F=(\sigma_{2.3},m_{\sigma-},m_{\sigma+})$,

  4. 4.

$q_{i,j}\sim\mathcal{N}(0,1)$,

  5. 5.

$r_{i,k,l}\sim\mathcal{N}(0,3)$,

  6. 6.

$\sigma_{\varpi,0}\sim\mathcal{N}(-4.6,1.5)$,

  7. 7.

$\ln\sigma_{\mu,0,2}\sim\mathcal{N}(0.5,0.5)$,

  8. 8.

$\ln\vartheta_{2}\sim\mathcal{N}(-4.6,1.5)$,

  9. 9.

$\ln L_{0}\sim\mathcal{N}(1.1,0.6)$,

  10. 10.

$\boldsymbol{\tilde{S}}\sim\mathcal{N}(\boldsymbol{0},\tau\times\boldsymbol{1})$,
$\tau\sim\mathcal{C}(0,1)$ (a unit Cauchy prior),

  11. 11.

and when required $e\sim\mathcal{N}(e_{0},0.05e_{0})$ where $e_{0}$ is from
Wang & Chen (2019).

The LMC fits from Appendix C have been used as weak priors on the slopes
$(b_{m},c_{m})$ and error model parameters
$F=(\sigma_{2.3},m_{\sigma-},m_{\sigma+})$. For $(b_{m},c_{m})$ and
$\sigma_{2.3}$ a generous $5$ and $3$ times the LMC fit uncertainty is used
respectively as the prior width. Instead of performing the integration in
equation (2), the logarithm of the true parallax of each star minus the
zeropoint offset in magnitude, $-\ln
s_{i}-0.2\ln(10)(m_{\mathrm{abs}}(P_{i})-m_{\mathrm{abs,LMC}}(P_{i}))$, is
sampled (accounting for the additional Jacobian factor of $s$ due to sampling
in $\ln s$). This combination of parameters minimises the correlations in the
likelihood leading to more efficient sampling.

##  5 Results

Figure 6: Difference between parallaxes from the fitted O-rich Mira variable
$W_{Ks,J-Ks}$ period–luminosity relations (fourth row of Table 1) and the
zeropoint-corrected Gaia DR3 parallaxes (using the Lindegren et al., 2021b,
with an additional colour-dependent term). The median and uncertainty for 30
(15) equally-populated bins for 5(6)-parameter astrometric solutions are shown
as blue circles (orange squares). The annotation in each panel shows the
standard deviation of the estimates over the typical error (i.e. a measure of
any additional bias). Figure 7: Results of fitting the Gaia EDR3 parallax
zeropoint and the parallax uncertainty scaling factor. The results use the
Wesenheit $W_{Ks,J-Ks}$ magnitude relation (fourth row of Table 1). The top
(bottom) row corresponds to Gaia EDR3 5(6)-parameter solutions. The left plots
show the fitted colour-dependent parallax zeropoint variation in addition to
that reported by Lindegren et al. (2021b) at
$\nu_{\mathrm{eff}}=1.25\,\mu\mathrm{m}^{-1}$. The black points show the mean
difference between the corrected DR3 parallaxes and the parallax computed from
the period–luminosity relation for the data. The middle two panels show the
parallax uncertainty scaling factor as a function of $\nu_{\mathrm{eff}}$ and
$G$ (models from El-Badry et al. (2021) and Maíz Apellániz (2022) are shown).
The right panel shows the distribution of the parallax residual between the
period–luminosity relation and the corrected Gaia DR3 astrometry divided by
the combined error with (grey) and without (blue) the parallax scaling factor.
Figure 8: Comparison between the derived period–luminosity relation (using the
Wesenheit index $W_{Ks,J-Ks}$) and (i) VLBI parallax measurements for O-rich
Mira variables (orange crosses) and (ii) likely O-rich Mira variable globular
cluster members (dots coloured by their metallicity). The solid blue line and
shaded region gives the Milky Way model and its scatter, whilst the dashed
line shows the LMC measurement. The numbers show the median offsets (data –
model, as shown in the lower panel) with respect to the MW (and LMC) relations
for VLBI measurements. The uncertainty is the inverse-variance-weighted error
from the VLBI parallaxes, the photometric uncertainties and the scatter about
the period–luminosity relation (due to using single-epoch observations).

The results of the period–luminosity relation fitting are presented in Table 1
and the associated parameters for the Gaia EDR3 systematics in Table 2. The
default base parallax zeropoint model is option (i) from Section 4.4 that
primarily uses the correction from Lindegren et al. (2021b). As previously
reported elsewhere (see Iwanek et al., 2021a), the gradients, $b$ and $c$,
steepen for longer wavelengths. The scatter $\sigma_{2.3}$ also decreases with
wavelength. The Wesenheit models typically agree very well with those computed
using the single band models (e.g. $W_{Ks,J-Ks}=K_{s}-e(J-K_{s})$ for $a$
gives $a=-7.41$ compared to $a=-7.40$) suggesting circumstellar dust in the
O-rich Mira variables is unimportant (if it has a similar reddening law to the
interstellar medium). This tallies with the results of Bladh et al. (2015) who
showed using a grid of theoretical models that circumstellar dust around
O-rich stars is mostly transparent in optical and near-infrared bands. When
comparing the Milky Way results to linear fits of the LMC period–luminosity
relation (see Table 5), consistently fainter zeropoints (higher $a$) of the
Milky Way period–luminosity relation are found ($\Delta J,\Delta H,\Delta
K_{s})_{\log_{10}P=2.3}=(0.19,0.19,0.11)\,\mathrm{mag}$) but these differences
are well within the $0.5\,\mathrm{mag}$ prior width. Typically the gradients
($b$ and $c$) are found to be steeper for the Milky Way relations but it
doesn’t appear the broad LMC prior is causing any tension (possibly for $c$
for the $J$ and $H$ relations as illustrated in Fig. 18 although this may be
more linked to the selection of LMC sources). Note that due to the gradient
differences, the magnitude difference between the Milky Way and LMC
period–luminosity relations decreases with increasing period. As evidenced in
Fig. 18, literature quadratic model fits to the LMC Mira variables show
smaller offsets with respect to the Milky Way linear fits particularly around
the characteristic $200$ day period. It could be that more flexible models
produce less tension between the two period–luminosity relations. The
differences between the LMC and Milky Way relations in the context of their
population differences are discussed further in Section 5.2.

Fig. 6 shows the residuals of the parallaxes predicted from the $W_{Ks,J-Ks}$
relation from Table 1 compared to the zeropoint-corrected Gaia EDR3
parallaxes. We see in general the satisfactory agreement demonstrating the
quality of the period–luminosity relation. However, residuals and trends
remain. The left panel of Fig. 7 shows the fitted Gaia EDR3 zeropoint term for
this model. For 5-parameter solutions small corrections ($\lesssim
5\,\mu\mathrm{as}$) are required on top of the Lindegren et al. (2021b)
corrections. For 6-parameter solutions however, larger corrections are
required that typically increase as the sources get redder. This implies the
recommended zeropoint corrections evaluated at
$\nu_{\mathrm{eff}}=1.25\,\mu\mathrm{m}^{-1}$ do not apply well to redder
sources with 6-parameter solutions and appear to overcorrect the parallaxes.
Similar behaviour is found for the other models shown in Table 2. Fig. 7 shows
the factor by which the parallax uncertainties must be inflated to account for
the observed spread about the period–luminosity relation. In agreement with
previous work (e.g. El-Badry et al., 2021; Maíz Apellániz, 2022;
Andriantsaralaza et al., 2022) an inflation of the parallax uncertainties is
required. The behaviour is relatively flat with colour (although increases
quite steeply for very red sources with 6-parameter solutions). For
5-parameter solutions the factor is around $1.3$ for brighter $(G\sim 9)$ and
fainter ($G\sim 16$) sources but for more intermediate ($G\sim 12$ as reported
in Table 2) the factor increases to around $1.6$. This behaviour mirrors that
found by El-Badry et al. (2021) using wide binaries although larger factors
are found that are more consistent with the results of Maíz Apellániz (2022).
A fit using only five-parameter solutions from Gaia produces very similar
results for the Gaia systematic parameters and the period–luminosity relations
suggesting although the six-parameter solutions appear more biased, they are
not affecting the overall fit too strongly.

As shown in Fig. 6, some residuals in the fits remain, particularly as a
function of $G$ and on-sky location. In Section 5.2 possible population
effects producing such residuals are discussed. However, particularly in the
case of the residuals with $G$ where there are features around $G\approx 13$,
some level of residual at the $10\,\mu\mathrm{as}$ level appears to arise from
the Gaia EDR3 zeropoint model. The Groenewegen (2021) and Maíz Apellániz
(2022) zeropoint corrections have been used as variants of the base model. As
seen in Table 1, this can produce changes in the period–luminosity zeropoint
of $\sim 0.1\,\mathrm{mag}$. However, both of these alternatives also produce
larger residual features with $G$. The residual scatter is quantified using
the inverse-variance-weighted bin-to-bin scatter in the mean divided by the
mean uncertainty in the mean residual in each bin ($\sigma/\epsilon$). For the
five-parameter solutions binned as a function of $G$, the base $K_{s}$ model
produces $\sigma/\epsilon=1.6$ for the Lindegren et al. (2021b) model whilst
this inflates to $\sigma/\epsilon=2.3$ and $\sigma/\epsilon=2.4$ for
Groenewegen (2021) and Maíz Apellániz (2022) models respectively. The largest
problems occur around $G\approx 12-13$. As noted previously, simultaneously
fitting the magnitude (and on-sky dependence) of the parallax zeropoint was
found to be degenerate with parameters of the period–luminosity relation. A
future approach should adopt a more flexible model for the parallax zeropoint
constrained to be small by a careful choice of prior.

Table 1 also displays results for the Yuan et al. (2013) extinction law. As
with the case using the Wang & Chen (2019) extinction law, the Wesenheit
magnitude zeropoint is very similar ($\lesssim 0.01\,\mathrm{mag}$) to that
computed using the single band results suggesting the adopted extinction law
doesn’t change the conclusions significantly. The sensitivity to the RUWE cut
(by default $1.4$) has been investigated. Relaxing to RUWE $<2$ produces a
slightly steeper fainter $K_{s}$ relation that is consistent with the RUWE
$<1.4$ relation for $P>200$ day but deviates slightly at the shorter period
end. Many of the higher RUWE stars are located near the midplane and so
potentially are affected by high source density. Results are also reported for
C-rich Mira variables. As done in Appendix C for the C-rich LMC Mira
variables, a quadratic period–luminosity relation
$m_{\mathrm{abs}}(P)=a+b(\log_{10}P-2.3)+c(\log_{10}P-2.3)^{2}$ with a linear
scatter $\sigma_{\mu}(P)=\sigma_{2.3}+m_{\sigma-}(\log_{10}P-2.3)$ is used.
C-rich Mira variables are typically not employed as distance indicators due to
their larger scatter in the period–luminosity relation compared to the O-rich
Mira variables. Here it is found that in the Wesenheit magnitude $W_{Ks,J-Ks}$
the C-rich Mira variables at short periods ($\lesssim 300\,\mathrm{days}$) are
$\sim 0.4\,\mathrm{mag}$ brighter than the O-rich relations (also seen in the
LMC, Appendix C) and the scatter is comparable to that of the O-rich Mira
variables. At longer periods ($\gtrsim 400\,\mathrm{days}$) the
period–luminosity relation flattens (or possibly even turns over, see Appendix
C).

Table 1: Period–luminosity relations for O-rich Mira variables. The period–luminosity relations have the form $a+b(\log_{10}P-2.3)$ for $\log_{10}P\leq 2.6$ and $a+0.3b+c(\log_{10}P-2.6)$ for $\log_{10}P>2.6$ with scatter $\sigma=\sigma_{2.3}+m_{\sigma-}(\log_{10}P-2.3)$ for $\log_{10}P\leq 2.6$ and $\sigma=\sigma_{2.3}+0.3m_{\sigma-}+m_{\sigma+}(\log_{10}P-2.6)$ for $\log_{10}P>2.6$ (note for the C-rich relation a quadratic relation $a+b(\log_{10}P-2.3)+c(\log_{10}P-2.3)^{2}$ with a linear scatter $\sigma=\sigma_{2.3}+m_{\sigma-}(\log_{10}P-2.3)$ for all periods is used instead). Here $P$ is in days. $L_{0}$ is the logarithm of the mean of the exponential of the distance prior scalelength in kpc. The first section of rows show results for the 2MASS $JHK_{s}$ bands and using the Wesenheit indices $W_{x,y-x}=x-e(y-x)$. All of these models use the default setup correcting the Gaia EDR3 parallaxes using the Lindegren et al. (2021b) zeropoints evaluated at $\nu_{\mathrm{eff}}=1.25\,\mu\mathrm{m}^{-1}$ as a base model and fitting for an additional colour-dependent term. The second section shows model variations: (i) using the Yuan et al. (2013) extinction coefficients for the $W_{Ks,J-Ks}$ relation, (ii) allowing the extinction coefficient $e$ to vary for the $W_{Ks,J-Ks}$ relation, (iii) using stars with RUWE $<2$ for the $K_{s}$ relation, (iv) using the Groenewegen (2021) parallax zeropoint correction as a base model for the $K_{s}$ relation and (v) using the Maíz Apellániz (2022) parallax zeropoint correction evaluated at $\nu_{\mathrm{eff}}=1.25\,\mu\mathrm{m}^{-1}$ as a base model for the $K_{s}$ relation. The final section gives the $W_{Ks,J-Ks}$ relation for C-rich stars using the Lindegren et al. (2021b) zeropoints as a base model. Band/Model | $a$ | $b$ | $c$ | $\ln\sigma_{2.3}$ | $m_{\sigma-}$ | $m_{\sigma+}$ | $L_{0}$ | $e$
---|---|---|---|---|---|---|---|---
$J$ | $-5.66\pm 0.02$ | $-3.56\pm 0.06$ | $-2.42\pm 0.48$ | $-1.75\pm 0.06$ | $-0.02\pm 0.04$ | $2.35\pm 0.30$ | $1.51\pm 0.05$ | $-$
$H$ | $-6.46\pm 0.02$ | $-3.84\pm 0.06$ | $-4.08\pm 0.48$ | $-1.77\pm 0.07$ | $0.06\pm 0.04$ | $1.86\pm 0.25$ | $1.51\pm 0.05$ | $-$
$K_{s}$ | $-6.85\pm 0.02$ | $-4.22\pm 0.06$ | $-5.52\pm 0.47$ | $-1.84\pm 0.07$ | $0.04\pm 0.04$ | $1.58\pm 0.28$ | $1.51\pm 0.05$ | $-$
$W_{Ks,J-Ks}$ | $-7.40\pm 0.02$ | $-4.52\pm 0.06$ | $-7.06\pm 0.45$ | $-1.90\pm 0.06$ | $0.01\pm 0.04$ | $1.42\pm 0.24$ | $1.50\pm 0.05$ | $0.47$
$W_{Ks,H-Ks}$ | $-7.40\pm 0.02$ | $-4.76\pm 0.06$ | $-7.81\pm 0.48$ | $-1.91\pm 0.07$ | $-0.04\pm 0.04$ | $0.97\pm 0.22$ | $1.50\pm 0.05$ | $1.47$
$W_{H,J-H}$ | $-7.34\pm 0.02$ | $-4.11\pm 0.07$ | $-6.13\pm 0.49$ | $-1.85\pm 0.05$ | $0.05\pm 0.03$ | $1.55\pm 0.24$ | $1.50\pm 0.05$ | $1.17$
Yuan $e_{i}$ $W_{Ks,J-Ks}$ | $-7.74\pm 0.02$ | $-4.66\pm 0.06$ | $-8.04\pm 0.46$ | $-1.85\pm 0.07$ | $0.02\pm 0.03$ | $1.38\pm 0.24$ | $1.51\pm 0.05$ | $0.74$
Free $e_{i}$ $W_{Ks,J-Ks}$ | $-7.36\pm 0.02$ | $-4.51\pm 0.06$ | $-6.95\pm 0.46$ | $-1.87\pm 0.05$ | $0.02\pm 0.04$ | $1.39\pm 0.25$ | $1.50\pm 0.05$ | $0.45\pm 0.02$
RUWE $<2$ $K_{s}$ | $-6.80\pm 0.02$ | $-4.37\pm 0.06$ | $-5.11\pm 0.38$ | $-1.88\pm 0.06$ | $-0.01\pm 0.03$ | $1.72\pm 0.27$ | $1.61\pm 0.04$ | $-$
G21 $K_{s}$ | $-6.73\pm 0.02$ | $-4.12\pm 0.06$ | $-5.53\pm 0.47$ | $-1.79\pm 0.06$ | $0.06\pm 0.04$ | $1.56\pm 0.29$ | $1.42\pm 0.05$ | $-$
MA22 $K_{s}$ | $-6.76\pm 0.02$ | $-4.18\pm 0.06$ | $-5.54\pm 0.46$ | $-1.83\pm 0.06$ | $0.04\pm 0.03$ | $1.61\pm 0.28$ | $1.43\pm 0.05$ | $-$
C-rich $W_{Ks,J-Ks}$ | $-7.73\pm 0.09$ | $-4.00\pm 0.51$ | $0.59\pm 1.07$ | $-1.72\pm 0.15$ | $0.48\pm 0.04$ | $-$ | $2.64\pm 0.12$ | $0.47$
Table 2:  Parallax zeropoint and uncertainty model results for the models shown in Table 1. $\varpi_{i,\mathrm{zp}}$ gives the sky-averaged Gaia EDR3 parallax zeropoint in addition to the assumed model for $i$-parameter astrometric solutions (at $\nu_{\mathrm{eff}}=1.15\,\mu m^{-1}$ for 5-parameter solutions and $\nu_{\mathrm{eff}}=1.05\,\mu m^{-1}$ for 6) in units of $\mu\mathrm{as}$. $f_{i,\varpi}$ gives the scaling of the parallax errors for the $i$-parameter solutions (again at the representative colours and magnitude $G=12$). The assumed base parallax zeropoint model is by default the Lindegren et al. (2021b) correction at $\nu_{\mathrm{eff}}=1.25\,\mu m^{-1}$ except for G21 that uses the Groenewegen (2021) parallax zeropoint correction and MA22 that uses Maíz Apellániz (2022) parallax zeropoint correction evaluated at $\nu_{\mathrm{eff}}=1.25\,\mu\mathrm{m}^{-1}$. Band/Model | $\varpi_{5,\mathrm{zp}}$ | $\varpi_{6,\mathrm{zp}}$ | $f_{5,\varpi}$ | $f_{6,\varpi}$
---|---|---|---|---
$J$ | $-2\pm 1$ | $21\pm 2$ | $1.52\pm 0.04$ | $1.62\pm 0.02$
$H$ | $-2\pm 1$ | $22\pm 3$ | $1.56\pm 0.04$ | $1.61\pm 0.02$
$K_{s}$ | $-3\pm 1$ | $17\pm 3$ | $1.58\pm 0.04$ | $1.61\pm 0.02$
$W_{Ks,J-Ks}$ | $-4\pm 2$ | $15\pm 2$ | $1.61\pm 0.04$ | $1.60\pm 0.02$
$W_{Ks,H-Ks}$ | $-6\pm 2$ | $10\pm 3$ | $1.60\pm 0.03$ | $1.59\pm 0.02$
$W_{H,J-H}$ | $-5\pm 1$ | $18\pm 2$ | $1.58\pm 0.05$ | $1.62\pm 0.02$
Yuan $e_{i}$ $W_{Ks,J-Ks}$ | $-6\pm 1$ | $10\pm 3$ | $1.62\pm 0.04$ | $1.59\pm 0.02$
Free $e_{i}$ $W_{Ks,J-Ks}$ | $-4\pm 1$ | $15\pm 3$ | $1.60\pm 0.04$ | $1.60\pm 0.03$
RUWE $<2$ $K_{s}$ | $-4\pm 1$ | $15\pm 2$ | $1.61\pm 0.03$ | $1.59\pm 0.02$
G21 $K_{s}$ | $-25\pm 1$ | $0\pm 3$ | $1.57\pm 0.05$ | $1.59\pm 0.02$
MA22 $K_{s}$ | $-1\pm 1$ | $9\pm 3$ | $1.60\pm 0.04$ | $1.61\pm 0.02$
C-rich $W_{Ks,J-Ks}$ | $26\pm 9$ | $16\pm 15$ | $1.49\pm 0.30$ | $0.35\pm 0.18$

###  5.1 Comparison with VLBI parallaxes

An alternative to the astrometric distances of Mira variables from Gaia are
interferometric measurements from very long-baseline interferometry (VLBI). As
VLBI is able to resolve AGB stars, any systematics from photocentre wobble are
minimal (see Section 3). In combination with Hipparcos parallaxes, Whitelock
et al. (2008) used the available VLBI measurements to calibrate the $K$-band
period–luminosity relation. Since then, several more AGB stars have had VLBI
measurements. Andriantsaralaza et al. (2022) has inspected the Gaia DR3
astrometry of AGB stars with VLBI measurements. Fig. 8 displays the absolute
$W_{Ks,J-Ks}$ measurements against period for the recent VLBI compilations of
AGB stars from Xu et al. (2019) and VERA Collaboration et al. (2020),
preferentially using the results from VERA Collaboration et al. (2020) in the
case of duplicates. The periods are from VSX (Watson et al., 2006) and
magnitudes from 2MASS. Only O-rich Mira variables as defined by the selection
in Section 2 are displayed. FV Boo is removed as it appears to be a clear
outlier as noted by Kamezaki et al. (2016) and there are concerns it displays
additional variability due to potentially being in a binary system (Kamezaki
et al., 2016). The inverse-variance-weighted offset of the absolute Wesenheit
magnitudes computed using VLBI parallaxes with respect to the
period–luminosity relation is $(0.12\pm 0.11)\,\mathrm{mag}$. Here the error
is the inverse-variance-weighted error from the photometric uncertainties, the
VLBI parallax uncertainties and the scatter model due to using single epoch
observations. Although the measurements are consistent, the VLBI measurements
are slightly fainter than the Gaia-derived Milky Way trend, possibly as they
are a dustier or a more metal-rich population compared to the Gaia-selected
O-rich Mira variables (also seen in Whitelock et al., 2008). A concern is that
many of the 2MASS measurements are saturated for these bright stars. Whitelock
et al. (2000) and Whitelock et al. (2008) provide $JHK$ measurements in the
SAAO system. Transformation to the 2MASS system is not simple for these very
red sources but using the relations in Koen et al. (2007) the offset with
respect to the derived period–luminosity relation is $(0.15\pm
0.05)\,\mathrm{mag}$. However, it should be noted that Koen et al. (2007) find
brighter stars appear to have larger differences between SAAO $K$ and 2MASS
$K_{s}$ ($K_{s}$ smaller than $K$) which could explain some of this
difference.

###  5.2 Population variations

It has been found that the Milky Way O-rich Mira variable relations derived
here are typically slightly fainter than those derived for the LMC (see
Appendix C) particularly at the short period end due to a steeper gradient.
One interpretation of this result is that there is variation of the O-rich
Mira period–luminosity relation with stellar population, in particular with
the age and metallicity of the population. Typically, it has been found that
population effects are quite minimal for the Mira variables, particularly in
the near- and mid-infrared ($K_{s}$, $[3.6]$ and $[4.5]$, Whitelock et al.,
2008; Goldman et al., 2019; Menzies et al., 2019) or using bolometric
magnitudes (e.g. Andriantsaralaza et al., 2022). However, there are
suggestions from theoretical results that there can be more significant
variations in the period–luminosity relations (Wood, 1990; Qin et al., 2018)
particularly for the bluer bands, $J$ and $H$, that are also investigated
here.

####  5.2.1 Comparison with theoretical models

Fundamentally, it is expected that a given mass and radius combination will
give rise to the same fundamental period. Wood (1990) demonstrated using a
linear calculation how the period of a Mira variable is related to the
luminosity $L$, metallicity $Z$ and mass as $M$ as $P\propto
L^{1.59}Z^{0.46}M^{-1.55}$. If it is assumed that an AGB star will only
pulsate with Mira-like oscillations when it reaches a certain radius (or
narrow radial range) for its given mass, this gives us a relationship between
bolometric magnitude $M_{\mathrm{bol}}$ and metallicity at fixed radius
$\Delta M_{\mathrm{bol}}=0.72\Delta\log_{10}Z$ (see also figure 12 of
Trabucchi et al., 2019, for a similar calculation with a very similar result).
As noted by Wood (1990), the corresponding change in near-infrared magnitudes
with metallicity is smaller than the change in bolometric magnitude. Assuming
Mira variables of fixed radius but different metallicities are black-bodies
with varying effective temperatures $\Delta\log_{10}T_{\mathrm{eff}}\approx
0.072\Delta\log_{10}Z$ the magnitude differences are $(\Delta M_{J},\Delta
M_{H},\Delta M_{Ks})=(0.68,0.52,0.42)\Delta\log_{10}Z$. Taking the typical
$Z_{\mathrm{LMC}}=0.5Z_{\mathrm{MW}}$, the magnitude differences are $(\Delta
M_{J},\Delta M_{H},\Delta M_{Ks})=(0.20,0.16,0.13)$ in rough agreement with
the zeropoint differences found.

It is anticipated that linear calculations will differ most strongly from non-
linear calculations in the computation of period at a given mass and radius
(Trabucchi et al., 2021a) making these arguments valid irrespective or whether
linear or non-linear calculations are considered. However, Trabucchi et al.
(2019) has shown that, particularly for the fundamental mode, the composition
(metallicity, C/O ratio) can affect the period at fixed mass and radius. For
instance, making a star more metal-rich (increasing from typical LMC to
typical Milky Way metallicity) or making a star carbon-rich (increasing C/O
from $0.55$ to $\sim 3$) decreases the period by $\sim 10\,\mathrm{per\,cent}$
(for a linear calculation). Therefore, period is not solely a function of mass
and radius. In a similar vein, Feast (1996) has questioned the validity of the
assumption that a star of given mass reaches Mira-like oscillations at fixed
radius independent of its metallicity as it is related to the mass loss. For a
given initial mass and metallicity, an AGB star could reach the Mira pulsation
stage with a different mass-radius combination that produces a similar period.
However, there is evidence to suggest metallicity-dependence on mass loss is
not a significant effect (see Höfner & Olofsson, 2018, for a summary).

Using $P\propto L^{1.59}Z^{0.46}M^{-1.55}$ and the period-mass-radius
relation, the dependence of the effective temperature can be derived as
$T_{\mathrm{eff}}\propto P^{-0.1}Z^{-0.073}M^{0.014}$ demonstrating that at
fixed period the effective temperature is a weak function of the mass and more
dependent upon metallicity. This then suggests even when the mass evolution at
a given metallicity is poorly known, the metallicity of a Mira variable of
fixed period will be related to its effective temperature and hence infrared
colours (this is corroborated by the fuller calculation of Qin et al., 2018,
that is considered later and that shows $J$, $H$ and $K_{s}$ at fixed period
all have similar age dependence such that the gradient of $J-K_{s}$ with age
is $\lesssim 0.002\,\mathrm{mag}/\mathrm{Gyr}$). Using the blackbody model
from before, the colour difference is found to be $\Delta(J-K_{s})\approx
0.26\Delta\log_{10}Z=-3.56\log_{10}T_{\mathrm{eff}}$. This is in agreement
with PARSEC isochrones (Bressan et al., 2012; Marigo et al., 2017) which
suggest $\mathrm{d}(J-K_{s})/\mathrm{d}[\mathrm{Fe}/\mathrm{H}]\approx 0.2$.
For the LMC sample, the mean colour $(J-K_{s})\approx 1.11$ at
$\log_{10}P=2.3$ whilst for the Milky Way sample it is $\sim 1.2$, which using
the simplistic approach would translate into a $\sim 0.4\,\mathrm{dex}$
metallicity shift.

It seems from simple considerations that the derived differences between the
LMC and Milky Way relations are consistent with linear pulsation calculations.
However, the Wood (1990) formulae have been criticized by Feast (1992) as they
fail to simultaneously explain the period-colour relation in the Milky Way/LMC
and the period–metallicity relation observed in globular cluster Mira
variables (Feast & Whitelock, 2000a). Fig. 8 displays possible globular
cluster members taken from the main Milky Way sample defined as within $3$
half-light radii of a known globular cluster (Harris, 2010) with proper
motions in each component consistent at the $4\sigma$ level with those
determined by Baumgardt & Vasiliev (2021). It is clear this generous cross-
match introduces a couple of non-members. A globular cluster
period–metallicity gradient is visible where there is a collection of metal-
poor stars at around $140$ day periods and a collection of more metal-rich
stars at $300$ day periods. This is slightly puzzling but it should be noted
that some globular clusters show Mira variables with a range of periods
(Matsunaga & IRSF/SIRIUS Team, 2007) suggesting we are seeing the effects of
age-metallicity correlations and/or the impact of multiple populations in
globular clusters.

Figure 9: Differences between MW and LMC O-rich Mira variable
period–luminosity relations at $\log_{10}P/\,\mathrm{day}=2.4$ in the
$JHK_{s}$ bands compared to the linear pulsation models reported by Qin et al.
(2018). A set of models with different age differences (as numbered in Gyr)
and $[\mathrm{Fe}/\mathrm{H}]$ (coloured) are shown. The upwards diagonal
sequences depict $\Delta M_{H}$ as a function of $\Delta M_{J}$ whilst the
downwards diagonal sequences depict $\Delta M_{Ks}$. Given the measured
differences the MW O-rich Mira variable population is found to be $(0.4\pm
0.1)\,\mathrm{dex}$ more metal-rich and $(4.9\pm 0.3)\,\mathrm{Gyr}$ older
than the LMC population considered.

The previous arguments explained in simple terms why both magnitude and colour
differences with varying metallicity at fixed period are to be expected for
Mira variables. This can be elucidated further with a more sophisticated
model. Qin et al. (2018) have used the linear pulsation models from Wood &
Olivier (2014) combined with a relation for mass as a function of age,
metallicity and helium abundance from Nataf et al. (2012) and the bolometric
corrections from Casagrande & VandenBerg (2014) to derive gradients of
$JHK_{s}$ magnitude with these quantities at fixed period ($\log_{10}P=2.4$,
although they report similar gradients for other periods in the near infrared
bands). These authors caution that the models are approximate and do not seem
to explain the differences between Mira variables in the Galactic bulge and
the LMC. Indeed, at fixed age and helium abundance, the models predict
brighter $K_{s}$ with metallicity in contrast to the previous discussion.
Nonetheless, in the absence of other models, they are used here. Again,
although the period for a given mass and radius combination is affected by the
linear approximation (e.g. Trabucchi et al., 2021a), the gradient of magnitude
with age and metallicity at fixed period is more related to the gross stellar
evolutionary properties. The models from Qin et al. (2018) are used to infer
the age and metallicity difference between the Milky Way population and LMC
population (see Appendix C) as shown in Fig. 9. Here it is assumed the helium
abundance is similar in both systems. The combination of $J$ and $H$
differences provides little leverage for breaking age/metallicity differences
but when combined with the comparatively smaller $K_{s}$ difference the LMC
O-rich Mira variable population is found to be younger by $(4.9\pm
0.3)\,\mathrm{Gyr}$ and more metal-poor by $(0.4\pm 0.1)\,\mathrm{dex}$,
somewhat consistent with expectation. There is evidence for a gap in the star
formation history of the LMC and an increase in the star formation rate in the
last $\sim 1\,\mathrm{Gyr}$ based on the properties of its star clusters
(Jensen et al., 1988), its chemical evolution (e.g. Hasselquist et al., 2021)
and its photometrically-derived star formation history (Javiel et al., 2005).

Further evidence for variation in the zeropoint with metallicity (or more
generally stellar population) comes from the globular clusters. Fig. 8
demonstrates that there is a weak tendency for the globular cluster members to
get brighter as a function of metallicity relative to the LMC and Milky Way
relations (or putting it another way, the globular clusters alone suggest a
flatter period–luminosity slope). The lack of metal-rich shorter-period and
metal-poor longer-period globular cluster Mira variables makes this conclusion
somewhat uncertain. Using a globular cluster-calibrated period–luminosity
relation, Feast et al. (2002) find a distance modulus for the LMC $\sim
0.1\,\mathrm{mag}$ further than modern estimates suggest and Whitelock et al.
(2008) find the $K_{s}$ period–luminosity relation for globular cluster
members is $\sim 0.1$ brighter than the LMC relation (using the Pietrzyński et
al., 2019, LMC distance modulus), but in both cases the uncertainties were
large. Finally, in Appendix C the period–luminosity relations for the
Sagittarius dwarf spheroidal galaxy (Sgr dSph) and the Small Magellanic Cloud
(SMC) are estimated. It is found that typically the (relatively few) O-rich
Mira variables in these systems are slightly brighter than their presumably
more metal-rich counterparts in the LMC in all bands particularly for periods
greater than $250$ days (corroborating the results of Ita et al., 2004). The
steep period–luminosity relations typically found for the SMC mean for stars
with periods less than $200$ days the SMC Mira variables are fainter than
those in the LMC but these stars are comparatively rare.

####  5.2.2 Population gradients within the samples

We have seen how differences in period–luminosity relations between systems
can be explained by population differences. However, the populations in the
LMC and Milky Way are not homogeneous so similar gradients should be observed
within these systems.

Fig. 7 shows the variation of the zeropoint-corrected Gaia EDR3 parallax
residual with respect to the estimates from the $W_{Ks,J-Ks}$ model of Table
1. We see there is a tendency for the outer parts of the Galaxy to have larger
Gaia parallaxes (smaller distances) than the period–luminosity relations
suggest. This implies that for the outer disc the absolute $W_{Ks,J-Ks}$ needs
to be fainter. Using the Qin et al. (2018) relations, inside-out formation (a
negative age gradient with radius, Frankel et al., 2019; Grady et al., 2019)
would imply $W_{Ks,J-Ks}$ gets brighter with Galactocentric radius, but a
negative radial metallicity gradient produces the opposite effect although
with a too weak $0.03\,\mathrm{mag}/\mathrm{dex}$ gradient. Neither age nor
metallicity effects appear to explain the observations, although the exact
slope reported by Qin et al. (2018) depends on the uncertain bolometric
corrections for cool stars (Casagrande & VandenBerg, 2014) and Qin et al.
(2018) themselves find inconsistencies between the theoretical models and the
expectations for Mira variables in the Galactic bulge. The effect we are
seeing could be driven by C-rich contamination that is more prevalent in the
outer-disc. There are some very red stars ($H-K_{s}>0.7$) even after
extinction correction. Typically removal of these redder sources makes the
long period end of the period–luminosity relation brighter (note the bias in
Fig. 6 at long periods which is somewhat alleviated by removing very dusty
sources) but the trends with Galactocentric radius remain. A further cause
could be incorrect extinction correction but there is no trend in the parallax
residuals against extinction. It is clear from Fig. 7 that systematic trends
in $G$ and on-sky position are present (the inner and outer Galaxy samples
have different mean $G$ magnitudes) and so potentially the cause of the
Galactocentric radius trend is remaining systematics in the Gaia parallaxes
and not due to any population differences.

As previously highlighted, the metallicity of giant stars correlates well with
their colour (Qin et al. 2018 suggest that colours at fixed period are
insensitive to age variations, $\lesssim 0.002\,\mathrm{mag}/\mathrm{Gyr}$,
and nearly completely depend upon helium abundance and metallicity). Here, the
impact of a colour term in the period–luminosity relations is investigated.
Table 1 gives the result of fitting the Wesenheit magnitude
$W_{Ks,J-Ks}=K_{s}-e(J-K_{s})$ with $e$ a free parameter finding $e=(0.45\pm
0.02)$ fully consistent with the estimate from interstellar extinction
considerations (0.47, Wang & Chen, 2019). This gives no evidence that there is
additional colour dependence and in turn metallicity dependence to the O-rich
period–luminosity relation. However, this simple approach uses the extincted
$J$ and $K_{s}$ magnitudes in the modelling. Instead including an additional
extinction-corrected colour term $b_{JK}(J-K_{s})$ in the $K_{s}$
period–luminosity relation, the best-fitting gradient is found as
$b_{JK}=(0.34\pm 0.05)$ giving evidence that the period–luminosity relation is
fainter for redder (more metal-rich) stars. However, the remaining colour-
magnitude-spatial correlations in the Gaia zeropoints make this conclusion
uncertain.

As discussed in Appendix C, there is also evidence in the LMC sample for a
metallicity gradient to the period–luminosity relation with more metal-rich
stars being fainter although this interpretation is somewhat complicated by
age-metallicity correlations. However, again assuming colours are age-
insensitive, the $(J-K_{s})$ colour gradient to the $K_{s}$ period–luminosity
relation is $b_{JK}=(0.45\pm 0.07)$ or using
$\mathrm{d}(J-K_{s})/\mathrm{d}[\mathrm{Fe}/\mathrm{H}]\approx 0.2$ the
gradient with metallicity is $(0.09\pm 0.02)\,\mathrm{dex}^{-1}$. This is in
rough agreement with the differences found between the MW and LMC systems as a
whole and consistent with the population gradient in the Milky Way sample.

A further check of metallicity dependence of the period–luminosity relation is
through analysis of the Galactic bulge Mira variables (Groenewegen &
Blommaert, 2005; Qin et al., 2018). The period–luminosity relation can be
calibrated under the assumption that the spatial distribution peaks around the
now well-determined distance of Sgr A* (Gravity Collaboration et al., 2021).
However, these bulge stars are more sensitive to extinction assumptions and
modelling the distance distribution requires good knowledge of the selection
function. Finally, in the Galactic disc, the period–luminosity relation could
be inspected as a function of kinematics which acts as a proxy for
age/metallicity. Alvarez et al. (1997) reported differences in the
period–luminosity relation for different kinematically-defined populations
using Hipparcos data. Both of these avenues require further investigation that
is deferred to future work. In conclusion, there is evidence from both the
mean difference between the LMC and Milky Way and from differences within the
LMC and Milky Way samples of a metallicity gradient to the period–luminosity
relations for O-rich Mira variables with the more metal-rich stars
intrinsically fainter than the metal-poor as expected from theoretical
studies.

##  6 Consequences for the Hubble constant

Our period–luminosity relations for Milky Way O-rich Mira variables provide
useful anchors for the Type Ia supernova Hubble diagram and in turn a
measurement of the Hubble constant. Currently the only SNIa host with observed
Mira variables is NGC 1559 (Huang et al., 2020) so Mira-based Hubble constant
measurements are limited primarily by the uncertainty on the properties of the
single supernova. However, over the coming years more observations of Mira
variables in other SN Ia host galaxies are expected, so reducing the sources
of uncertainty in the period–luminosity calibrations will become increasingly
important. Here measurements of the Hubble constant are provided largely
following the analysis of Huang et al. (2020) but replacing their
period–luminosity relations with those derived here. In addition to the Milky
Way relations, the LMC period–luminosity relations and Mira variables in the
water maser host galaxy NGC 4258 are used as further anchors.

NGC 1559 hosted the Type Ia supernova SN 2005df with peak magnitude
$m_{B}=(12.14\pm 0.11)\,\mathrm{mag}$ (Scolnic et al., 2018). Given a distance
modulus to NGC 1559, $\mu_{1559}$, the Hubble constant is estimated as

$\log_{10}H_{0}=\tfrac{1}{5}(m_{B}+5a_{B}+25)-\tfrac{1}{5}\mu_{1559},$ (10)

where $a_{B}=(0.71273\pm 0.00176)$ is the SNIa magnitude-redshift intercept as
measured by Riess et al. (2016). It is beyond the scope of this work to
combine the Type Ia supernovae modelling with the anchors in a probabilistic
model as done by Riess et al. (2022a) but the adopted $a_{B}$ encompasses the
range of fits from Riess et al. (2022a) and alters $H_{0}$ by $\sim
0.2\,\mathrm{km\,s}^{-1}\mathrm{kpc}^{-1}$. The model for the Mira variables
in NGC 1559 as presented by Huang et al. (2020) is first described and then
used to derive the estimate of $H_{0}$.

###  6.1 Basic model and data

The NGC 1559 Mira variables are taken from Huang et al. (2020) and the NGC
4258 Mira variables are from Huang et al. (2018). For both samples, mean
magnitudes (and for NGC 1559 uncertainties) in the Hubble WFC3 $F160W$ band
are provided along with period estimates. Both samples are defined to have
peak-to-trough $F160W$ amplitude between $0.4$ and $0.8$ (to reduce C-rich
contamination as discussed later). NGC 4258 has an additional colour cut
($F125W-F160W<1.3$ equivalent to $J-H<2.2$ using the colour transformations
from the X-Shooter spectra as described below) which is relatively mild as for
the LMC Mira variable sample it only removes $2$ of $907$ Mira variables with
$P<300$ days (independent of whether extinction corrections are applied). For
the NGC 4258 sample, there are further cuts on $F814W$ detection and
variability to define a ‘Silver’ and ‘Gold’ sample respectively. For the NGC
1559 sample these colour and variability cuts are not possible due to the lack
of multiband data. However, as a quality cut, sources in NGC 1559 with
crowding corrections $>0.25\,\mathrm{mag}$ are removed.

The $F160W$ magnitudes are corrected for foreground extinction of
$E(B-V)=0.0298$ for NGC 1559 and $E(B-V)=0.0163$ for NGC 4258 in Schlegel et
al. (1998) units using the extinction coefficients from Wang & Chen (2019) and
the $F160W$ uncertainties are broadened by a $16\,\mathrm{per\,cent}$
uncertainty in $E(B-V)$ and a $2.5\,\mathrm{per\,cent}$ uncertainty in the
$F160W$ coefficient (the systematic uncertainty on the derived NGC 1559 and
NGC 4258 distance moduli arising from the uncertainty in the extinction is
$\sim 0.002\,\mathrm{mag}$ so negligible compared to other sources of
uncertainty). This ignores any extinction within the systems. The
uncertainties on the periods of the Mira variables are ignored as they are not
provided and for near-linear models period uncertainties are equivalent to an
additional intrinsic magnitude spread (for approximately constant period
uncertainties).

For each galaxy’s Mira variable sample, a two-component Gaussian mixture model
is fitted to the residuals of the $F160W$ magnitudes with respect to the
period–luminosity relation (shifted by the distance modulus $\mu$) as

$\begin{split}p(F160W|P)=\sum_{j=1}^{j=2}\vartheta_{j}\mathcal{N}(F160W|&a+b(\log_{10}P-2.3)+\mu_{j},\\\
&f^{2}\sigma_{F160W}^{2}+\sigma_{0,j}^{2}).\end{split}$ (11)

All considered Mira variables have $P<400\,\mathrm{day}$ so only a linear
model is considered. The mixture model allows for a contribution from outliers
that do not follow a tight period–luminosity relation. An initial
consideration is that the Milky Way (and LMC) period–luminosity relations are
derived in the 2MASS $JHK_{s}$ bands whilst the extragalactic Mira variable
observations have been made in the Hubble WFC3 $F160W$ band (effective
wavelength of $1.528\,\mu\mathrm{m}$ compared to $J$ of $1.235\,\mu\mathrm{m}$
and $H$ of $1.662\,\mu\mathrm{m}$, Huang et al., 2018, 2020). Following Huang
et al. (2020) a colour term is used to convert 2MASS $H$ magnitudes into
$F160W$ magnitudes. $43$ stars in the O-rich Mira sample with periods $<400$
days are taken from the second release of the X-Shooter Spectral Library
(Gonneau et al., 2020). Using the filters provided by the SVO filter service
(Rodrigo et al., 2012; Rodrigo & Solano, 2020), the expected magnitudes of
these stars in the 2MASS filters and $F160W$ are found. The expected $J$ and
$H$ 2MASS magnitudes are on average $0.07$ and $0.13$ mag brighter than
measured in agreement with the comparison from Gonneau et al. (2020). The
broad-band colours are extinction corrected using the procedure described in
Section 4.2 using the interpolated $F160W$ coefficient
$A_{F160W}/A_{V}=0.1556$ from Wang & Chen (2019). The relationship between the
$F160W$ and 2MASS bands is found to be $F160W=H+c_{JH}(J-H)$ with
$c_{JH}=(0.379\pm 0.012)$ which agrees well with the colour coefficient from
Huang et al. (2020). Using Table 1, this implies a period–luminosity relation
for the $F160W$ band of

$F160W_{\mathrm{abs}}(P)=(-6.16\pm 0.02)+(-3.73\pm 0.05)(\log_{10}P-2.3),$
(12)

for O-rich Mira variables with $P<400\,\mathrm{days}$. The unknown $F160W$
period–luminosity relation is modelled probabilistically by allowing the
parameters $h=(a,b)$ and $c_{JH}$ to vary and including a ‘prior’ term of the
form

$p(\bar{h}|h,c_{JH})=-\tfrac{1}{2}(h-\bar{h})^{\mathrm{T}}\Sigma_{h}^{-1}(h-\bar{h}).$
(13)

Here $\bar{h}=\langle h_{H}+c_{JH}(h_{J}-h_{H})\rangle$,
$\Sigma_{f,a,a}=\mathrm{Var}(a_{H})$ as the uncertainty in the zeropoint is
assumed to be wholly driven by distance uncertainties,
$\Sigma_{h,a,b}=\mathrm{Cov}(a_{H}+c_{JH}(a_{J}-a_{H}),b_{H}+c_{JH}(b_{J}-b_{H}))$
and $\Sigma_{h,b,b}=\mathrm{Var}(b_{H}+c_{JH}(b_{J}-b_{H})$. A prior
$c_{JH}\sim\mathcal{N}(0.379,0.012^{2})$ is adopted along with flat priors on
$h$. $(a_{i},b_{i})$ are from the Milky Way fits, the LMC fits or a
combination of both. The models of the scatter about the period–luminosity
relation from the fits of the Milky Way and LMC data are not used as the NGC
1559 and NGC 4258 data are multi-epoch mean magnitudes whilst for the Milky
Way only single-epoch data are available. Therefore, a simple constant scatter
about the period–luminosity relation is adopted. In total, nine parameters are
fitted for: two distance moduli $\mu_{j}$ and widths $\sigma_{0,j}$ of the
Gaussian components (no uncertainties are available for the NGC 4258 sample so
the width models both intrinsic and observational spread), their relative
weight ($\vartheta_{1}=1-\vartheta_{2}$), a scaling ($f$) of the reported
uncertainties ($\sigma_{F160W}$), the colour term ($c_{JH}$) and the two
parameters $h=(a,b)$ of the (linear) period–luminosity relation. Logarithmic
priors are used for all intrinsically positive parameters. The condition
$\sigma_{0,1}<\sigma_{0,2}$ is to identify the outlier as the 2nd component.
Further priors are adopted on the mixing parameter
$\ln\vartheta_{1}\sim\mathcal{N}(0,1)$ (with $0<\vartheta_{1}<1$ and
$\vartheta_{1}+\vartheta_{2}=1$) and the error scaling $\ln
f\sim\mathcal{N}(0,1)$. The model is sampled from using emcee (Foreman-Mackey
et al., 2013).

The period–luminosity relation anchors (as defined in equation (13)) are
selected as the O-rich Mira variable period–luminosity relations reported in
Table 1 for the Milky Way and Table 5 for the LMC. Using the Milky Way
relation for the NGC 4258 sample gives $\mu_{4258}=(29.34\pm 0.05)$,
$\mu_{4258}=(29.36\pm 0.04)$ and $\mu_{4258}=(29.36\pm 0.05)$ for the
‘Bronze’, ‘Silver’ and ‘Gold’ samples respectively in good agreement with
(though slightly lower than) the water maser distance of
$\mu_{4258,\mathrm{maser}}=(29.398\pm 0.032)$. Due to the similarity of the
results, from now on the ‘Bronze’ sample is used. The agreement with the water
maser distance suggests the level of C-rich contamination is low in NGC 4258
and the metallicities of the Mira variables in the Milky Way are similar to
those in NGC 4258. If instead the O-rich period–luminosity relation for the
central LMC sample from Yuan et al. (2017b) is used (as given in Table 5), it
is found that $\mu_{4258}=(29.53\pm 0.06)$ which is $\sim 2\sigma$ higher than
the water maser distance modulus suggesting the metallicities of the LMC Mira
variables are lower than those in NGC 4258 (Bresolin, 2011).

Initially the NGC 1559 sample is assumed to be purely composed of O-rich Mira
variables and this assumption is relaxed below. These examples are
illustrative as it is expected that the selection of Mira variables will
introduce C-rich contamination, so the results should not be taken as
realistic estimates of $H_{0}$. In the top section of Table 4 the results
using the Milky Way O-rich relation and the LMC O-rich relation are reported.
For the Milky Way O-rich relation a $\sim
5\,\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$ higher $H_{0}$ is found than when
using the LMC O-rich relation (here the LMC relation for the Yuan et al. 2017b
subsample is used which is $\sim 0.05\,\mathrm{mag}$ fainter than the
relations for the full LMC sample) due to the different period–luminosity
zeropoints (driven by population effects). When combining both in the
modelling, different $\mu_{j}$, $\sigma_{0,j}$, $f$ and $\vartheta_{j}$ are
used for each galaxy, and for NGC 4258 the prior
$\mu_{1}\sim\mathcal{N}(29.398,0.032^{2})$ is used (Reid et al., 2019). For
the pure O-rich period–luminosity relation case, the combined NGC 4258, MW and
LMC fits give an average value of $H_{0}$ between the estimates from the MW
and LMC alone.

###  6.2 C-rich contamination

Figure 10: LMC Mira variable amplitudes: $G$-band Fourier amplitude, $\Delta
G_{\mathrm{Fourier}}$, or amplitude measured from the Gaia DR3 uncertainties,
$\Delta G$ (fainter points), against the approximate $F160W$ amplitude $\Delta
F160W$ (computed from the modelled $J$ and $H$ amplitude scaled by the ratio
of the primary period amplitude to the total amplitude from the OGLE $I$-band
light curves). The blue circles are classified as O-rich whilst the red
squares are C-rich using the classification from Lebzelter et al. (2018). The
colour shades correspond to the period. All stars have periods less than
$400\,\mathrm{days}$. The approximate trends of the two types are shown by the
yellow dashed lines. For the full sample, C-rich Mira variables make up
$61\,\mathrm{per\,cent}$ of the sample, whilst for the selection $0.4<\Delta
F160W<0.8$ they make up $32\,\mathrm{per\,cent}$ or $30\,\mathrm{per\,cent}$
using the approximate $\Delta F160W$ computed from $\Delta G$ (yellow dashed
lines). The vertical dashed line is $\Delta G=\Delta
G_{\mathrm{thresh}}=0.865\,\mathrm{mag}$ which defines the lower boundary for
Mira variables from Gaia photometry.

The first set of models ignored the selection of the Mira variables simply
assuming that the samples were fair representations of the O-rich Mira
population. However, there is significant but uncertain contamination from
C-rich Mira variables in these samples. Limiting to periods less than $400$
days mitigates the impact of C-rich contamination considerably but some
contamination remains that typically makes the mean magnitude at fixed period
fainter, but also flattens the period–luminosity relation and increases the
scatter at longer periods. As C-rich Mira variables have higher near-infrared
amplitudes than O-rich Mira variables in this period range, Huang et al.
(2020) imposed a cut of $0.4<\Delta F160W<0.8$ on their sample to isolate the
O-rich Mira variables. Here $\Delta F160W$ is the peak-to-trough amplitude
over a single cycle and neglects any longer term periodic trends common for
Mira variables. This cut reduces C-rich contamination but some contaminants
remain. First, the impact of this cut on the MW and LMC samples is estimated
and then models of the period–luminosity relation with appropriate C-rich
contamination levels are generated.

The LMC Mira variable sample from Yuan et al. (2017b) have well-sampled I-band
light curves from OGLE for which Soszyński et al. (2009) have provided
amplitudes $I_{1}$ and $I_{2}$ for two identified periods. Yuan et al. (2017b)
used the $I$-band light curves to model the more sparsely sampled $JHK_{s}$
light curves reporting the mean, maximum and minimum $JHK_{s}$ magnitudes. The
single-cycle amplitude in $JHK_{s}$ is approximated as e.g. $\Delta
J=(J_{\mathrm{min}}-J_{\mathrm{max}})I_{1}/(I_{1}+I_{2})$ where the ratio
between the amplitudes of the two periodic trends is assumed to be similar in
all (near-infrared) bands. Using these approximations, the single cycle
$\Delta F160W$ is found as $\Delta H+0.38(\Delta J-\Delta H)$ (see previous
subsection). For the main Milky Way sample of Mira variables from Gaia, only
amplitude indicators in the Gaia passbands are available. Gaia’s observing
window ($22$ months) is relatively short compared to long period trends in
Mira variables so the Gaia amplitudes are assumed to correspond most closely
with single-cycle amplitudes. Fig. 10 displays $\Delta G$ as defined in
equation (1) against $\Delta F160W$ for the LMC sample from Yuan et al.
(2017b). We see $\Delta G$ correlates with $\Delta F160W$ but follows
different relations for O-rich ($\Delta F160W\approx 0.4+0.185(\Delta G-0.7)$)
and C-rich ($\Delta F160W\approx 0.4+0.617(\Delta G-0.63)$) Mira variables. As
shown by Iwanek et al. (2021a), O-rich Mira variables have a steeper fall-off
in amplitude with wavelength than C-rich Mira variables. For the full LMC
sample, $61\,\mathrm{per\,cent}$ of the Mira variables with periods
$<400\,\mathrm{days}$ are C-rich, whilst restricting to $0.4<\Delta F160W<0.8$
(using the infrared amplitudes) reduces this to $32\,\,\mathrm{per\,cent}$.
These numbers are in good agreement with those reported by Huang et al.
(2020). Using the approximate $\Delta F160W$ computed from $\Delta G$,
$30\,\,\mathrm{per\,cent}$ of selected stars are C-rich, thus validating the
approximate relations. Repeating this analysis for the Milky Way Mira sample,
$2.0\,\mathrm{per\,cent}$ of the sample without spatial cuts is C-rich which
reduces to $0.6\,\mathrm{per\,cent}$ using $0.4<\Delta F160W<0.8$, whilst
removing the bulge and $|b|<3\,\mathrm{deg}$ results in a similar reduction
from $3.4\,\mathrm{per\,cent}$ to $0.8\,\mathrm{per\,cent}$. Clearly even with
the $\Delta F160W$ cut, for more metal-poor systems the C-rich contamination
can be significant.

Period–luminosity relations like those in Section 4 and Appendix C have been
fitted to the contaminated LMC and Milky Way Mira variables with $P<400$ days
and $0.4<\Delta F160W<0.8$ (using the previously derived approximate
relations). Instead of a linear relation, a quadratic period–luminosity
relation of the form

$m_{\mathrm{abs,contam}}(P)=a+b(\log_{10}P-2.3)+b_{2}(\log_{10}P-2.3)^{2},$
(14)

is used due to increasing C-rich contamination with increasing period causing
a down-turning for periods $>300$ days (also exhibited by a contaminated SMC
sample). The spread about the period–luminosity relation is modelled as

$\sigma_{\mu,\mathrm{contam}}(P)=\sigma_{2.3}+\begin{cases}0,&\mathrm{if}\log_{10}P<2.3,\\\
p_{\sigma}(\log_{10}P-2.3)\\\
+q_{\sigma}(\log_{10}P-2.3)^{2},&\mathrm{otherwise},\end{cases}$ (15)

to capture the sharp increase in the scatter for periods $>300$ days due to
the C-rich contamination. The results of these fits to the Milky Way stars
with $|b|>3\,\mathrm{deg}$ and outside the bulge region, and to the central
LMC sample of Mira variables with measurements from Yuan et al. (2017b) are
given in Table 3.

The modelling of the previous section has been repeated with these C-rich
contaminated models (using $f=(a,b,b_{2})$) and the results are reported in
Table 4. This modelling assumes the contamination level is the same in NGC
1559 as in the MW and/or LMC. For the MW, the C-rich contamination is so low
that adopting the contaminated model makes a very small difference to $H_{0}$.
However, for the LMC C-rich contamination changes $H_{0}$ by $\sim
10\,\mathrm{km\,s}^{-1}\mathrm{\,Mpc}^{-1}$. This shows that the effects of
C-rich contamination are comparable, if not larger, than the effects of
age/metallicity on the period–luminosity relations. Inclusion of the NGC 4258
measurements (already argued to have a low C-rich contamination level based on
comparison with the water maser distance) brings $H_{0}$ down to values more
consistent with the MW O-rich model or the MW (weakly) C-rich contaminated
model.

Table 3: C-rich contaminated period–luminosity relations. Quadratic relations of the form $m_{\mathrm{abs,contam}}(P)=a+b(\log_{10}P-2.3)+b_{2}(\log_{10}P-2.3)^{2}$ have been fitted to Mira variable samples defined by $0.4<\Delta F160W<0.8$. C/(O+C) is the fraction of C-rich Mira variables in each system (also called $\eta$ in the modelling). The $F160W$ relations are derived using $F160W=H+(0.38\pm 0.01)(J-H)$. System/Band | C/(O+C) | $a$ | $b$ | $b_{2}$
---|---|---|---|---
MW $J$ | $0.008$ | $-5.65\pm 0.02$ | $-3.34\pm 0.08$ | $-1.57\pm 0.40$
MW $H$ | $0.008$ | $-6.45\pm 0.02$ | $-3.61\pm 0.08$ | $-1.52\pm 0.39$
MW $F160W$ | $0.008$ | $-6.15\pm 0.02$ | $-3.51\pm 0.06$ | $-1.55\pm 0.28$
LMC $J$ | $0.324$ | $-5.89\pm 0.03$ | $-2.19\pm 0.15$ | $+5.71\pm 0.97$
LMC $H$ | $0.324$ | $-6.62\pm 0.03$ | $-2.91\pm 0.12$ | $+2.31\pm 0.78$
LMC $F160W$ | $0.324$ | $-6.34\pm 0.02$ | $-2.64\pm 0.09$ | $+3.59\pm 0.62$

###  6.3 Variable C-rich contamination

Table 4: Hubble constant measurements using Mira variables in the SNIa host galaxy NGC 1559 and a range of different anchors. $H_{0}$ in units of $\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$. $\mu_{1559}$ is the distance modulus of NGC 1559, $\mu_{4258}$ is the distance modulus of NGC 4258, $a$ is the zeropoint of the $F160W$ period–luminosity relation evaluated at $\log_{10}P=2.3$, $b$ is the slope of the $F160W$ period–luminosity relation with $\log_{10}P$ and when given $b_{2}$ the quadratic term (i.e the relation is $a+b(\log_{10}P-2.3)+b_{2}(\log_{10}P-2.3)^{2}$). For reference, the distance modulus to the water maser in NGC 4258 is $\mu_{4258,\mathrm{maser}}=(29.398\pm 0.032)\,\mathrm{mag}$ (Reid et al., 2019), the $H_{0}$ measurement from Planck Collaboration et al. (2014) is $(67.4\pm 0.5)\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$, the recent Cepheid-based $H_{0}$ estimate from Riess et al. (2022a) is $(73.04\pm 1.04)\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$ or $(73.01\pm 0.99)\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$ for those Cepheids in clusters (Riess et al., 2022b), the recent tip of the giant branch $H_{0}$ estimate from Freedman (2021) is $69.8\pm 0.6\mathrm{(stat)}\pm 1.6\mathrm{(sys)}\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$ and the combination of Cepheid-based and tip of the giant branch from Riess et al. (2022a) is $(72.53\pm 0.99)\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$. The top section uses pure O-rich Mira variable period–luminosity relations whilst the middle section uses Mira variable period–luminosity relations for contaminated samples. These are illustrative limits and the results should not be considered as recommended measurements. The final section uses a variable C-rich fraction for each system with the recommended final measurement in bold. Anchor | $H_{0}$ | $\mu_{1559}$ | $\mu_{4258}$ | $a$ | $b$ | $b_{2}$
---|---|---|---|---|---|---
MW O-rich | $71.6\pm 4.0$ | $31.43\pm 0.05$ | $-$ | $-6.16\pm 0.04$ | $-3.73\pm 0.10$ | $-$
LMC O-rich | $66.7\pm 3.8$ | $31.58\pm 0.06$ | $-$ | $-6.35\pm 0.06$ | $-3.51\pm 0.15$ | $-$
NGC 4258 + MW +LMC O-rich | $69.9\pm 3.9$ | $31.48\pm 0.05$ | $29.40\pm 0.03$ | $-6.22\pm 0.04$ | $-3.64\pm 0.08$ | $-$
MW C-rich contam. | $71.2\pm 4.1$ | $31.44\pm 0.06$ | $-$ | $-6.15\pm 0.04$ | $-3.51\pm 0.11$ | $-1.51\pm 0.56$
LMC C-rich contam. | $76.5\pm 4.5$ | $31.28\pm 0.07$ | $-$ | $-6.34\pm 0.06$ | $-2.67\pm 0.18$ | $3.28\pm 1.09$
NGC 4258 + MW + LMC C-rich contam. | $70.4\pm 4.0$ | $31.47\pm 0.05$ | $29.39\pm 0.03$ | $-6.22\pm 0.04$ | $-3.36\pm 0.09$ | $-1.19\pm 0.45$
MW + LMC variable contam. | $74.6\pm 4.5$ | $31.34\pm 0.07$ | $-$ | $-6.14\pm 0.05$ | $-3.54\pm 0.11$ | $-1.72\pm 0.57$
NGC 4258 + MW + LMC variable contam. | $73.7\pm 4.4$ | $31.37\pm 0.07$ | $29.38\pm 0.03$ | $-6.15\pm 0.04$ | $-3.54\pm 0.12$ | $-1.69\pm 0.58$
Figure 11: Period–magnitude diagrams for Mira variables in NGC 1559 (left) and
NGC 4258 (right – smaller faint points are non-‘Gold’ Mira variables). The
black line shows the best-fitting period–luminosity relation (as given in the
left panel) shifted by the distance estimate for each galaxy and evaluated at
the estimated C-rich contamination fraction (as given in each panel) with the
grey interval showing the $1\sigma$ uncertainty. The period–luminosity
relation has been anchored using the Gaia EDR3 Milky Way results from this
paper together with the LMC Mira variables anchored with the detached
eclipsing binary distance from Pietrzyński et al. (2019) and NGC 4258 Mira
variables anchored using the distance to the NGC 4258 maser from Reid et al.
(2019). The resulting Hubble constant using SN 2005df is $(73.7\pm
4.4)\,\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$.

In reality, the C-rich contamination of the NGC 4258 and NGC 1559 sample is
unknown and the best choice of period–luminosity relation lies somewhere
between the pure O-rich and C-rich contaminated cases. To incorporate this
into the modelling, the NGC 4258 and NGC 1559 samples are considered to have
individual unknown C-rich contamination factors, $\eta=N_{C}/(N_{C}+N_{O})$,
where $N_{i}$ is the number of $i$-rich stars. The parameters of the quadratic
period–luminosity relations, $a,b$ and $b_{2}$, are assumed to vary linearly
with $\eta$. The Milky Way and LMC C-rich contaminated samples then give two
(probabilistic) points on this relation with $\eta_{\mathrm{MW}}=0.008$ and
$\eta_{\mathrm{LMC}}=0.324$ respectively (the binomial uncertainty in these
numbers is not considered). In this way, the range of different environments
has been reduced down to a single parameter, $\eta$. As discussed in Section
5.2, at the most basic level both age and metallicity play a role in
determining the period–luminosity relation so a single parameter is an
oversimplification. Both age and metallicity also govern the degree of C-rich
contamination. Age largely influences the periods of Mira variables and so
over a limited period range $\eta$ can be considered as a proxy for
metallicity (Brewer et al., 1995; Hamren et al., 2015) and age is then
considered as affecting the period distribution rather than the shape of the
period–luminosity relation. As $\eta$ is increased (metallicity decreased),
the period–luminosity relation of the O-rich Mira variables becomes brighter,
the C-rich contamination increases and the contaminated period–luminosity
relation gets fainter and flatter (as C-rich Mira variables are more common
for younger, longer period populations). In the absence of clear O-rich/C-rich
discrimination on a star-by-star basis, C-rich contamination can be measured
from the shape of the period–luminosity relation, and in turn the metallicity
of the environment measured and the period–luminosity zeropoint more precisely
known.

This approach is incorporated into the modelling by taking $h=(a,b,b_{2})$ as
the period–luminosity relation parameters for a completely uncontaminated
sample, $\eta=0$, and introducing parameters of the gradients of $h$ with
$\eta$, $g=\mathrm{d}h/\mathrm{d}\eta$. The term in equation (13) is then
adjusted to instead be

$p(\bar{h}_{i}|\eta_{i},h,g,c_{JH})=-\tfrac{1}{2}(\bar{h}_{i}-\eta_{i}g-h)^{\mathrm{T}}\Sigma_{hi}^{-1}(\bar{h}_{i}-\eta_{i}g-h).$
(16)

for system $i$ (Milky Way and LMC) with measured $\bar{h}_{i}=\langle
h_{H}+c_{JH}(h_{J}-h_{H})\rangle$ and covariance $\Sigma_{hi}$ as before. For
the NGC 1559 Mira variables, the period–luminosity relation parameters are
$f+\eta_{1559}g$ (and similar for NGC 4258). The prior
$\eta\sim\mathcal{N}(0,0.3^{2})$ is used as it is anticipated the C-rich
contamination in both NGC 1559 and NGC 4258 is lower than that in the LMC. In
theory, the scatter parameters $\sigma_{0,j}$ could also be made to vary with
$\eta$ in a similar way. This may give more handle on $\eta$ as the scatter
about the period–luminosity relation increases substantially with increased
C-rich contamination. However, mean $F160W$ magnitudes for NGC 4258 and NGC
1559 are considered, whilst only single-epoch data are available for the Milky
Way sample making it awkward to estimate the expected mean scatter at fixed
contamination.

This procedure is quite similar to that of Huang et al. (2020), who used a
fixed gradient linear model fitted to the mean magnitude in a set of period
bins to find the variation in the zeropoint. The zeropoint variation was then
matched on to the corresponding variation in the LMC scaling by an unknown
contamination factor, $\alpha$, that gave the fraction of the zeropoint shift
between the contaminated and uncontaminated LMC sample that must be applied to
the sample.

The variable C-rich contamination models are fitted to NGC 1559 alone, and NGC
1559 and NGC 4258 together and the results are reported in Table 4. There is
good agreement between the two models. In agreement with the previous O-rich
models of NGC 4258, it is found that $\eta_{4258}<0.13$ at
$95\,\mathrm{per\,cent}$ confidence. For NGC 1559, it is found that
$\eta_{1559}=(0.17\pm 0.07)$ or $\eta_{1559}/\eta_{\mathrm{LMC}}=(0.52\pm
0.22)$, in very good agreement with Huang et al. (2020) who find the
adjustment of the zeropoint must be a fraction $(0.58\pm 0.18)$ of the LMC
zeropoint adjustment. The final estimate of $H_{0}=(73.7\pm
4.4)\,\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$ is in agreement with the analysis
of Huang et al. (2020) who found $H_{0}=(73.3\pm
4.0)\,\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$. Despite also using the Milky Way
Mira variables in the analysis, the uncertainty here is
$0.4\,\mathrm{km\,s}^{-1}\mathrm{Mpc}^{-1}$ larger. This is probably because
of the difference in the handling of the C-rich contamination and the
marginalization over the period–luminosity relation gradient. If instead the
parallax zeropoint corrections from Groenewegen (2021) or Maíz Apellániz
(2022) as described in Section 4.4 are used, larger $H_{0}$ of $(76.3\pm 4.3)$
and $(75.8\pm 4.2)\,\mathrm{km\,s}^{-1}\mathrm{kpc}^{-1}$ respectively are
found. However, these zeropoint models are disfavoured as they lead to larger
parallax residuals as a function of $G$. Nonetheless, these estimates are
within the reported uncertainties, and point to the importance of an improved
understanding of the Gaia parallax systematics for refining these estimates.
In accord with many recent estimates of the Hubble constant based on near-
Universe tracers, the estimate is higher than the Planck Collaboration et al.
(2014) estimate but only at the $\sim 1.4\sigma$ level. Although some
uncertainty arises from the modelling of the Mira variable period–luminosity
relation and contamination, the dominant uncertainty is from the peak
luminosity of SN2005df so further measurements of Mira variables in Type Ia
host galaxies are required.

##  7 Conclusions

Preliminary period–luminosity calibrations have been presented for Milky Way
O-rich Mira variables in the 2MASS $JHK_{s}$ bands using astrometric data from
Gaia Data Release 3. The relations have been derived using a probabilistic
model incorporating a flexible distance prior and models for the Gaia EDR3
parallax zeropoint and uncertainty underestimates. Period–luminosity relations
have been estimated for $JHK_{s}$ magnitudes and extinction-free Wesenheit
indices, and also estimated for the C-rich Mira population. The corresponding
relations for the Large and Small Magellanic Clouds, and the Sagittarius dwarf
spheroidal galaxy have also been estimated.

The Mira variables provide an interesting regime for testing the Gaia
astrometry. Although the large angular size is a potential concern for a
handful of the nearest brightest Mira variables, it appears the Gaia EDR3
astrometry is accurate for these very red, bright stars. A full investigation
of the theoretical expected astrometric performance of Gaia EDR3 for AGB stars
has been performed indicating that, despite the intrinsic photocentre wobble
of these stars, on average the measured parallaxes should be unbiased but the
uncertainties are likely underestimated for $\varpi\gtrsim 0.5\,\mathrm{mas}$
and $G\lesssim 11$. It is found that the parallax zeropoint corrections from