# A VLBA Trigonometric Parallax for RR Aql and the Mira PL Relation

Yan Sun Shanghai Astronomical Observatory, Chinese Academy of Sciences
80 Nandan Road, Shanghai 200030, China  University of Chinese Academy of
Sciences
No.19 (A) Yuquan Rd. Shijingshan, Beijing, 100049, China  Bo Zhang Shanghai
Astronomical Observatory, Chinese Academy of Sciences
80 Nandan Road, Shanghai 200030, China  Mark J. Reid  Center for Astrophysics
$|$ Harvard & Smithsonian
60 Garden Street, Cambridge, MA 02138, USA  Shuangjing Xu  Shanghai
Astronomical Observatory, Chinese Academy of Sciences
80 Nandan Road, Shanghai 200030, China  Korea Astronomy and Space Science
Institute, 776 Daedeok-daero, Yuseong-gu, Daejeon 34055, Republic of Korea
Shiming Wen  Shanghai Astronomical Observatory, Chinese Academy of Sciences
80 Nandan Road, Shanghai 200030, China  Jingdong Zhang  Shanghai Astronomical
Observatory, Chinese Academy of Sciences
80 Nandan Road, Shanghai 200030, China  University of Chinese Academy of
Sciences
No.19 (A) Yuquan Rd. Shijingshan, Beijing, 100049, China  Xingwu Zheng  School
of Astronomy and Space Science, Nanjing University
22 Hankou Road, Nanjing 210093, China

(Received December 24, 2021; Accepted April 22, 2022)

###### Abstract

We report VLBA observations of 22 GHz H2O and 43 GHz SiO masers toward the
Mira variable RR Aql. By fitting the SiO maser emission to a circular ring, we
estimate the absolute stellar position of RR Aql and find agreement with
$Gaia$ astrometry to within the joint uncertainty of $\approx 1$ mas. Using
the maser astrometry we measure a stellar parallax of 2.44 $\pm$ 0.07 mas,
corresponding to a distance of 410${}^{+12}_{-11}$ pc. The maser parallax
deviates significantly from the $Gaia$ EDR3 parallax of 1.95 $\pm$ 0.11 mas,
indicating a $3.8\sigma$ tension between radio and optical measurements. This
tension is most likely caused by optical photo-center variations limiting the
Gaia astrometric accuracy for this Mira variable. Combining infrared
magnitudes with parallaxes for RR Aql and other Miras, we fit a period-
luminosity relation using a Bayesian approach with MCMC sampling and a strong
prior for the slope of –3.60 $\pm$ 0.30 from the LMC. We find a $K$-band zero-
point (defined at logP(days) = 2.30) of –6.79 $\pm$ 0.15 mag using VLBI
parallaxes and –7.08 $\pm$ 0.29 mag using $Gaia$ parallaxes. The $Gaia$ zero-
point is statistically consistent with the more accurate VLBI value.

astrometry – masers – stars: individual (RR Aql) – stars: Mira-type –
techniques: high angular resolution

††journal: ApJ††facilities: VLBA††software: DiFX (Deller et al., 2007), AIPS
(van Moorsel et al., 1996)

##  1 Introduction

Mira variables in the Large Magellanic Could (LMC) show a well defined period-
luminosity relation (PLR) in the infrared (Glass & Lloyd Evans, 1981;
Whitelock et al., 2008; Yuan et al., 2017; Iwanek et al., 2021a). Assuming the
PLR slope from the LMC, zero-points of the PLRs for Miras have also been
determined in other galaxies, e.g., M 33 (Yuan et al., 2018), NGC 4258 (Huang
et al., 2018) and NGC 1559 (Huang et al., 2020). Compared with Cepheids, Miras
are brighter at infrared wavelengths (Feast & Whitelock, 2014) and have a
larger range of ages and greater numbers (Eyer & Cuypers, 2000). Miras can
also provide distances rivaling those of Cepheids (Iwanek et al., 2021a).
Thus, Miras will be important in the era of large space infrared telescopes,
e.g., the James Webb Space Telescope (JWST). However, there are some issues
regarding the use of the PLR of Miras which need to be addressed. First, more
observational evidence that the PLR of Miras in different galaxies is
insensitive to metallicity variations is desirable. Second, if the PLR of
Miras is universal, the zero-point for the PLR should be better determined.
Third, using the PLR of LMC Miras to estimate distances to other galaxies
includes a component of systematic error due to the assumed distance to the
LMC; currently the most accurate distance modulus of $18.477\pm 0.026$ mag was
published by Pietrzyński et al. (2019).

The $Gaia$ Early Data Release 3 (EDR3) provides parallaxes for more than 1000
Miras with formal uncertainties of $\leq$ 0.5 mas. This provides an
opportunity to derive a Galactic PLR for Miras independently of the LMC. Thus,
it is important to check the $Gaia$ parallaxes with independent methods of
comparable or better accuracy, e.g., using Very Long Baseline Interferometry
(VLBI). As shown in Xu et al. (2019), there are significant discrepancies
between the parallaxes from $Gaia$ and VLBI for Miras. One likely reason for
these discrepancies involves optical photo-center variations, owing to a small
number of very large convective cells (Chiavassa et al., 2011) or hot spots
(Lacour et al., 2009). Since Miras typically have diameters of $\approx 2$ AU
(Willson, 1981), which imply angular sizes of twice their parallaxes, photo-
center variations of, say, 10% could lead to astrometric noise of 20% of a
parallax. To further complicate this problem, Miras typically have periods
near one year (Willson, 1981), and thus photocenter variations can easily
correlate with a yearly parallax signature.

In order to determine the PLR for Galactic Miras, accurate parallaxes for 5
Miras have been obtained with Very Long Baseline Array (VLBA) observations and
17 have been obtained with the VLBI Exploration of Radio Astrometry (VERA), as
part of a key project started in 2004 (Nakagawa, 2009; VERA Collaboration et
al., 2020). Since the total number of VLBI parallaxes for Miras is modest,
adding more is important. RR Aql is an oxygen rich Mira with a period of 396
days with strong circumstellar OH, H2O and SiO masers (Benson et al., 1990).
Although its trigonometric parallax has been estimated optically by
$Hipparcos$ (van Leeuwen, 2007) and in the radio by observations of its OH
masers (Vlemmings & van Langevelde, 2007), the uncertainties in these
measurements are too large to be useful as a check on $Gaia$. In this paper,
we present VLBI parallax measurements of 22 GHz H2O and 43 GHz SiO maser
emission toward RR Aql using the VLBA. In Section 2, we describe the
observations and analysis procedures. In Section 3, we present the positions
of both H2O and SiO maser spots in RR Aql relative to extragalactic quasars at
different epochs and the derived parallaxes and proper motions. In Section 4,
we compare the absolute positions of $Gaia$ and our VLBI results. Then in
Section 5, we fit a PLR using the VLBI and $Gaia$ parallaxes separately, using
a Bayesian approach with MCMC sampling which avoids hand editing of data by
decreasing the penalty in the likelihood function for outliers compared to
using a Gaussian data uncertainty. Finally, we summarize our studies and
discuss the future outlook for the PLR of Miras in Section 6.

##  2 Observation and Data reduction

We conducted phase-referenced observations of RR Aql ($\alpha_{\rm J2000}$ =
19h57m36${\hbox to0.0pt{.\hss}}^{\rm s}$0334/36${\hbox to0.0pt{.\hss}}^{\rm
s}$0377, $\delta_{\rm J2000}$ = –01°53′12${\hbox
to0.0pt{.\hss}}{{}^{\prime\prime}}$157/12${\hbox
to0.0pt{.\hss}}{{}^{\prime\prime}}$203 for H2O/SiO masers, see Table 1 for
details.) relative to a background quasar J1947–0103 ($\alpha_{\rm J2000}$ =
19h47m43${\hbox to0.0pt{.\hss}}^{\rm s}$7837, $\alpha_{\rm J2000}$ =
–01°03′24${\hbox to0.0pt{.\hss}}{{}^{\prime\prime}}$528) at 22 and 43 GHz
under the National Radio Astronomy Observatory’s111The National Radio
Astronomy Observatory is a facility of the National Science Foundation
operated under cooperative agreement by Associated Universities, Inc. VLBA
program BZ069. These observations included six epochs spanning about one year
from 2017/10/29 to 2018/11/02 at intervals of roughly two months. An
observation on 2018/08/21 failed with no detections of any quasars or masers
on all of the interferometer baselines; the reason for the failure is unknown.
Each epoch used a 6-h track which included multiple observing blocks of 25 min
of fast-switching between RR Aql and a background quasar J1947–0103222Another
quasar, J2001–0042, was also observed, but proved too weak to use.. Scan
durations for the maser source were 60 and 40 seconds, yielding typical on-
source times of 40 and 20 seconds at 22 and 43 GHz, respectively. For better
interferometer uv-coverage at each band, we alternated the fast-switching
blocks at 22 and 43 GHz. The data were taken in left- and right-circular
polarizations with four adjacent bands of 16 MHz bandwidth, with the maser
emission placed in the third band.

We observed three International Celestial Reference Frame (ICRF) sources:
3C345, 1638+572 and 1740+521 at the beginning, middle and end of the
observations, in order to calibrate electronic delay and phase offsets among
different bands. In order to calibrate tropospheric delays, we adopted the
methods described in Reid et al. (2009). We included three 0.5-h “geodetic
blocks” at the beginning, middle and end of tracks. Geodetic blocks were
recorded in left-circular polarization with eight 16-MHz bands spanning 480
MHz between 23.5 and 24.0 GHz. The eight bands were spaced in a “minimum
redundancy configuration” to uniformly sample frequency differences.

The data were correlated in two passes using the DiFX software correlator
(Deller et al., 2007) in Socorro, NM. One pass generated 32 spectral channels
for all the data, and another pass generated 1000 spectral channels (in “zoom”
mode) only for the single (dual-polarized) frequency band containing the
masers, yielding a spectral channels spaced by 8 kHz, corresponding to
velocity resolutions of 0.11 and 0.06 km s-1 for H2O and SiO maser emissions,
respectively.

Data calibration and imaging were performed using the Astronomical Image
Processing System (AIPS; van Moorsel et al. 1996). For phase calibration, we
corrected for the effects of diurnal feed rotation, errors in the Earth
Orientation Parameters (EOP) and the positions of the masers and the
background source adopted in the correlation. We removed instrumental delays
using the data from 3C345 and calibrated ionospheric and tropospheric delays
using the Global Ionospheric Model, provided by the International GNSS
Service, and data from the geodetic blocks (Reid & Honma, 2014). For amplitude
calibration, we adopted the VLBA’s recommended system equivalent flux
densities (SEFD) for each antenna and feed, using gain curves and system
temperatures.

Figure 1: Vector-averged interferometer spectra of H2O (left panel) and SiO
(right panel) maser emission toward RR Aql from the Kitt Peak to Los Alamos
baseline at five epochs.

As shown in Figure 1, we found strong H2O and SiO maser emissions with peaks
at the velocity of Local Standard Rest (VLSR) values of around 29 and 31 km
s-1 for all epochs. In this paper, we measure parallax using the maser
channels centered at 28.8 and 28.9 km s-1 for H2O masers and 30.9 km s-1 for
the SiO masers. Since the water maser emissions were heavily resolved on the
longer baselines, we phase referenced those data using the quasar J1947–0103
and only used the data from the VLBA inner-five antennas: Fort Davis (FD),
Kitt Peak (KP), Los Alamos (LA), Owens Valley (OV) and Pie Town (PT). These
antennas are located close to the center of the array, with relatively short
baselines. However, on 2017/10/29, there were clock jumps at the PT antenna,
which were evident in plots of interferometric delay versus time using the
strongest calibrator, and we discarded the data involving PT. In addition, the
PT antenna was unavailable for the observation on 2017/12/31 owing to a
maintenance issue. For the SiO masers, we used the 30.9 km s-1 maser spot as
the interferometer phase reference (i.e., “inverse” phase-referencing).

After calibration, we combined the four dual-polarized bands and imaged the
continuum emission of the background source using the AIPS task IMAGR, which
Fourier transforms the visibility data to make a brightness image. The
restoring beam for each epoch was calculated using uniform weighting. In
individual spectral channels, we estimated peak brightnesses, positions and
angular sizes using the task JMFIT, which fits Gaussian brightness
distributions to a portion of an image. Table 1 lists basic information for
the observed sources, and Table 2 gives the measured residual position offsets
for the maser spots relative to the quasar J1947–0103 used for parallax
fitting. Figure 3 shows an example image of the H2O maser spots at VLSR = 28.8
km s-1 and 28.9 km s-1 toward RR Aql and the background source J1947–0103.

Table 1: Positions and Brightnesses

Source | R.A. (J2000) (K / Q) | Dec. (J2000)(K / Q) |  $S_{p}$ (K / Q) | $\theta_{sep}$ | P.A. | VLSR (K / Q) | Beam (K / Q)
---|---|---|---|---|---|---|---
Name | (h m s) | (° ′ ″) | (Jy/beam) | (°) | (°) | (km s-1) | (mas mas °)
RR Aql | 19 57 36.0334 / 36.0377 | –01 53 12.157 / 12.203 | 6 / 13 | … | … | 28.9 / 30.9 | 3.49$\times$1.3 @ $+$22 / 1.85$\times$0.6 @ $+$23
J1947–0103 | 19 47 43.7837 | –01 03 24.528 | 0.167 / 0.042 | 2.60 | +159 | … | 3.48$\times$1.2 @ $+$22 / 1.74$\times$0.6 @ $+$23

Note. — K and Q denote the observing bands, $\sim 22$ GHz for the H2O maser
and $\sim 43$ GHz for the SiO maser, respectively. The second and third
columns list the absolute positions of the H2O maser spot at VLSR = 28.9 km
s-1 and of the SiO maser spot at VLSR = 30.9 km s-1 on 2018/01/09. The fourth
and seventh columns give the peak brightnesses ($S_{p}$) and VLSR of the maser
spots. The fifth and sixth columns give the separations ($\theta_{sep})$ and
position angles (P.A.) east of north of the quasar relative to the maser. The
last column gives the FWHM size and P.A. of the Gaussian restoring beam at the
first epoch. The quasar position is from http://astrogeo.org.

Table 2: Residual position differences between maser spots and J1947–0103 used
for the parallax fits.

Epoch | Maser | VLSR | East offset | North offset
---|---|---|---|---
(yyyy/mm/dd) | Species | (km s-1) | (mas) | (mas)
2017/10/29 | H2O | 28.9 | 0.043$\pm$0.080 | -0.051$\pm$0.133
2017/12/31 |  |  | –2.242$\pm$0.060 | –8.407$\pm$0.102
2018/03/29 |  |  | –4.543$\pm$0.129 | –17.096$\pm$0.180
2018/05/21 |  |  | –8.113$\pm$0.111 | –24.495$\pm$0.194
2018/11/02 |  |  | –22.480$\pm$0.186 | –46.289$\pm$0.238
2017/10/29 | H2O | 28.8 | –1.067$\pm$0.061 | –6.565$\pm$0.113
2017/12/31 |  |  | –3.112$\pm$0.052 | –15.165$\pm$0.094
2018/03/29 |  |  | –5.850$\pm$0.138 | –24.830$\pm$0.215
2018/05/21 |  |  | –9.058$\pm$0.131 | –31.472$\pm$0.238
2018/11/02 |  |  | –23.535$\pm$0.199 | –53.031$\pm$0.242
2017/10/29 | SiO | 30.9 | 23.945$\pm$0.038 | 50.273$\pm$0.069
2017/12/31 |  |  | 21.800$\pm$0.054 | 41.297$\pm$0.097
2018/03/29 |  |  | 19.286$\pm$0.052 | 30.119$\pm$0.107
2018/05/21 |  |  | 15.468$\pm$0.132 | 26.109$\pm$0.243
2018/11/02 |  |  | 1.125 $\pm$0.069 | 0.736$\pm$0.153

Note. — The first column lists the observing epochs. The second and third
column list the maser species and their VLSR. The fourth and fifth columns
give position offsets relative to $\alpha_{\rm J2000}$ = 19h57m36${\hbox
to0.0pt{.\hss}}^{\rm s}$0337, $\delta_{\rm J2000}$ = –01°53′12${\hbox
to0.0pt{.\hss}}{{}^{\prime\prime}}$148 for the H2O masers on 2017/10/29, and
$\alpha_{\rm J2000}$ = 19h57m36${\hbox to0.0pt{.\hss}}^{\rm s}$0364,
$\delta_{\rm J2000}$ = –01°53′12${\hbox to0.0pt{.\hss}}{{}^{\prime\prime}}$242
for the SiO maser on 2018/11/02. Uncertainties given are statistical only,
based on thermal noise in the images, and do not include systematic effects
such as uncompensated atmospheric delays.

Figure 2:  Figure 3:  Images of the H2O maser spots at VLSR=28.8 km s-1(left
panel) and VLSR=28.9 km s-1 (middle panel) toward RR Aql and the background
radio source J1947–0103 (right panel) on 2018/03/29 using the inner-five
antennas of the VLBA. The restoring beams are shown in the lower left corner
of each panel. Contour levels are spaced linearly by 0.25 Jy beam-1 for RR Aql
and 0.030 Jy beam-1 for J1947–0103.

##  3 Results

###  3.1 Maser spots identification for astrometry

####  3.1.1 Water masers

As shown in Figure 4, the phase variations of the reference source J1947–0103
on 2018/03/29, as an example, are relatively slow; thus the phases on the
baselines of inner-five VLBA antennas can be confidently connected without
2$\pi$-ambiguities between adjacent scans separated by 2 min.

Figure 4: Visibility phase versus observing time of calibrator J1947–0103 at
22 GHz on 2018/03/29. The two-letter station codes of the VLBA inner-five
antennas are labeled on the top-right corner in each panel. FD, KP, LA, OV and
PT indicate Fort Davis, Kitt Peak, Los Alamos, Owens Valley and Pie Town,
respectively. The reference station was KP. The interval between adjacent data
points is 2 min. Figure 5:  Spatial distribution of the H2O maser emission at
VLSR = 28.9 km s-1 toward RR Aql from VLBA observations at five epochs. The
observing epoch is labeled next to each maser spot. Contour levels are integer
multiples of 10% of the peak brightness. Note that on 2018/03/29 two maser
features are detected.

Water maser emission was detected at VLSR of 27.9 km s-1 from 2017/10/29 to
2018/03/29. However, by 2018/05/21 its brightness fell below 0.4 Jy beam-1 and
was not detected. This is not unusual for maser emission from Miras, which
also vary over the stellar period. For parallax measurement, we focus on
analyzing the stronger water maser spots at VLSR of 28.8 and 28.9 km s-1,
which were clearly detected at all epochs. Figure 5 shows the evolution of the
28.9 km s-1 maser spot relative to the calibrator J1947–0103. On 2018/03/29 we
detected strong emission at (E,N) of $\approx(-5,-20)$ mas (left panel) and
weaker emission at $\approx(80,-55)$ mas (right panel). Clearly the weaker
spot is unrelated to those tracking smoothly toward the south-south-west over
the five epochs. Therefore, we excluded the weaker spot from our parallax
fitting. On 2018/11/02 the maser emission displays a head-tail like structure.
This is probably due to blending two maser spots separated by several mas.
Parallax fitting with one or the other spot position clearly shows that the
stronger “head” spot is the continuation of the previous four epochs.

####  3.1.2 SiO masers

As shown in the right panel of Figure 1, the SiO masers have more features
than the H2O masers: there are at least a handful of SiO features at VLSR
ranging from 23 to 31 km s-1. After phase-referencing to the maser spot at
VLSR of 30.9 km s-1, we made the velocity integrated images shown in Figure 6.
There are clear indications of a clumpy ring-like structure at the first three
epochs; in later epochs the western part of the ring weakens and vanishes.
Fortunately, maser feature A was detected at all epochs, and we used this
feature to determine the parallax of RR Aql.

Figure 6:  Spatial distributions of SiO maser spots toward RR Aql from VLBA
observations at five epochs. Maser features are designated with letters. The
long-lived feature A was used for parallax estimation. The dashed circles
denote rings fitted to the SiO maser emissions.  Figure 7:  East (upper panel)
and north (lower panel) offsets relative to J1907–0103 versus time used for
parallax fitting for two H2O maser spots (blue squares and red circles) and
one SiO maser spot (purple triangle), with individual constant offset and
proper motion removed. The VLSR of each maser spot is labeled on the left-top
corner in the upper panel. The best-fitting model (dotted lines) denote the
common parallax curve. Positions are slightly offset in time for clarity.
Table 3: Parallax and proper motion for RR Aql

Method | VLSR | Parallax | $\mu_{x}$ | $\mu_{y}$ | $\Delta x$ | $\Delta y$ | $\chi^{2}_{\nu}$ | $\sigma_{x}$ | $\sigma_{y}$
---|---|---|---|---|---|---|---|---|---
| (km s-1) | (mas) | (mas yr-1) | (mas yr-1) | (mas) | (mas) |  | (mas) | (mas)
H2O | 28.9 | 2.44 $\pm$ 0.07 | –22.29 $\pm$ 0.17 | –45.70 $\pm$ 0.85 | –8.94 $\pm$ 0.05 | –22.43 $\pm$ 0.30 | 0.83 | 0.00 | 0.70
H2O | 28.8 | 2.41 $\pm$ 0.09 | –22.29 $\pm$ 0.23 | –45.91 $\pm$ 0.46 | –9.96 $\pm$ 0.08 | –29.39 $\pm$ 0.17 | 0.96 | 0.09 | 0.31
SiO | 30.9 | 2.44 $\pm$ 0.08 | –22.68 $\pm$ 0.18 | –48.47 $\pm$ 1.46 | 14.86 $\pm$ 0.07 | 26.33 $\pm$ 0.52 | 0.98 | 0.13 | 1.14
Combined |  | 2.44 $\pm$ 0.07 |  |  |  |  | 0.98 | 0.00 | 0.00
$Gaia$ EDR3 |  | 1.95 $\pm$ 0.11 | –20.42 $\pm$ 0.12 | –47.69 $\pm$ 0.08 |  |  |  |  |

Note. — Absolute proper motions are defined as $\mu_{x}$ =
$\mu_{\alpha}\cos{\delta}$ and $\mu_{y}$ = $\mu_{\delta}$. Position offsets
$\Delta x$ and $\Delta y$ are for 2018/05/01 (the midpoint between the first
and last epochs) and are relative to the positions listed in Table 1.
$\chi^{2}_{\nu}$ is the reduced $\chi^{2}$ of post-fit residuals, $\sigma_{x}$
and $\sigma_{y}$ are error floors in $x$ and $y$, respectively. ‘H2O’ and
‘SiO’ denote the maser species; ‘Combined’ means the parameters derived by
combing the data of the two H2O and one SiO maser spots. The last row lists
the parameters from $Gaia$ EDR3.

###  3.2 Parallax and Proper Motion

As described in Section 3.1, we used two H2O and one SiO maser spot detected
in all five epochs spanning one year. The positions of maser spots relative to
the background source were modeled by the parallax sinusoid in both
coordinates (determined by one parameter) and a linear proper motion in each
coordinate, and best parameter values were obtained by variance-weighted
least-squares fitting. The formal position uncertainties listed in Table 2
reflect thermal (random) noise in the interferometric images and do not
include systematic uncertainties. To estimate and allow for systematic
position errors caused by residual atmospheric delays and variable structure
of the maser spots and/or the background quasar, we added “systematic errors”
(error floors) for each coordinate in quadrature with the formal measurement
errors and adjusted these to reach a reduced $\chi^{2}$ for the post-fit
residuals near unity for each coordinate.

Figure 7 shows the positions of the two H2O spots and one SiO spot relative to
the background source J1947–0103, with superposed curves showing the best fit
parallax sinusoid. Note that the error bars in the bottom panel are larger
than those in the top panel, indicating the “error floors” for the north-south
offsets were larger than those for the east-west offsets. This occurs because
the interferometer beams are larger in the north-south direction and
uncompensated atmospheric delays are more correlated with them. Table 3 lists
the individual fits (for each maser spot) and a combined fit (using all three
spots), allowing for different offsets and proper motions among them.

For the combined fits using both H2O and SiO masers, to avoid a single error
floor down-weighting the better H2O maser data, we combined measurement
uncertainties with the error floors used in the individual fits (multiplied by
$0.8$) in revised data files, in order to get a reduced chi-square near unity.
The parallax fitting results are in good agreement within their joint
measurement uncertainties. Considering that the dominant error is from
uncompensated atmospheric delays, which are essentially the same for all spots
and both maser species, we conservatively inflated the formal parallax
uncertainty by $\sqrt{3}$ (for the 3 spots combined). We adopt the parallax
from the combined fits, which is 2.44 $\pm$ 0.07 mas, corresponding to a
distance of 410${}^{+12}_{-11}$ pc. We find that the maser parallax deviates
significantly from the $Gaia$ EDR3 parallax of 1.95 $\pm$ 0.11 mas, indicating
a $3.8\sigma$ tension between radio and optical astrometry. We discuss
possible reasons for this discrepancy in the Section 4.

##  4 Stellar Position of RR Aql

###  4.1 Relative positions of masers respect to the central star of RR Aql

The absolute stellar positions of Miras determined from VLBI observations can
be used to check the optical $Gaia$ positions, which might be corrupted by
non-uniform circumstellar dust extinction and photo-center variation.
Additionally, knowing the position of masers relative to the central stars can
be crucial to understanding gas kinematics in the circumstellar material.

As shown in Figure 6, the sky distribution of SiO masers from our VLBA
observations at the first three epochs show a clumpy ring-like structure. This
has previously been noted for H2O maser emission (e.g., from W Hya by Reid &
Menten 1991) and SiO maser emission (e.g., from TX Cam by Kemball & Diamond
1997). We fit the central position (assumed to be the stellar position) and
the radius of a ring to the maser positions using a least-square method
described by Yoon et al. (2018), which assumes the distribution of SiO masers
lie at a common radius from the central star. The fitted results are shown in
Table 4.

Figure 8 shows relative positions of masers respect to the central star on
2017/10/29. Adopting a stellar radius ($R_{\ast}$) of 3.8 $\pm$ 0.3 mas for RR
Aql from multi-epoch mid-infrared interferometric observations by Karovicova
et al. (2011), we conclude that the SiO masers are distributed near
2$R_{\ast}$, while the H2O maser emission at VLSR= 28.8 km s-1 and 28.9 km s-1
is located $\approx 18$$R_{\ast}$ north-west of the central star. This is
consistent with previously reported circumstellar maser distributions by Reid
& Menten (1997) and Richards et al. (2010).

Table 4: Center position and radius of SiO maser ring

Epoch | $\alpha_{\rm J2000}$ | $\sigma_{\alpha}$ | $\delta_{\rm J2000}$ | $\sigma_{\delta}$ | Radius
---|---|---|---|---|---
(yyyy/mm/dd) | (h m s) | (mas) | (° ′ ″) | (mas) | (mas)
2017/10/29 | 19 57 36.03728 |  $\pm$0.3 | –01 53 12.1908 |  $\pm$1.2 | 10.0$\pm$0.1
2017/12/31 | 19 57 36.03715 |  $\pm$0.3 | –01 53 12.1997 |  $\pm$1.2 | 9.9$\pm$0.1
2018/03/29 | 19 57 36.03797 |  $\pm$0.3 | –01 53 12.2100 |  $\pm$1.2 | 10.2$\pm$0.2

Note. — Columns 2 though 5 list the absolute positions and uncertainties in
Right Ascension and Declination of the center of the SiO maser ring. Column 6
list the best-fit radii of SiO maser rings.

Figure 8:  The sky distribution SiO and H2O masers with respect to the central
star of RR Aql on 2017/10/29. The dotted line represents the fitted SiO maser
ring. The circles in dash-dotted lines indicate multiples of five stellar
radii.

###  4.2 Comparison of the absolute position of RR Aql from VLBA and $Gaia$

In order to evaluate differences between optical photo-centers and radio
measurements, we compare the absolute positions of the center of the SiO ring
with stellar positions from $Gaia$ (see Table 4). At the first three epochs
the SiO masers have a well-constrained ring structure, which give very
accurate positions for the central star. Table 5 shows the differences between
our radio and the $Gaia$ optical stellar positions (using $Gaia$’s five
astrometric parameters listed in EDR3); joint uncertainties are estimated by
adding the individual values in quadrature. There are no statistically
significant differences in either the Right Ascension or Declination
components. We note that uncertainty in the optical to radio celestial frame-
link also contributes to the position differences; there is approximately a 1
mas difference in the orientation of radio and optical frame when measured at
brighter magnitudes between the $Gaia$ DR2 and VLBI data of 26 radio stars
(Lindegren, 2020).

The $Gaia$ EDR3 “astrometric excess noise” value, which measure the
disagreement between the observations and the best-fitting astrometric model,
is 0.84 mas for RR Aql. As pointed out by Lindegren et al. (2021), this noise
not only includes the effect of the photo-center variability, but can also
come from modelling errors (e.g., from excess attitude noise of $Gaia$
satellite). However, assuming modelling errors are small, it is likely that
photo-center variability would be near the astrometric excess noise value,
which is $\approx$ 15% of its angular diameter. This is possible, since as
shown in the optical interferometric image of the Mira variable $\chi$ Cyg by
Lacour et al. (2009), a single hot spot can offset to the stellar center by up
to 30% of its diameter.

Table 5: Differences in positions of RR Aql from $Gaia$ EDR3 and VLBA

Epoch | R.A. Difference | Decl. Difference | 2D Difference
---|---|---|---
(yyyy/mm/dd) | (mas) | (mas) | (mas)
2017/10/29 | –0.22 $\pm$ 0.41 | 0.59 $\pm$ 1.19 | 0.63$\pm$1.26
2017/12/31 | 0.58 $\pm$ 0.42 | 1.00 $\pm$ 1.20 | 1.16$\pm$1.27
2018/03/29 | –0.49 $\pm$ 0.44 | 0.64 $\pm$ 1.20 | 0.81$\pm$1.28

Note. — The position differences are calculated by subtracting the values of
VLBA measurements from those of $Gaia$ EDR3, based on fitting an ring to the
SiO maser spot distributions. Column 1, 2 and 3 list the Right Ascension,
Declination differences and the 2-dimensional position differences.

##  5 Mira Period-Luminosity relation

Figure 9: Period analysis result (left panel) and phase-folded light curve
(right panel) of the visual data from 1904 November to 2021 August for RR Aql.
A pulsating period of 396 days is found.

As pointed out by Wood (2000), Miras are radial fundamental-mode pulsators,
while semi-regular variables (SRVs) can pulsate in the fundamental mode and/or
in overtones. Pulsation periods of SRVs can be hard to determine accurately,
since as their name indicates they often are not regular. Hence, we focus only
on Miras in this paper. For RR Aql, in order to determine its location on the
PLR diagram, we estimated a primary period of 396 $\pm$ 5 days (see the left
panel of Figure 9) using the visual data from 1904 November to 2021 August and
the data analysis tool VStar provided by the American Association of Variable
Star Observers (AAVSO)333https://www.aavso.org, using the Date Compensated
Discrete Fourier Transform (DCDFT) algorithm (Ferraz-Mello, 1981), which
corresponds to a maximum likelihood sinusoidal-plus-constant regression curve-
fitting. Compared to a DFT, the DCDFT provides a better estimation of the
power spectrum of a time series with uneven spacing. The right panel of Figure
9 shows a phase-folded visual light curve of RR Aql using the period of 396
days.

Figure 10:  The light curve at optical $V$ (lower panel) and a model near-
infrared $K$ band (upper panel) for RR Aql (here only the data from 1996 to
2005 are presented for clarity). The single-epoch 2MASS datum is shown in red.
The solid line shows the modeled optical light curve estimated by VStar. The
dashed line shows the near-infrared light curve converted from the modeled
optical light curve. The dash-dotted line indicates the mean magnitude
estimated by the light curve correction method for RR Aql.

The apparent magnitude at an effective wavelength of $2.2$ $\mu$m for RR Aql
($m_{K}$= 0.84 $\pm$ 0.26 mag) comes from the Two Micron All Sky Survey
(2MASS) point-source catalog. Since, Miras vary by upwards of 1 mag at 2$\mu$m
and 2MASS measurements are only at one epoch, this leads to extra uncertainty
in the mean magnitude owing to periodic variations. In order to better
estimate a mean apparent magnitude at near-infrared $K$ band, we adopted a
method similar to those demonstrated in Yuan et al. (2018) and Iwanek et al.
(2021b). This involved converting a modeled optical light curve to a near-
infrared one using an amplitude ratio ($A_{K}/A_{V}$=0.21 $\pm$ 0.09) and a
phase-lag ($\Delta\Phi_{KV}$=0.10 $\pm$ 0.07) between the two (Iwanek et al.,
2021b); then we added an infrared phase-adjusted constant (–0.44 mag) to the
measured magnitude. The mean magnitude for RR Aql using the light curve
correction method is 0.40 $\pm$ 0.33 mag. Figure 10 shows an example of the
light curve correction to mean magnitude for RR Aql. We estimate the
uncertainty of $m^{\rm mean}_{K}$ by taking into account the uncertainties
from both $m_{K}^{\rm 2MASS}$ and the amplitude ratio and phase-lag for light
curve correction.

Table 6 lists all Miras with both VLBI and $Gaia$ parallaxes so that we could
fit for PLR parameters separately for comparison. Where available, we adopted
mean magnitudes $m_{K}^{\rm Kagoshima}$ of time-monitored $m_{K}$ measurements
from the Kagoshima 1-m telescope, which used between 18 and 38 epochs per
star. Otherwise we used $m_{K}^{\rm mean}$ estimated with one-epoch 2MASS data
by the light curve correction method as done for RR Aql. Absolute magnitudes
of Miras at near-infrared $K$-band ($M_{K}$) were estimated by adopting the
parallax results from $Gaia$ and VLBI maser observations, separately. The
uncertainties in the absolute magnitudes, $M_{K}$ , takes into account
uncertainties of both $m_{K}$ and the parallax measurements used to convert
from apparent to absolute magnitude. All values used for Miras with both VLBI
and $Gaia$ parallaxes to fit Galactic PLRs are listed in Table 6.

It is clear from Figure 11 that, while a reasonable PLR is evident using
either VLBI or $Gaia$ parallaxes, there are some clear outliers. Rather than
remove outliers “by eye”, which has the potential for significant bias, in
order to fit a PLR we use a Bayesian approach with MCMC sampling which
intrinsically handles outliers without bias. Specifically we adopt the
“conservative formulation” of Sivia & Skilling (2006), for which the
probability distribution function for the uncertainty of a datum, $\sigma$, is
given by

${\rm prob}(\sigma|\sigma_{0},I)=\sigma_{0}/\sigma^{2}~{}~{}~{},$

where $\sigma_{0}$ is a (Gaussian) uncertainty of a “good” datum such that
$\sigma\geq\sigma_{0}$. Marginalizing over the unknown $\sigma$, leads to a
likelihood function proportional to $\ln\bigl{[}{1-e^{-R^{2}/2}\over
R^{2}}\bigr{]}$, where $R=(D-M)/\sigma_{0}$ and $D$ and $M$ are the datum and
model values. Unlike a Gaussian likelihood, this function has Lorentzian-like
tails, which assigns a significant probability to multi-$\sigma_{0}$ outliers.
As prior information for our fits, we adopt a slope of –3.60 $\pm$ 0.30 and a
zero point defined at logP(days) = 2.30 of –6.90 $\pm$ 1.00 mag. Our fitted
Mira PLR parameters are listed in Table 7. We note that the VLBI parallaxes
yield nearly a factor of two smaller uncertainty than the $Gaia$ parallaxes.

Table 6: Miras with both VLBI and $Gaia$ parallaxes

Star | $\Pi_{\rm VLBI}$ | $\Pi_{Gaia}$ |  $\Pi_{\rm VLBI}$-$\Pi_{Gaia}$ | Period | $m_{K}^{\rm 2MASS}$ | $m_{K}^{\rm Kagoshima}$ | $m_{K}^{\rm mean}$ | $\sigma^{\Pi}_{M_{K}}$ | $M^{\rm VLBI}_{K}$ | $M^{\rm Gaia}_{K}$ | Maser | Array | Ref.
---|---|---|---|---|---|---|---|---|---|---|---|---|---
| (mas) | (mas) | (mas) | (day) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) |  |  |
Y Lib | 0.86$\pm$0.05 | 0.83$\pm$0.08 | 0.03 $\pm$ 0.10 | 277 | 3.08$\pm$0.29 | 3.25$\pm$0.16 |  | 0.13 | –7.08 $\pm$ 0.20 | –7.15 $\pm$ 0.27 | H2O | VERA | k
X Hya | 2.07$\pm$0.05 | 2.53$\pm$0.11 | –0.46 $\pm$ 0.12 | 300 | 1.13$\pm$0.28 |  |  1.03 $\pm$ 0.28 | 0.05 | –7.39 $\pm$ 0.28 | –6.95 $\pm$ 0.30 | H2O | VERA | e
R UMa | 1.97$\pm$0.05 | 1.75$\pm$0.09 | 0.22 $\pm$ 0.10 | 302 | 1.11$\pm$0.20 | 1.19$\pm$0.02 |  | 0.06 | –7.34 $\pm$ 0.06 | –7.60 $\pm$ 0.11 | H2O | VERA | g
FV Boo | 0.97$\pm$0.06 | 1.01$\pm$0.09 | –0.04 $\pm$ 0.11 | 313 | 3.84$\pm$0.27 | 2.91$\pm$0.09 |  | 0.13 | –7.15 $\pm$ 0.16 | –7.06 $\pm$ 0.21 | H2O | VERA | i
R Cnc | 3.84$\pm$0.29 | 3.94$\pm$0.18 | –0.10 $\pm$ 0.34 | 357 | –0.71$\pm$0.17 |  | –0.78 $\pm$ 0.17 | 0.16 | –7.86 $\pm$ 0.24 | –7.80 $\pm$ 0.20 | H2O | VERA | e
S CrB | 2.39$\pm$0.17 | 2.60$\pm$0.11 | –0.21 $\pm$ 0.20 | 360 | –0.17$\pm$0.17 |  | –0.11 $\pm$ 0.18 | 0.15 | –8.22 $\pm$ 0.24 | –8.04 $\pm$ 0.20 | OH | VLBA | c
T Lep | 3.06$\pm$0.04 | 3.09$\pm$0.10 | –0.03 $\pm$ 0.11 | 368 | –0.27$\pm$0.35 |  | –0.57 $\pm$ 0.38 | 0.03 | –8.14 $\pm$ 0.38 | –8.12 $\pm$ 0.39 | H2O | VERA | e
S Ser | 1.25$\pm$0.04 | 0.77$\pm$0.13 | 0.48 $\pm$ 0.14 | 372 | 1.68$\pm$0.18 |  |  1.35 $\pm$ 0.24 | 0.07 | –8.17 $\pm$ 0.25 | –9.22 $\pm$ 0.44 | H2O | VERA | e
R Peg | 2.76$\pm$0.28 | 2.63$\pm$0.12 | 0.13 $\pm$ 0.30 | 378 | 0.38$\pm$0.36 |  |  0.68 $\pm$ 0.39 | 0.22 | –7.12 $\pm$ 0.45 | –7.22 $\pm$ 0.40 | H2O | VERA | e
R Hya | 7.93$\pm$0.18 | 6.74$\pm$0.46 | 1.19 $\pm$ 0.50 | 380 | –2.66$\pm$0.19 |  | –2.44 $\pm$ 0.22 | 0.05 | –7.94 $\pm$ 0.23 | –8.30 $\pm$ 0.27 | H2O | VERA | e
R Aqr | 4.59$\pm$0.24 | 2.59$\pm$0.33 | 2.00 $\pm$ 0.41 | 390 | –1.60$\pm$0.33 |  | –1.85 $\pm$ 0.35 | 0.11 | –8.54 $\pm$ 0.37 | –9.78 $\pm$ 0.45 | SiO | VERA | e
W Leo | 1.03$\pm$0.02 | 0.88$\pm$0.11 | 0.15 $\pm$ 0.11 | 392 | 1.80$\pm$0.23 |  |  1.47 $\pm$ 0.27 | 0.04 | –8.47 $\pm$ 0.27 | –8.81 $\pm$ 0.38 | H2O | VERA | e
RR Aql | 2.45$\pm$0.08 | 1.95$\pm$0.11 | 0.50 $\pm$ 0.14 | 396 | 0.84$\pm$0.26 |  |  0.40 $\pm$ 0.33 | 0.07 | –7.65 $\pm$ 0.34 | –8.15 $\pm$ 0.35 | H2O | VLBA | d
U Her | 3.76$\pm$0.27 | 2.36$\pm$0.08 | 1.40 $\pm$ 0.28 | 406 | –0.63$\pm$0.16 |  | –0.55 $\pm$ 0.17 | 0.16 | –7.67 $\pm$ 0.23 | –8.69 $\pm$ 0.18 | OH | VLBA | c
SY Scl | 0.75$\pm$0.03 | 0.52$\pm$0.12 | 0.23 $\pm$ 0.13 | 411 | 2.65$\pm$0.24 |  |  2.90 $\pm$ 0.27 | 0.09 | –7.72 $\pm$ 0.28 | –8.50 $\pm$ 0.57 | H2O | VERA | e
R Cas | 5.67$\pm$1.95 | 5.74$\pm$0.20 | –0.07 $\pm$ 1.96 | 430 | –1.40$\pm$1.00 |  | –1.90 $\pm$ 1.03 | 0.75 | –8.13 $\pm$ 1.27 | –8.10 $\pm$ 1.03 | OH | VLBA | a
AP Lyn | 2.00$\pm$0.04 | 2.02$\pm$0.12 | –0.02 $\pm$ 0.13 | 433 | 0.98$\pm$0.20 | 0.60$\pm$0.01 |  | 0.04 | –7.90 $\pm$ 0.05 | –7.87 $\pm$ 0.13 | H2O | VERA | h
U Lyn | 1.27$\pm$0.06 | 1.01$\pm$0.08 | 0.26 $\pm$ 0.10 | 434 | 1.53$\pm$0.23 | 1.15$\pm$0.09 |  | 0.10 | –8.33 $\pm$ 0.14 | –8.82 $\pm$ 0.20 | H2O | VERA | f
BX Cam | 1.73$\pm$0.03 | 1.76$\pm$0.10 | –0.03 $\pm$ 0.11 | 486 | 0.91$\pm$0.19 | 1.21$\pm$0.11 |  | 0.04 | –7.60 $\pm$ 0.12 | –7.55 $\pm$ 0.17 | H2O | VERA | h
V837 Her | 1.09$\pm$0.02 | 0.18$\pm$0.10 | 0.91 $\pm$ 0.10 | 520 | 2.06$\pm$0.27 | 1.68$\pm$0.03 |  | 0.04 | –8.13 $\pm$ 0.05 | –12.08$\pm$ 1.25 | H2O | VERA | h
UX Cyg | 0.54$\pm$0.06 | 0.70$\pm$0.09 | –0.16 $\pm$ 0.11 | 565 | 1.40$\pm$0.19 |  |  1.79 $\pm$ 0.26 | 0.24 | –9.55 $\pm$ 0.35 | –8.98 $\pm$ 0.39 | H2O | VLBA | b
OZ Gem | 0.81$\pm$0.04 | 0.45$\pm$0.33 | 0.35 $\pm$ 0.33 | 598 | 3.00$\pm$0.35 | 2.65$\pm$0.16 |  | 0.11 | –7.81 $\pm$ 0.19 | –9.05 $\pm$ 1.55 | H2O | VERA | j

Note. — Column 1 lists the names of Miras, columns 2 and 3 list VLBI and
$Gaia$ parallaxes, respectively, and column 4 lists VLBI and $Gaia$ parallax
differences. Column 5 lists the periods of magnitude variation at optical
wavelength from AAVSO, except for V837 Her and AP Lyn, for which periods are
from Chibueze et al. (2020). Apparent magnitudes at infrared $K$-band are
listed in column 6, 7 and 8: $m_{K}^{\rm 2MASS}$ values are from the 2MASS
point-source catalog (Skrutskie et al., 2006), $m_{K}^{\rm Kagoshima}$ are
from the time-monitored observations with the Kagoshima University 1-m
telescope and $m_{K}^{\rm mean}$ are the mean apparent magnitudes estimated
using the light curve correction method. Adopted to calculate the absolute
magnitude $M_{K}$ (listed in column 10 and 11) is $m_{K}^{\rm Kagoshima}$ if
available, otherwise $m_{K}^{\rm mean}$. The uncertainties of $M_{K}$ owing to
VLBI parallax uncertainties are listed in column 9. The maser species are
listed in column 12. The last two columns list the VLBI array for parallax
measurements and references, respectively. The $Gaia$ parallaxes are from Gaia
Collaboration et al. (2021). The VLBI parallaxes are from (a) Vlemmings et al.
(2003), (b) Kurayama et al. (2005), (c) Vlemmings & van Langevelde (2007), (d)
this paper, (e) VERA Collaboration et al. (2020), (f) Kamezaki et al. (2016b),
(g) Nakagawa et al. (2016), (h) Chibueze et al. (2020), (i) Kamezaki et al.
(2016a), (j) Urago et al. (2020), (k) Chibueze et al. (2019), (l) Kamezaki et
al. (2014).

Figure 11:  The PL($M_{K}$) relation for Galactic Miras (dashed line)
established in this paper using an “outlier tolerant” Bayesian MCMC approach,
fitted to data from VLBI (left panel) and $Gaia$(right panel) parallaxes. The
location of RR Aql on the PL diagram is given by a filled circle with its
name.

As shown by Iwanek et al. (2021a) and Iwanek et al. (2021b), a sample of
O-rich Miras with a large period range clearly shows a non-linear PLR, while a
sample restricted to log $P$ $\approx$ 2.0 - 2.8 can be well described by a
linear PLR. To compare our Galactic PLR with an even smaller range of log $P$
from $\approx$ 2.4 - 2.8 to the existing PLRs in literature, we list only the
latest determined linear PLR of O-rich Miras measured at near-infrared
$K$-band for different galaxies in Table 7. Yuan et al. (2017) determined a
linear PLR based on about 160 O-rich Miras with period $\textless$ 400 days at
near-infrared wavelengths, and we adjusted the zero-point using the updated
distance modulus of 18.477 $\pm$ 0.026 mag for LMC from Pietrzyński et al.
(2019). We took account of uncertainties from both the photometric zero-point
($\sim$ 0.02 mag) and distance modulus ($\sim$ 0.03 mag), which were not
included in Yuan et al. (2017). The linear PLR for M 33 (Yuan et al., 2018) is
from near-infrared light curves of 1169 O-rich Miras with periods $\textless$
400 days assuming a slope of –3.77, which is the same as the LMC value from
Yuan et al. (2017). In addition, assuming a distance modulus of 18.477 $\pm$
0.026 mag from Pietrzyński et al. (2019), Iwanek et al. (2021b) determined the
PLR of Miras in LMC at near-infrared $K$-band using a sample of 29 O-rich
Miras.

Since we rely on a strong prior for the slope of the PLR, we now only compare
zero-points for the LMC and M 33 and with Milky Way values based on VLBI and
Gaia calibrations. It is clear that the PLR zero-points of LMC and M 33
derived by Yuan et al. (2017, 2018) differ by only $0.05$ mag and are
consistent with each other. The zero-points for the Milky Way based on VLBI
and Gaia parallaxes are statistically consistent with each other, with the
VLBI value being more accurate. Compared with our Galactic PLR using VLBI
parallaxes, the LMC zero-point derived by Yuan et al. (2017) differs by 0.13
mag. The PLR zero-point derived by Iwanek et al. (2021b) has a larger
uncertainty and differs by 0.12 mag. All these zero-point differences are
consistent within their joint uncertainties.

We derive “parallaxes” (from PLR-based distances) of 2.05 $\pm$ 0.31 mas and
2.45 $\pm$ 0.41 mas for RR Aql, using the LMC $K$-band PLRs of Yuan et al.
(2017) and Iwanek et al. (2021b), respectively. While the latter PLR-based
parallax is close to our VLBI parallax of 2.44 $\pm$ 0.07 mas, both are
consistent within their uncertainties.

The scatter of PLR post-fit residuals listed in Table 7 for Galactic Miras is
larger than for the LMC or M 33. Two likely reasons include 1) the Miras in
other galaxies are at nearly the same distance, which removes distance
uncertainty as a source of scatter, and 2) more than half of the Galactic
Miras listed in Table 6 have only a single epoch of 2MASS magnitudes and these
were measured with saturated detectors. The 2MASS apparent magnitude for these
Miras at near-infrared $K$-band are brighter than $4^{th}$ mag, which is the
saturation limit for the 51 ms exposures, leading to larger magnitude
uncertainty. Fluxes for these Miras were estimated by template fitting to the
unsaturated scattered light in the wings of the saturated star image
(Skrutskie et al., 2006). An examination of Table 6, suggests a typical
uncertainty of the mean apparent 2MASS magnitudes is $\approx 0.3$ mag.
Removing a scatter of this magnitude from the total scatter of 0.45 mag for
the VLBI calibrated PLR, we find a residual scatter of 0.34 mag (ie,
$\sqrt{0.45^{2}-0.3^{2}}$). This residual scatter is currently unaccounted
for, but could come from slightly underestimated parallax uncertainties and/or
intrinsic brightness variations over different stellar cycles.

Table 7: Coefficients of PLR of Miras at near-infrared $K$-band in different
galaxies

Galaxies | Slope ($a$) | Zero-point ($b$) | Scatter | $N$ | References
---|---|---|---|---|---
|  | (mag) | (mag) |  |
LMC | –3.77 $\pm$ 0.07 | –6.92 $\pm$ 0.04 | 0.12 | 158 | Yuan et al. (2017)
M 33 | [–3.77] | –6.97 $\pm$ 0.01 | 0.21 | 1169 | Yuan et al. (2018)
LMC | –3.30 $\pm$ 0.46 | –6.67 $\pm$ 0.06 | 0.27 | 29 | Iwanek et al. (2021b)
Milky Way | –3.59 $\pm$ 0.29 | –6.79 $\pm$ 0.15 | 0.45 | 22 | This paper using VLBI parallaxes
Milky Way | –3.63 $\pm$ 0.30 | –7.08 $\pm$ 0.29 | 0.94 | 22 | This paper using $Gaia$ parallaxes

Note. — PLRs are defined as $M=a({\rm logP(days)}-2.30)+b$ for Miras at near-
infrared $K$-band, matching those of Yuan et al. (2017, 2018) and Iwanek et
al. (2021b). Columns 2 and 3 give the PLR slopes and zero-points. The Milky
Way fits used a strong prior of –3.60 $\pm$ 0.30 for the slope and a very weak
prior of –6.90 $\pm$ 1.00 mag for the zero-points. Columns 4 and 5 list the
standard deviation of the post-fit residuals in magnitudes and the number of
Miras fitted. The results for the LMC by Yuan et al. (2017) were based on
Miras with periods $\textless~{}400$ days, assuming distance modulus of 18.493
$\pm$ 0.048 mag (Pietrzyński et al., 2013); we have adjusted the zero-point
and its uncertainty to an updated distance modulus of 18.477 $\pm$ 0.026 mag
(Pietrzyński et al., 2019). The slope of the PLR for M 33 by Yuan et al.
(2018) was fixed at –3.77 to match the LMC slope determined by Yuan et al.
(2017). The PLR for the LMC by Iwanek et al. (2021b) assumed a distance
modulus of 18.477 $\pm$ 0.026 mag (Pietrzyński et al., 2019); the zero-point
uncertainty includes both the statistical and systematic errors.

##  6 Conclusions

We have measured the trigonometric parallax and proper motion of the Mira
variable RR Aql using VLBA observations of its water and SiO masers. Our
parallax is 2.44 $\pm$ 0.07 mas, which differs significantly from the $Gaia$
EDR3 parallax of 1.95 $\pm$ 0.11 mas. We fitted the positions of SiO masers to
a circle and obtained the absolute position of the central star. Comparing the
absolute positions of RR Aql from VLBA and $Gaia$, we find agreement within
their joint uncertainties.

Using 22 Miras, including RR Aql, which have both VLBI and $Gaia$ parallaxes,
we have derived separately PLRs of Miras at near-infrared K band. The post-fit
residuals using VLBI parallaxes are about half as large as those using $Gaia$
parallaxes, suggesting that VLBI parallaxes are more accurate than $Gaia$
parallaxes for Miras. Thus, VLBI astrometry should be effective to validate
future $Gaia$ results for Miras.

Compared with the PLRs for Miras in other galaxies, the Galactic PLR based on
the VLBI and $Gaia$ parallaxes are statistically consistent with PLRs in LMC
and M 33 at near-infrared $K$-band. Since Miras vary by about 1 mag at
infrared $K$-band, absolute magnitudes of Miras are uncertain by about
$\pm~{}0.4$ mag when photometry from only a single epoch is used. This often
dominates over parallax uncertainty when estimating absolute magnitudes. Thus,
priority should be given to monitoring surveys of Miras at infrared $K$-band
in order to obtain accurate measurements of mean $m_{K}$. Lacking multi-epoch
magnitudes, combining optical light curves with infrared magnitudes can help
to reduce the uncertainty due to the periodic variation.

In the era of space infrared telescopes, e.g., JWST, one could obtain
extremely accurate apparent magnitudes for Miras. In addition, future
astrometric results of $Gaia$, which could be validated by VLBI observations,
are expected to greatly increase the number of Miras with accurate parallaxes.
More reliable Galactic PLRs of Miras will be available, leading to Miras being
widely used as distance indicators.

This work was partly supported by the National Science Foundation of China
under grant U1831136 and U2031212, and the Key Laboratory for Radio Astronomy,
Chinese Academy of Sciences.

## References

  * Benson et al. (1990) Benson, P. J., Little-Marenin, I. R., Woods, T. C., et al. 1990, ApJS, 74, 911, doi: 10.1086/191526
  * Chiavassa et al. (2011) Chiavassa, A., Pasquato, E., Jorissen, A., et al. 2011, A&A, 528, A120, doi: 10.1051/0004-6361/201015768
  * Chibueze et al. (2019) Chibueze, J. O., Omodaka, T., Urago, R., et al. 2019, PASJ, 71, 92, doi: 10.1093/pasj/psz075
  * Chibueze et al. (2020) Chibueze, J. O., Urago, R., Omodaka, T., et al. 2020, PASJ, 72, 59, doi: 10.1093/pasj/psaa056
  * Deller et al. (2007) Deller, A. T., Tingay, S. J., Bailes, M., & West, C. 2007, PASP, 119, 318, doi: 10.1086/513572
  * Eyer & Cuypers (2000) Eyer, L., & Cuypers, J. 2000, in Astronomical Society of the Pacific Conference Series, Vol. 203, IAU Colloq. 176: The Impact of Large-Scale Surveys on Pulsating Star Research, ed. L. Szabados & D. Kurtz, 71–72
  * Feast & Whitelock (2014) Feast, M., & Whitelock, P. A. 2014, in IAU Symposium, Vol. 298, IAU Symposium, ed. S. Feltzing, G. Zhao, N. A. Walton, & P. Whitelock, 40–52, doi: 10.1017/S1743921313006182
  * Ferraz-Mello (1981) Ferraz-Mello, S. 1981, AJ, 86, 619, doi: 10.1086/112924
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Glass & Lloyd Evans (1981) Glass, I. S., & Lloyd Evans, L. 1981, Nature, 291, 303, doi: 10.1038/291303a0
  * Huang et al. (2018) Huang, C. D., Riess, A. G., Hoffmann, S. L., et al. 2018, ApJ, 857, 67, doi: 10.3847/1538-4357/aab6b3
  * Huang et al. (2020) Huang, C. D., Riess, A. G., Yuan, W., et al. 2020, ApJ, 889, 5, doi: 10.3847/1538-4357/ab5dbd
  * Iwanek et al. (2021a) Iwanek, P., Soszyński, I., & Kozłowski, S. 2021a, ApJ, 919, 99, doi: 10.3847/1538-4357/ac10c5
  * Iwanek et al. (2021b) Iwanek, P., Kozłowski, S., Gromadzki, M., et al. 2021b, ApJS, 257, 23, doi: 10.3847/1538-4365/ac1797
  * Kamezaki et al. (2014) Kamezaki, T., Kurayama, T., Nakagawa, A., et al. 2014, PASJ, 66, 107, doi: 10.1093/pasj/psu112
  * Kamezaki et al. (2016a) Kamezaki, T., Nakagawa, A., Omodaka, T., et al. 2016a, PASJ, 68, 75, doi: 10.1093/pasj/psw068
  * Kamezaki et al. (2016b) —. 2016b, PASJ, 68, 71, doi: 10.1093/pasj/psv030
  * Karovicova et al. (2011) Karovicova, I., Wittkowski, M., Boboltz, D. A., et al. 2011, A&A, 532, A134, doi: 10.1051/0004-6361/201116527
  * Kemball & Diamond (1997) Kemball, A. J., & Diamond, P. J. 1997, ApJ, 481, L111, doi: 10.1086/310664
  * Kurayama et al. (2005) Kurayama, T., Sasao, T., & Kobayashi, H. 2005, ApJ, 627, L49, doi: 10.1086/432051
  * Lacour et al. (2009) Lacour, S., Thiébaut, E., Perrin, G., et al. 2009, ApJ, 707, 632, doi: 10.1088/0004-637X/707/1/632
  * Lindegren (2020) Lindegren, L. 2020, A&A, 637, C5, doi: 10.1051/0004-6361/201936161e
  * Lindegren et al. (2021) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2, doi: 10.1051/0004-6361/202039709
  * Nakagawa (2009) Nakagawa, A. 2009, in AGB Stars and Related Phenomena, ed. T. Ueta, N. Matsunaga, & Y. Ita, 58
  * Nakagawa et al. (2016) Nakagawa, A., Kurayama, T., Matsui, M., et al. 2016, PASJ, 68, 78, doi: 10.1093/pasj/psw069
  * Pietrzyński et al. (2013) Pietrzyński, G., Graczyk, D., Gieren, W., et al. 2013, Nature, 495, 76, doi: 10.1038/nature11878
  * Pietrzyński et al. (2019) Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Nature, 567, 200, doi: 10.1038/s41586-019-0999-4
  * Reid et al. (2009) Reid, M. J., Braatz, J. A., Condon, J. J., et al. 2009, ApJ, 695, 287, doi: 10.1088/0004-637X/695/1/287
  * Reid & Honma (2014) Reid, M. J., & Honma, M. 2014, ARA&A, 52, 339, doi: 10.1146/annurev-astro-081913-040006
  * Reid & Menten (1991) Reid, M. J., & Menten, K. 1991, in Astronomical Society of the Pacific Conference Series, Vol. 16, Atoms, Ions and Molecules: New Results in Spectral Line Astrophysics, ed. A. D. Haschick & P. T. P. Ho, 375–379
  * Reid & Menten (1997) Reid, M. J., & Menten, K. M. 1997, ApJ, 476, 327, doi: 10.1086/303614
  * Richards et al. (2010) Richards, A. M. S., Assaf, K., Bains, I., et al. 2010, in 10th European VLBI Network Symposium and EVN Users Meeting: VLBI and the New Generation of Radio Arrays
  * Sivia & Skilling (2006) Sivia, D. S., & Skilling, J. 2006, Data Analysis–A Bayesian Tutorial, 2nd edn. (Oxford Science Publications)
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Urago et al. (2020) Urago, R., Yamaguchi, R., Omodaka, T., et al. 2020, PASJ, 72, 57, doi: 10.1093/pasj/psaa024
  * van Leeuwen (2007) van Leeuwen, F. 2007, A&A, 474, 653, doi: 10.1051/0004-6361:20078357
  * van Moorsel et al. (1996) van Moorsel, G., Kemball, A., & Greisen, E. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 37
  * VERA Collaboration et al. (2020) VERA Collaboration, Hirota, T., Nagayama, T., et al. 2020, PASJ, 72, 50, doi: 10.1093/pasj/psaa018
  * Vlemmings & van Langevelde (2007) Vlemmings, W. H. T., & van Langevelde, H. J. 2007, A&A, 472, 547, doi: 10.1051/0004-6361:20077897
  * Vlemmings et al. (2003) Vlemmings, W. H. T., van Langevelde, H. J., Diamond, P. J., Habing, H. J., & Schilizzi, R. T. 2003, A&A, 407, 213, doi: 10.1051/0004-6361:20030766
  * Whitelock et al. (2008) Whitelock, P. A., Feast, M. W., & van Leeuwen, F. 2008, MNRAS, 386, 313, doi: 10.1111/j.1365-2966.2008.13032.x
  * Willson (1981) Willson, L. A. 1981, in Astrophysics and Space Science Library, Vol. 89, IAU Colloq. 59: Effects of Mass Loss on Stellar Evolution, ed. C. Chiosi & R. Stalio, 353–355, doi: 10.1007/978-94-009-8500-1_52
  * Wood (2000) Wood, P. R. 2000, PASA, 17, 18, doi: 10.1071/AS00018
  * Xu et al. (2019) Xu, S., Zhang, B., Reid, M. J., Zheng, X., & Wang, G. 2019, ApJ, 875, 114, doi: 10.3847/1538-4357/ab0e83
  * Yoon et al. (2018) Yoon, D.-H., Cho, S.-H., Yun, Y., et al. 2018, Nature Communications, 9, 2534, doi: 10.1038/s41467-018-04767-8
  * Yuan et al. (2017) Yuan, W., Macri, L. M., He, S., et al. 2017, AJ, 154, 149, doi: 10.3847/1538-3881/aa86f1
  * Yuan et al. (2018) Yuan, W., Macri, L. M., Javadi, A., Lin, Z., & Huang, J. Z. 2018, AJ, 156, 112, doi: 10.3847/1538-3881/aad330


# Quantum Hoare Type Theory: Extended Abstract

Kartik Singhal John Reppy University of Chicago {ks<EMAIL_ADDRESS>
###### Abstract

As quantum computers become real, it is high time we come up with effective
techniques that help programmers write correct quantum programs. In classical
computing, formal verification and sound static type systems prevent several
classes of bugs from being introduced. There is a need for similar techniques
in the quantum regime. Inspired by Hoare Type Theory [NMB08] in the classical
paradigm, we propose Quantum Hoare Types by extending the Quantum IO Monad
[AG09] by indexing it with pre- and postconditions that serve as program
specifications. In this paper, we introduce Quantum Hoare Type Theory (QHTT),
present its syntax and typing rules, and demonstrate its effectiveness with
the help of examples.

QHTT has the potential to be a unified system for programming, specifying, and
reasoning about quantum programs. This is a work in progress.111An updated
report is available [Sin20] at the time of this publication.

##  1 Introduction

It is difficult to reason about the correctness of quantum programs. Sound
static type systems help prevent a huge class of bugs from occurring but since
the realm of quantum programming is still new there is not a lot of consensus
on what kind of types make the most sense. Further, it is unclear how much
they help programmers reason about the semantic properties associated with the
quantum algorithms that they are implementing.

Recent work [HM18, HM19] as part of EPiQC222EPiQC: Enabling Practical-Scale
Quantum Computation: https://epiqc.cs.uchicago.edu has identified several
classes of bugs in quantum programs and proposed approaches to tackle them.
The technique that holds the most promise is assertion checking using
preconditions and postconditions. But assertions are usually checked
dynamically during runtime, which can be wasteful of precious quantum
computing or simulation resources.

As programming languages researchers, we think it will be better to encode
such assertions into a static type system both for formal verification and to
aid the programmers in writing correct programs from the start. Inspired by
the use of Hoare triples in the verification of imperative programs and
building on the idea of Hoare Type Theory (HTT) for classical programming
languages [NMB08], we extend the Quantum IO Monad interface [AG09] and propose
Quantum Hoare Type Theory aimed at enabling both sound static type checking
and formal verification of quantum programs.

The main idea is that instead of just indexing the QIO monad type with the
type of computation result, we can also index it with preconditions and
postconditions so as to integrate Hoare-style reasoning into the type system
itself. The resulting Quantum Hoare Type, $\\{P\\}x{:}A\\{Q\\}$, specifies the
preconditions, $P$, that hold on the quantum state before execution; the
result, $x$, and its type, $A$; and, the postconditions, $Q$, that are true
for the quantum state after successful execution. In this way the effectful
quantum fragment of the program is effectively encapsulated inside the Quantum
Hoare monad. Our theory, like HTT, allows the usual equational reasoning for
the pure classical fragments of the program and uses syntax-directed type
checking for generating strongest postconditions for the quantum fragment.

This paper is organized as follows. We discuss some related work and
background in the next section, section 2; specifically, Hoare Type Theory in
the classical setting (section 2.1) and then the Quantum IO Monad interface
(section 2.2). Then we introduce our contributions in merging these ideas
together for specification and enforcement of useful semantic properties in
section section 3 along with its syntax (section 3.1), typing rules (section
3.2) and, some examples (section 3.3). Finally, we conclude and share ideas
for future work in the final section (section 4).

##  2 Related Work and Background

Previous work, such as Proto-Quipper [Ros15, RS17] and QWire [PRZ17, RPZ18,
Ran18], utilize a linear type system and dependent types to enforce a small
subset of semantic properties, such as the no-cloning restriction and whether
a unitary gate is of the right dimension. These advances in quantum type
systems, although helpful, still fall short in encoding and enforcing even
more useful properties that one would like to be able to express for the
purpose of verification.

Our approach builds upon previous work in reasoning about quantum programs
such as Quantum Weakest Preconditions [DP06] and Quantum Hoare Logic [Yin12]
in the spirit of Hoare [Hoa69] and Dijkstra [Dij76] but attempts to bring
those reasoning techniques into the type system. The hope is that programmers
will be able to encode some of the semantic properties that they expect of
their programs as specifications in their code and type checking will ensure
correctness of some of those properties. In the classical setting, Hoare Type
Theory [NMB08] accomplishes exactly this goal. Our attempt is to merge these
ideas for the quantum realm.

In the rest of this section, we provide background on the core ideas from
existing literature that form the foundation for our work — Hoare Type Theory
(HTT) and the Quantum IO Monad interface. We assume background in basics of
quantum computing and refer the reader to Nielsen’s excellent series of essays
[MN19] for a first introduction and to the standard textbook [NC10] for
advanced material.

###  2.1 Hoare Type Theory

A Hoare type, $\Delta.\Psi.\\{P\\}\;\texttt{x}:A\;\\{Q\\}$, encodes
preconditions and postconditions in the same spirit as Hoare triples to allow
both specification and verification of effectful classical programs. It can be
read as ‘for a stateful computation executed in a heap that satisfies
precondition $P$, return a value of type $A$ in a heap that satisfies
postcondition $Q$.’ The contexts, $\Delta$ and $\Psi$, contain the variables
and heap variables that may appear in both $P$ and $Q$. In this presentation,
we will omit the contexts in the Hoare type when they are unneeded. We show an
example to demonstrate the expressiveness of a Hoare type — the alloc
primitive from HTT [NMB08]: alloc : ∀α. Πx:α. {emp} y:nat {y ↦α x}

This type specifies that alloc is a polymorphic function that takes as input a
variable, $x$, of any simple type, $\alpha$, that is executed in an empty heap
(meaning it does not affect existing heap), returns a new location bound by a
fresh variable, $y$, of type nat and initializes it with the supplied value
$x$ of type $\alpha$.

Type checking in HTT involves generation of strongest postconditions at each
step of the program. Verification is a two-step process: the first phase does
basic type checking and verification condition generation which is decidable,
the second phase needs to show the validity of the generated verification
conditions, which can be undecidable. This second phase can be deferred to an
automated theorem prover.

We provide more details when we introduce QHTT in section section 3. An
accessible introduction to Hoare Type Theory is available in the lecture notes
by Perconti [Per12].

###  2.2 Quantum IO Monad

The Quantum IO (QIO) monad [AG09] is a purely functional interface for quantum
programming that provides a separation between the unitary (reversible) and
non-unitary (irreversible) fragments of quantum computation. It provides
isolation of quantum effects inside a monad similar to what we need to provide
in our Quantum Hoare monad.

QIO interface was developed as a library for the Haskell programming language
but its design was influenced [Gre10] by the category FQC of finite quantum
computations (where computations are interpreted as superoperators) as
explored by the authors in previous work [GA08]. This is the most we will say
about the relation of QIO monad to category theory in this paper.

Here we discuss the relevant bits from the QIO interface. QIO provides a type
for referring to qubits, Qbit; a type for unitary operations, U; and, the QIO
type operator which is a monad indexed by the result of the quantum
computation. The primitive quantum operations defined using the QIO monad
below (Haskell syntax) are:

⬇

mkQbit :: Bool $\rightarrow$ QIO Qbit \-- initialization

applyU :: U $\rightarrow$ QIO () \-- apply a unitary to quantum state

measQbit :: Qbit $\rightarrow$ QIO Bool \-- measurement

Further, U is defined as a monoid with neutral element, mempty, and operation,
mappend, that encodes reversible operations on quantum state. The core unitary
primitives that we will need are:

⬇

rot :: Qbit $\rightarrow$ ((Bool, Bool) $\rightarrow$ C) $\rightarrow$ U

cond :: Qbit $\rightarrow$ (Bool $\rightarrow$ U) $\rightarrow$ U

rot takes a qubit and a two-by-two complex-valued unitary matrix (represented
as a function from the matrix indices to $\mathbb{C}$) and lets one define
arbitrary rotation on a single qubit. rot can then be used to define standard
gates such as the Hadamard, H, and the Pauli gates, X, Y and Z. cond is the
unitary conditional that is used to perform branching unitary operations. For
the present paper, we use a simpler conditional, ifQ, that is defined as
follows:

⬇

ifQ :: Qbit $\rightarrow$ U $\rightarrow$ U

ifQ q u $=$ cond q (\x $\rightarrow$ if x then u else mempty)

That is, ifQ acts as the standard control operator in the quantum circuit
model that runs its second argument based on the truth value of its first
argument. Then, it is easy to define controlled gates such as CNOT with the
expression ifQ q1 (X q2) .

Quantum state in the QIO interface is modeled as a normalized vector that
stores pairs representing complex amplitudes associated with each basis state
(represented as a map from Qbit to Bool types). The operations defined over
the vector class ensure that the quantum state is kept normalized throughout.

QIO further ensures that the monoidal structure of unitary operations lets one
run a unitary operation over the complete quantum state. We elide details here
and refer the reader to [AG09] for more.

##  3 Quantum Hoare Type Theory

In this section, we introduce the Quantum Hoare Type inspired by the QIO monad
and its type theory (QHTT) by replacing the classical effectful portion of
Hoare Type Theory with quantum effects. The core idea is to encapsulate any
quantum effect inside a monadic Quantum Hoare type,
$\Delta.\Psi.\\{P\\}x{:}A\\{Q\\}$, and formalize reasoning of quantum effects
using strongest postconditions in a similar fashion as HTT in the classical
setting.

The next two sections discuss the syntax and typing of QHTT. The last section
show examples written in QHTT.

###  3.1 Syntax

Figure 1 shows the syntax of Quantum Hoare Type Theory. Our presentation
closely follows that of HTT [NMB08].

Types | $A,B,C$ | ::= |  $\textrm{{1}}\mid\textrm{{Bool}}\mid\textrm{{Qbit}}\mid\textrm{{U}}\mid\textrm{{Pure}}\mid A\otimes B\mid\Pi x{:}A.B\mid\Delta.\Psi.\\{P\\}x{:}A\\{Q\\}$
---|---|---|---
Assertions | $P,Q,R$ | ::= |  $\top\mid\bot\mid P\wedge Q\mid P\vee Q\mid P\supset Q\mid\neg P\mid\exists x{:}A.P\mid\forall x{:}A.P\mid\exists h{:}\,\mathrm{\texttt{heap}}\,.P\mid\forall h{:}\,\mathrm{\texttt{heap}}\,.P\mid\,\mathrm{\texttt{Id}}\,_{\\_}A(M,N)\mid\,\mathrm{\texttt{HId}}\,(H,G)\mid\,\mathrm{\texttt{indom}}\,(H,M)$
Q. Heaps | $H,G$ | ::= | $h\mid\,\mathrm{\texttt{empty}}\,\mid\,\mathrm{\texttt{upd}}\,(H,M,N)$
Elim terms | $K,L$ | ::= | $x\mid KM\mid M:A$
Intro terms | $M,N,O$ | ::= | $K\mid()\mid\lambda x.M\mid\,\mathrm{\texttt{do}}\,E\mid\,\mathrm{\texttt{true}}\,\mid\,\mathrm{\texttt{false}}\,$
Q. Commands | $c$ | ::= |  $\,\mathrm{\texttt{mkQbit}}\,M\mid\,\mathrm{\texttt{measQbit}}\,M\mid\,\mathrm{\texttt{applyU}}\,M\mid\,\mathrm{\texttt{if}}\,M\,\mathrm{\texttt{then}}\,N\,\mathrm{\texttt{else}}\,O$
Computations | $E,F$ | ::= | $\,\mathrm{\texttt{return}}\,M\mid x\leftarrow K;E\mid x\Leftarrow c;E\mid x=_{\\_}AM;E$
Variable context | $\Delta$ | ::= | $\cdot\mid x:A,\Delta$
Q. Heap context | $\Psi$ | ::= | $\cdot\mid h,\Psi$
Assertion context | $\Gamma$ | ::= | $\cdot\mid P,\Gamma$
Figure 1: Syntax of Quantum Hoare Type Theory

##### Types

We include primitive types for unit, booleans, qubits and unitary operations,
and a type Pure for representing pure quantum state vectors (complex valued
vectors in Hilbert space); these correspond to similar types in the QIO work.
There are also type constructors for pairs, dependent functions and the Hoare
type from HTT.

##### Assertions

Apart from the usual first order logic, we have assertions for reasoning about
propositional equality of terms, $\,\mathrm{\texttt{Id}}\,_{\\_}A(M,N)$, and
heaps, $\,\mathrm{\texttt{HId}}\,(H,G)$.

As done in HTT, some common convenience assertions can be defined using the
base primitives. Particularly, emp denotes that the current heap, say $h$, is
empty, that is, $\,\mathrm{\texttt{HId}}\,(h,\,\mathrm{\texttt{empty}}\,)$.
$M\mapsto N$ represents the only mapping in the (singleton) heap.
$M\hookrightarrow N$ says that looking up $M$ in the heap returns $N$ which
can be more explicitly written as $\,\mathrm{\texttt{seleq}}\,(H,M,N)$.

##### Quantum State

We use the terms quantum state and quantum heaps interchangeably throughout. A
quantum heap is a partial function from qubits (that can in turn be thought of
as locations, represented using natural numbers) to at most one quantum state
vector.

For the purpose of the type theory, heaps are functional, so that, for
example, $\,\mathrm{\texttt{upd}}\,(H,M,N)$ returns a new heap after updating
heap $H$ at location $M$ with $N$.

##### Terms

The terms are divided into introduction and elimination sorts for the purpose
of bidirectional typechecking like in HTT.

The pure fragment allows higher-order functions while the impure (effectful)
fragment is encapsulated within the monadic do $E$ construct and supports a
simple imperative quantum language inspired by QIO constructs. It is worth
noting that do $E$ represents suspended computation and is considered pure.

##### Quantum Commands

These are the quantum-specific effectful commands from QIO for initialization,
unitary application and measurement of qubits. We also support classical
control using the if $M$ then $N$ else $O$ construct similar to HTT.

##### Computations

The computations are the usual monadic return and three sequencing operations
that bind their result to a variable similar to HTT. The first, $x\leftarrow
K$, executes a suspended computation (such as applying a function that
includes a do $E$ in its body); the second, $x\Leftarrow c$ executes a
primitive command $c$; and the last, $x=_{\\_}AM$, is just syntactic sugar for
let-binding.

##### Contexts

Contexts are modelled exactly the same way as in HTT. The only difference is
that instead of classical heaps, we have quantum heaps.

###  3.2 Typing Rules

In this section, we provide some details about type checking in HTT and how it
is different in the quantum case.

####  3.2.1 Judgments

We explain the main judgments used in HTT that are required to understand the
typing rules in the next section. The type system is designed to be
bidirectional and syntax directed and hence includes separate judgments for
intro and elim terms. We are able to use them unchanged for QHTT. A recent
survey provides accessible introduction to bidirectional typing [DK19].

The type checking process also involves generating canonical forms for each
term. A canonical form in HTT means a beta-normal (containing no beta-redexes)
and eta-long (all intro terms are eta-expanded) form. We elide much detail
here.

As usual, elim terms synthesize types
($\,{\color[rgb]{0.85,0.0,0.0}{\Rightarrow}}\,$) and intro terms are checked
against the given type ($\,{\color[rgb]{0.0,0.0,0.85}\Leftarrow}\,$).
Canonical forms are always synthesized. For the pure fragment:

$\Delta\vdash K\,{\color[rgb]{0.85,0.0,0.0}{\Rightarrow}}\,A[N^{\prime}]$ | Elim term K has type A and canonical form $N^{\prime}$
---|---
$\Delta\vdash M\,{\color[rgb]{0.0,0.0,0.85}\Leftarrow}\,A[M^{\prime}]$ | Intro term M has type A and canonical form $M^{\prime}$

Note that the blue colored $\,{\color[rgb]{0.0,0.0,0.85}\Leftarrow}\,$ symbol
we use in the typing rules lives in a different syntactic category from the
$\Leftarrow$ symbol used in the syntax (fig. 1) that binds the output of
primitive commands.

The judgments for computations involve synthesizing the strongest
postconditions and checking whether a given postcondition applies:

$\Delta;P\vdash E\,{\color[rgb]{0.85,0.0,0.0}{\Rightarrow}}\,x:A.Q[E^{\prime}]$ | Computation E with precondition P has strongest postcondition Q
---|---
| and returns value x of type A. Its canonical form is $E^{\prime}$.
$\Delta;P\vdash E\,{\color[rgb]{0.0,0.0,0.85}\Leftarrow}\,x:A.Q[E^{\prime}]$ | Computation E with precondition P has postcondition Q
| and returns value x of type A. Its canonical form is $E^{\prime}$.

Even though, we have shown the full forms of these judgments above, we will
omit the canonical forms while presenting typing rules in the next section as
they crowd the rules and do not affect the insight to be gained from them.

Finally, the judgment $\Delta;\Psi;\Gamma_{\\_}1\Longrightarrow\Gamma_{\\_}2$
encodes the sequent calculus for the assertion logic. Recall that $\Delta$,
$\Psi$, and $\Gamma$ denote the variable context, the quantum heap context and
the assertion context, respectively.

The type system includes rules for primitive effectful commands and those for
structuring composition such as monadic unit and bind. We reuse most of the
HTT rules except those for primitives for quantum effects that need to specify
the strongest postconditions for each primitive, we discuss that next.

####  3.2.2 Strongest Postconditions

Here we show the strongest postconditions for the primitive quantum commands
of QHTT encoded in their typing rules. Given a quantum command for
initialization, unitary application or measurement, its strongest
postcondition is an assertion that most precisely captures the relationship
between the initial state and the modified state after the execution. These
rules along with the consq (consequent) rule of HTT (that ensures weakening of
the given strongest postcondition to an arbitrary postcondition) work together
to ensure composition of verification.

We need to use the relational composition connective from HTT ($P\circ Q$)
that captures how heap evolves with computation. It basically reads: $Q$ holds
of the current heap which is obtained after modification of a prior heap for
which P holds.

We also use HTT’s difference operator ($\multimap$) below that captures
changes to only the interesting fragment of the heap without modifying the
rest of the heap.

We can now state some typing rules:

##### Initialization

$x\Leftarrow\,\mathrm{\texttt{mkQbit}}\,M;E$

Δ⊢M ⇐ Bool
Δ, x : Qbit; P ∘(x ↦ state (M)) ⊢E ⇒ y:B . Q Δ; P ⊢x ⇐ mkQbit M; E ⇒ y:B.
(∃x:Qbit.Q) where x is fresh and state is a function that translates a
classical representation of quantum state (such as boolean here) to the
equivalent quantum state vector.

Initialization can only be performed for a term $M$ that can be type checked
as a Bool. Then, we look at the rest of the computation, $E$, in a context
that adds $x$ of type Qbit under a precondition that extends the previous
precondition, $P$, with the strongest postcondition for initialization, that
is, the newly bound variable $x$ points to the representation of a new qubit
state; this should synthesize the strongest postcondition (with respect to the
expanded precondition in the context), $Q$, for $E$.

To avoid dangling variables, we need to existentially quantify $x$ in the
postcondition, $Q$, of the conclusion.

##### Measurement

$x\Leftarrow\,\mathrm{\texttt{measQbit}}\,M;E$

Δ⊢M ⇐ Qbit
Δ; Ψ; P ⟹(M ↪-)
Δ, x:Bool; P ∘((M ↦-) ⊸ emp ) ⊢E ⇒ y:B . Q Δ; P ⊢x ⇐ measQbit M; E ⇒ y:B.
(∃x:Bool.Q)

The measurement rule can be read in a similar way as the previous rule. But in
measurement, we need to additionally prove that the verification condition
(written as a sequent), the precondition $P$ implies that the qubit $M$ is
allocated, holds in the heap context $\Psi$. The strongest postcondition here
is that the fragment of heap that only refers to location M becomes empty.

As may be apparent, initialization is analogous to alloc and measurement to
dealloc and lookup primitives of HTT. We are working out precise details for
unitary application, but at a high-level, they involve ensuring that the given
unitary term actually represents a unitary matrix and updating the quantum
state with the result of unitary application. This involves incorporating the
monoidal structure of unitaries into our theory (as is done in QIO). This step
also involves generating verification conditions that cannot be checked during
this first phase of typechecking. We are also considering more tractable
alternatives such as the sum-over-paths action semantics [Amy18] based on the
Feynman path integral formulation.333We ended up taking a different approach
in subsequent work [Sin20].

###  3.3 Sample programs

In this section, we show some simple examples from the QIO work [AG09]
translated in our language and annotated with their program specifications
using their Quantum Hoare types. An observation is that translating the
programs was a very simple process except for coming up with the right
specifications for the programs. We provide some commentary on assertions
specified for these programs.

At the end of this section, we show a sample verification using QHTT for Bell
pair generation.

#### Hello Quantum World

⬇

hqw : {emp} r : Bool {emp $/$\ Id(r, false)}

$=$ do q $\Leftarrow$ mkQbit false;

measQbit q

In this trivial program where we initialize a new qubit with false and then
immediately measure it, we assert that the result is equal to false. The
complete specification also implicitly says that the existing quantum state is
not affected in any way as both the pre- and postcondition include the emp
assertion.

#### Coin Toss

⬇

rnd : {emp} r : Bool {emp}

$=$ do q $\Leftarrow$ mkQbit false;

applyU (H q);

measQbit q

In the case of quantum coin toss, the result can be in either of its two
boolean values false or true, hence, we need not specify any special
postcondition as type checking the result will be sufficient for correctness.

#### Testing Bell Pair

⬇

testBell : {emp} (a, b) : (Bool, Bool) {emp $/$\ Id(a, b)} \-- a,b \in {true,
false}

$=$ do qa $\Leftarrow$ mkQbit false;

applyU (H qa);

qb $\Leftarrow$ mkQbit false;

applyU (ifQ qa (X qb));

(measQbit qa, measQbit qb)

Here, we are asserting that the returned booleans $a$ and $b$ hold the same
value, which is what we expect from the first Bell state prepared in this
program.

##### Verification

We would like to describe how QHTT lets us verify the correctness of the above
program based on the specification provided as the type of the function. For
illustration, we show the assertions annotated as comments as we step through
the testBell program below. We have elided specific details about unitary
application as they are still being worked out.

⬇

testBell : {emp} (a, b) : (Bool, Bool) {emp $/$\ Id(a, b)}

\-- P0: emp

$=$ do qa $\Leftarrow$ mkQbit false;

\-- P1: P0 \o (qa $\mapsto$ $|$0\$>$)

applyU (H qa);

\-- P2: P1 \o ((qa $\mapsto$ $|$0\$>$) $\multimap$ (qa $\mapsto$ $|$$+$\$>$))

qb $\Leftarrow$ mkQbit false;

\-- P3: P2 \o (qb $\mapsto$ $|$0\$>$)

applyU (ifQ qa (X qb));

\-- P4: P3 \o ((qa $\mapsto$ $|$$+$\$>$, qb $\mapsto$ $|$0\$>$) $\multimap$
(qa, qb) $\mapsto$ $|$\Phi$+$\$>$)

(measQbit qa, measQbit qb)

\-- P5: P4 \o ((qa $\mapsto$ -) $\multimap$ emp) \o ((qb $\mapsto$ -)
$\multimap$ emp)

In the assertion P2, we show that only the portion of the quantum state
referred to by the applyU command is affected. Similarly, in P4 we see that
only the relevant portion of the heap containing the control and target qubits
used in the CNOT operation are affected. In this case, we also needed to
combine the two qubit state as they are now entangled. The last assertion, P5,
merely encodes the strongest postcondition for measurement for each of the two
qubits.

Finally, the consq rule of HTT ensures that the assertion P5 leads to the
specified postcondition in the type of the program.

We show some two more examples in the appendix for a modular version of Bell
pair generation (Appendix A) that shows use of higher-order features of QHTT
and another for teleportation protocol (Appendix B). Both of these involve
specifying assertions over quantum state whose theory we still need to work
out.

##  4 Conclusions and Future Work

In this paper, we have described our ongoing work on Quantum Hoare Type
Theory, which modifies HTT for quantum computing using ideas from the Quantum
IO monad work. This approach has the potential to be a unified system for
programming, specification, and reasoning about quantum programs. This is
active work in progress and more details will be available in a forthcoming
technical report [Sin20].

There are several avenues of future work to explore:

##### Mechanization

There are multiple implementations of HTT in Coq such as Ynot [NMS+08]. We are
similarly working on mechanizing QHTT in Coq for higher assurance of the
usefulness and soundness of our theory. This will also enable us to extract
verified circuits in a lower level quantum language such as OpenQASM for
execution on real machines.

##### Quantum Assertion Logic

An obvious need is to come up with an assertion logic similar to Separation
Logic [Rey02] for quantum computing so as to be able to reason about only the
interesting portions of the quantum state while still ensuring correctness of
non-local effects such as entanglement. Various Quantum Hoare Logics that
exist [Yin12] currently do not support frame rules that provide Separation
Logic its power.444Since the time of this writing, we learned about Unruh’s
work [Unr19] that we use in subsequent work [Sin20].

##### Linearity

Peter Selinger and collaborators have recently proposed a linear dependent
typed version of Proto-Quipper (dubbed Proto-Quipper-D) [Sel20, FKS20]. It is
an interesting challenge to reconcile linearity in our theory based on their
proposal.

##### Circuits as Arrows

Further, Proto-Quipper treats quantum circuits as first class citizens of the
language. We would like to explore modifying our theory to treat Quantum Hoare
types as arrows instead of as monads as was suggested by Vizzotto et al
[VAS06]. It makes sense from the perspective of sequential composition as
arrows can have an arbitrary number of input/outputs as opposed to monads.

##### Behavioral Types

Another venue for exploration is to incorporate more precise types that can
distinguish between qubits in pure classical state vs. those in superposition
vs. those in entanglement [JP09] such as those inspired by the various quantum
resource theories or the Heisenberg representation of quantum mechanics
[RSSL19, RSSL20]. This may help us provide more specific postconditions that
quantum programmers expect to hold true of their programs.

##### Classical Effects

Finally, it will be an interesting challenge to reconcile both classical and
quantum effects together in a single grand unified theory for effects.

### Acknowledgments

We thank Robert Rand and the anonymous reviewers for their feedback on a
previous draft of this paper. This work is funded by EPiQC, an NSF Expedition
in Computing, under grant CCF-1730449.

## References

  * [1]
  * [AG09] T. Altenkirch & A.S. Green (2009): _The Quantum IO Monad_.  In: Semantic Techniques in Quantum Computation, Cambridge University Press, pp. 173–205, 10.1017/CBO9781139193313.006.  Available at http://www.cs.nott.ac.uk/~psztxa/g5xnsc/chapter.pdf.
  * [Amy18] M. Amy (2018): _Towards Large-scale Functional Verification of Universal Quantum Circuits_.  In: Proc. QPL ’18, pp. 1–21, 10.4204/EPTCS.287.1.
  * [Dij76] E.W. Dijkstra (1976): _A Discipline of Programming_.  Prentice-Hall, Englewood Cliffs, NJ.
  * [DK19] J. Dunfield & N. Krishnaswami (2019): _Bidirectional Typing_.  Available at https://arxiv.org/abs/1908.05839.  Submitted to ACM Computing Surveys.
  * [DP06] E. D’hondt & P. Panangaden (2006): _Quantum Weakest Preconditions_.  Math. Struct. Comput. Sci. 16(3), pp. 429–451, 10.1017/S0960129506005251.  Available at https://www.cs.mcgill.ca/~prakash/Pubs/weakest_mscs.pdf.
  * [FKS20] P. Fu, K. Kishida & P. Selinger (2020): _Linear Dependent Type Theory for Quantum Programming Languages: Extended Abstract_.  In: Proc. LICS ’20, p. 440–453, 10.1145/3373718.3394765.
  * [GA08] A.S. Green & T. Altenkirch (2008): _From Reversible to Irreversible Computations_.  Electron. Notes Theor. Comput. Sci. 210, pp. 65–74, 10.1016/j.entcs.2008.04.018.  Proc. QPL ’06.
  * [Gre10] A.S. Green (2010): _Towards a formally verified functional quantum programming language_.  Ph.D. thesis, University of Nottingham.  Available at http://eprints.nottingham.ac.uk/11457/.
  * [HM18] Y. Huang & M. Martonosi (2018): _QDB: From Quantum Algorithms Towards Correct Quantum Programs_.  In: 9th Workshop on Evaluation and Usability of Programming Languages and Tools (PLATEAU ’18), pp. 4:1–4:14, 10.4230/OASIcs.PLATEAU.2018.4.
  * [HM19] Y. Huang & M. Martonosi (2019): _Statistical Assertions for Validating Patterns and Finding Bugs in Quantum Programs_.  In: Proc. ISCA ’19, pp. 541–553, 10.1145/3307650.3322213.  Available at https://arxiv.org/abs/1905.09721.
  * [Hoa69] C.A.R. Hoare (1969): _An Axiomatic Basis for Computer Programming_.  Commun. ACM 12(10), pp. 576–580, 10.1145/363235.363259.
  * [JP09] P. Jorrand & S. Perdrix (2009): _Abstract Interpretation Techniques for Quantum Computation_.  In: Semantic Techniques in Quantum Computation, Cambridge University Press, pp. 206–234, 10.1017/CBO9781139193313.007.
  * [MN19] A. Matuschak & M.A. Nielsen (2019): _Quantum Computing for the Very Curious (and other essays)_.  Online.  Available at https://quantum.country.
  * [NC10] M.A. Nielsen & I.L. Chuang (2010): _Quantum Computation and Quantum Information_ , 10th anniversary edition.  Cambridge University Press, 10.1017/CBO9780511976667.
  * [NMB08] A. Nanevski, G. Morrisett & L. Birkedal (2008): _Hoare Type Theory, Polymorphism and Separation_.  J. Funct. Program. 18(5–6), pp. 865–911, 10.1017/S0956796808006953.  Available at https://software.imdea.org/~aleks/htt/jfpsep07.pdf.
  * [NMS+08] A. Nanevski, G. Morrisett, A. Shinnar, P. Govereau & L. Birkedal (2008): _Ynot: Dependent Types for Imperative Programs_.  In: Proc. ICFP ’08, pp. 229–240, 10.1145/1411204.1411237.  Available at https://software.imdea.org/~aleks/htt/ynot08.pdf.
  * [Per12] J.T. Perconti (2012): _Hoare Type Theory: Dependent Types for State_.  Available at http://www.ccs.neu.edu/home/amal/course/7480-s12/HTT-notes.pdf.
  * [PRZ17] J. Paykin, R. Rand & S. Zdancewic (2017): _QWIRE: A Core Language for Quantum Circuits_.  In: Proc. POPL ’17, pp. 846–858, 10.1145/3009837.3009894.  Available at https://jpaykin.github.io/papers/prz_qwire_2017.pdf.
  * [Ran18] R. Rand (2018): _Formally Verified Quantum Programming_.  Ph.D. thesis, University of Pennsylvania.  Available at https://repository.upenn.edu/edissertations/3175.
  * [Rey02] J.C. Reynolds (2002): _Separation Logic: A Logic for Shared Mutable Data Structures_.  In: Proc. LICS ’02, pp. 55–74, 10.1109/LICS.2002.1029817.  Available at https://www.cs.cmu.edu/~jcr/seplogic.pdf.
  * [Ros15] N.J. Ross (2015): _Algebraic and Logical Methods in Quantum Computation_.  Ph.D. thesis, Dalhousie University.  Available at http://arxiv.org/abs/1510.02198.
  * [RPZ18] R. Rand, J. Paykin & S. Zdancewic (2018): _QWIRE Practice: Formal Verification of Quantum Circuits in Coq_.  In: Proc. QPL ’17, pp. 119–132, 10.4204/EPTCS.266.8.
  * [RS17] F. Rios & P. Selinger (2017): _A Categorical Model for a Quantum Circuit Description Language (Extended Abstract)_.  In: Proc. QPL ’17, pp. 164–178, 10.4204/EPTCS.266.11.
  * [RSSL19] R. Rand, A. Sundaram, K. Singhal & B. Lackey (2019): _A Type System for Quantum Resources_.  Available at http://ks.cs.uchicago.edu/publication/quantum-resource-types/.  Draft.
  * [RSSL20] R. Rand, A. Sundaram, K. Singhal & B. Lackey (2020): _Gottesman Types for Quantum Programs_.  In: Proc. QPL ’20, this volume of EPTCS, Open Publishing Association.
  * [Sel20] P. Selinger (2020): _Dependently Typed Quantum Programming in Proto-Quipper_.  Available at https://popl20.sigplan.org/details/planqc-2020-papers/15/.  Invited talk at PLanQC ’20.
  * [Sin20] K. Singhal (2020): _Quantum Hoare Type Theory_.  Technical Report, University of Chicago.  Available at https://arxiv.org/abs/2012.02154.  Master’s paper.
  * [Unr19] D. Unruh (2019): _Quantum Hoare Logic with Ghost Variables_.  In: Proc. LICS ’19, pp. 1–13, 10.1109/LICS.2019.8785779.  Available at https://arxiv.org/abs/1902.00325.
  * [VAS06] J. Vizzotto, T. Altenkirch & A. Sabry (2006): _Structuring Quantum Effects: Superoperators as Arrows_.  Math. Struct. Comput. Sci. 16(3), pp. 453–468, 10.1017/S0960129506005287.  Available at https://arxiv.org/abs/quant-ph/0501151.
  * [Yin12] M. Ying (2012): _Floyd–Hoare Logic for Quantum Programs_.  ACM Trans. Program. Lang. Syst. 33(6):19, 10.1145/2049706.2049708.

##  Appendix A Modular Bell Pair

Writing a modular version of testBell.

###  A.1 Hadamard basis states

⬇

qplus : {emp} r : Qbit {Id(r, $|$$+$\$>$)}

$=$ do q $\Leftarrow$ mkQbit false;

applyU (H q);

return q

qminus : {emp} r : Qbit {Id(r, $|$-\$>$)}

$=$ do q $\Leftarrow$ mkQbit true;

applyU (H q);

return q

###  A.2 Creating entanglement

⬇

share : \Pi a : Qbit.

{a \in {$|$$+$\$>$, $|$-\$>$}} \-- a \in {b, c} is short for Id(a,b) \$/$
Id(a,c)

b : Qbit

{Id(a, b) $/$\ a \in {$|$0\$>$, $|$1\$>$}}

$=$ \a.do b $\Leftarrow$ mkQbit false;

applyU (ifQ a (X b));

return b

###  A.3 Bell pair

⬇

bell : {emp} (a, b) : (Qbit, Qbit) {Id(a, b) $/$\ a \in {$|$0\$>$, $|$1\$>$}}

$=$ do qa $\leftarrow$ qplus;

qb $\leftarrow$ share qa;

return (qa, qb)

###  A.4 Testing modular Bell pair

⬇

testBell : {emp} (a, b) : (Bool, Bool) {emp $/$\ Id(a, b)}

$=$ do (qa, qb) $\leftarrow$ bell;

(measQbit qa, measQbit qb)

##  Appendix B Teleportation

###  B.1 Alice’s circuit

⬇

alice : \Pi a : Qbit. \Pi e : Qbit

{(Id(a, -) $/$\ entangled(e)}

r : (Bool, Bool)

{emp}

$=$ \a.\e.do applyU (ifQ a (X e));

applyU (H a);

(measQbit a, measQbit e)

###  B.2 Bob’s circuit

⬇

bob : \Pi m1 m2 : Bool. \Pi e : Qbit

{entangled(e)}

r : Qbit

{(Id(r, -)}

$=$ \m1.\m2.\e.do if m1 then applyU (Z e) else ();

if m2 then applyU (X e) else ();

return e

###  B.3 Teleport

⬇

teleport : \Pi q : Qbit. x : Pure.

{Id(q, x)}

r : Qbit

{Id(r, x)}

$=$ \q.do (a, b) $\leftarrow$ bell;

(m1, m2) $\leftarrow$ alice q a;

tq $\leftarrow$ bob m1 m2 b;

return tq