# The First Photometric Study of AH Mic Contact Binary System

A. Poro  Astronomy Department of the Raderon Lab., BC, Burnaby, Canada  M.G.
Blackford  Variable Stars South (VSS), Congarinni Observatory, Congarinni,
NSW, 2447, Australia  S. Ranjbar Salehian  Binary Systems of South and North
(BSN) Project, Contact Systems Department, Iran  E. Jahangiri  Binary Systems
of South and North (BSN) Project, Contact Systems Department, Iran  M. Samiei
Dastjerdi  Binary Systems of South and North (BSN) Project, Contact Systems
Department, Iran  M. Gozarandi  Binary Systems of South and North (BSN)
Project, Contact Systems Department, Iran  R. Karimi  Binary Systems of South
and North (BSN) Project, Contact Systems Department, Iran  Department of
Physics, Kharazmi University, Tehran, Iran  T. Madayen  Astronomy and
Astrophysics Department, University of Toronto, Toronto, Canada  E. Bakhshi
Binary Systems of South and North (BSN) Project, Contact Systems Department,
Global Online Project, Iran  Department of Physics, Institute for Advanced
Studies in Basic Sciences (IASBS), Zanjan, Iran  F. Hedayati  Binary Systems
of South and North (BSN) Project, Contact Systems Department, Iran

###### Abstract

The first multi-color light curve analysis of the AH Mic binary system is
presented. This system has very few past observations from the southern
hemisphere. We extracted the minima times from the light curves based on the
Markov Chain Monte Carlo (MCMC) approach and obtained a new ephemeris. To
provide modern photometric light curve solutions, we used the Physics of
Eclipsing Binaries (Phoebe) software package and the MCMC approach. Light
curve solutions yielded a system temperature ratio of $0.950$, and we assumed
a cold star-spot for the hotter star based on the O’Connell effect. This
analysis reveals that AH Mic is a W-subtype W UMa contact system with a fill-
out factor of 21.3% and a mass ratio of 2.32. The absolute physical parameters
of the components are estimated by using the Gaia Early Data Release 3 (EDR3)
parallax method to be $M_{h}(M_{\odot})=0.702(26)$,
$M_{c}(M_{\odot})=1.629(104)$, $R_{h}(R_{\odot})=0.852(21)$,
$R_{c}(R_{\odot})=1.240(28)$, $L_{h}(L_{\odot})=0.618(3)$ and
$L_{c}(L_{\odot})=1.067(7)$. The orbital angular momentum of the AH Mic binary
system was found to be $51.866(35)$. The components’ positions of this system
are plotted in the Hertzsprung-Russell (H-R) diagram.

Techniques: photometric, Stars: binaries: eclipsing, Stars: individual: AH Mic

##  1 Introduction

Each component of a contact binary system has its own Roche lobe, with the
inner and outer critical surfaces sharing a common envelope. These stars are
generally late type F, G, or K type stars. The majority of contact systems are
solar-type dwarfs with active stellar populations, as evidenced by their
asymmetric light curves and spectral line profiles.

Contact binaries typically have orbital periods of less than one day, and the
Zhang & Qian (2020) investigation found that the lower limit of orbital period
is 0.15 days. The orbital period is frequently varied over time, due to mass
transfer from one component to the other, as well as mass and angular momentum
loss from the binary system. Although, according to Latković et al. (2021)
study, systems with orbital periods of more than 0.5 days and temperatures
greater than 7000 K probably have radiative envelopes and should not be
classified as W UMa-type contact binaries.

Contact binaries are debatable in terms of mass, and the lower limit of mass
ratio for these systems is presently under investigation (e.g. Yang & Qian
2015; Arbutina 2007; Li & Zhang 2006). Additionally, since many investigations
lack a spectroscopic method, determining the mass ratio is a challenging
issue.

W Ursae Majoris (W UMa)-type contact systems (EW) are classified into two
categories, W-type and A-type, based on their properties. The mass ratio in
A-subtype systems is often less than 0.5, showing weak or moderate activity.
In W-subtype systems the less massive component is hotter and the period
changes continuously over time (Binnendijk 1970). In comparison with the
W-subtypes, A-subtype systems are earlier spectral types with higher mass and
luminosity.

The purpose of this study is to present the first light curve analysis of the
AH Mic binary system and determine its properties. AH Mic is a binary system
from the southern hemisphere with a range of apparent magnitude of
$V=12.85-13.45$ and an orbital period of $0.3243344$ days (Watson et al.
2006). In all the available catalogs, this binary system is classified as a
EW-type system. This paper is arranged as follows: Section 2 contains
information about the multi-color photometric observations from an Australian
observatory as well as a data reduction technique. Section 3 explains how to
extract our times of minima, collect minima collection times from literature,
and compute a new ephemeris. The Phoebe Python code and the MCMC approach were
used to study the light curves of the AH Mic system, and the absolute
parameters were obtained in section 4. Section 5 has a summary and conclusion.

##  2 Observation and data reduction

AH Mic (RA. $21$ $04$ $57.8280$, Dec. $-40$ $33$ $06.084$ (J2000)) was
observed for six nights in August and September 2020 with a GSO 14-inch f/8
Ritchey Chretien telescope at the Congarinni Observatory in Australia
($152^{\circ}$ $52^{\prime}$ East and $30^{\circ}$ $44^{\prime}$ South). The
data was taken using a SBIG STT3200-ME CCD camera ($2184\times 1510$ pixels,
6.8 micron square) with $2\times 2$ binning and a CCD temperature of
$-15^{\circ}$C. These observations were made using the Astrodon Johnson-
Cousins $BVRI$ standard filters. The exposure time of these observations for
$BI$ filters is 120 seconds and for $VR$ filters is 60 seconds. GSC 7969-723
was selected as a comparison star and GSC 7969-1170 was chosen as a check star
with an appropriate apparent magnitude in comparison to AH Mic. The comparison
star was found at R.A. 21h 04m 49.91976s (J2000), Dec. $-40^{\circ}$
$39^{\prime}$ $42.8868"$ (J2000) with a $V=12.831$ magnitude, while the check
star was located at R.A. 21h 04m 50.96544s (J2000), Dec. $-40^{\circ}$
$37^{\prime}$ $18.0120"$ (J2000) with a $V=13.285$ magnitude, based on the
AAVSO Photometric All Sky Survey DR9 (APASS9) catalog. Figure 1 shows the
observed field-of-view for AH Mic with the comparison and the check stars.

During the observations, a total of 1417 images were acquired. The CCD image
processing was done with MaxIm DL software, which included dark, bias, and
flat-field corrections (George 2000). The airmass decreases the flux of stars
when light passes through the Earth’s atmosphere. We estimated the airmass
based on observatory coordinates and the position of the star in the sky
during the observations (Hiltner 1962). Finally, we used the AstroImageJ (AIJ)
software to normalize all of the ground-based data (Collins et al. 2017), and
apply airmass111This study’s data (after data reduction) is available in four
filters as a supplement file..

Figure 1: AH Mic, comparison star, and check star field-of-view. The circles
are solely used to show the position of each star, which is indicated by
bigger rings.

##  3 New ephemeris calculation

We extracted the nine times of minima in the light curves which four of them
are the primary minima. We used models based on the Gaussian and Cauchy
distributions to find the times of minima and the MCMC sampling methods to
estimate the uncertainty of the values (Poro et al. 2021). The Python code for
this extraction is implemented using the PyMC3 package (Salvatier et al.
2016). All times of minimum are expressed in the Barycentric Julian Date in
the Barycentric Dynamical Time ($BJD_{TDB}$), and table 1 shows all CCD times
of minima with their uncertainties of AH Mic system.

O-C variations are the differences between observed mid-eclipse times (O) and
their calculated values (C) based on a reference ephemeris. An observable
trend in O-C is caused by various factors and effects. Due to the small number
of minimum times and the short interval of observation (2013–2020), only a
linear fit can be considered for the AH Mic system. We fitted a line on all
times of minima to calculate a new ephemeris based on the MCMC approach using
the emcee package in Python (Foreman-Mackey et al. 2013) as shown in Figure 2.

We determined a new ephemeris for primary minimum as follows:

$BJD_{TDB}(MinI.)=(2456523.00559\pm 0.00068)+(0.32433351\pm 0.00000009)\times
E.$ (1)

Table 1: Available CCD times of minima for AH Mic. $Min.(BJD_{TDB})$ | Error | Epoch | O-C | Reference
---|---|---|---|---
$2456523.00509$ | $0.00040$ | $0$ | $0$ | Diethelm (2014)
$2457154.80792$ | $0.00027$ | $1948$ | $-0.0006$ | Jayasinghe et al. (2019)
$2459108.91600$ | $0.00017$ | $7973$ | $-0.0073$ | This study
$2459089.94354$ | $0.00010$ | $7914.5$ | $-0.0062$ | This study
$2459090.10578$ | $0.00007$ | $7915$ | $-0.0061$ | This study
$2459090.91664$ | $0.00008$ | $7917.5$ | $-0.0061$ | This study
$2459091.07740$ | $0.00013$ | $7918$ | $-0.0075$ | This study
$2459091.24157$ | $0.00016$ | $7918.5$ | $-0.0055$ | This study
$2459092.05120$ | $0.00050$ | $7921$ | $-0.0067$ | This study
$2459092.21189$ | $0.00092$ | $7921.5$ | $-0.0081$ | This study
$2459109.07933$ | $0.00020$ | $7973.5$ | $-0.0061$ | This study
Figure 2: The O-C diagram of AH Mic system.

##  4 Photometric Solutions

###  4.1 Light curve analysis

We utilized the latest Python version of the Phoebe software package (2.3.58)
to model the light curves of the AH Mic binary system. We employed the Phoebe
along with the emcee package (Foreman-Mackey et al. 2013), with the MCMC
approach for more precision in modeling, output results, and uncertainties. We
applied 32 walkers and 1000 iterations to each walker in the MCMC approach.
The Phoebe and MCMC approach then searched for the temperature ratio, orbital
inclination, and mass ratio. The bolometric albedo and gravity-darkening
coefficients were assumed to be $g_{h}=g_{c}=0.32$ (Lucy 1967) and
$A_{h}=A_{c}=0.5$ (Ruciński 1969). The stellar atmosphere was modeled using
the Castelli & Kurucz (2004) method, and the limb-darkening coefficients were
employed as a free parameter in the Phoebe. Different methods were used to
determine the effective temperature for the hotter component as a fixed
parameter. (1) The temperature of the Gaia
DR2222https://gea.esac.esa.int/archive/ for AH Mic system is
$5381_{\rm-(51)}^{+(74)}$ K; (2) The result of the $B-V$ index obtained from
this work’s light curves to be $0.73\pm 0.04$ after calibration (Høg et al.
2000) and as a result we found the effective temperature
$5470_{\rm-(119)}^{+(111)}$ K according to Flower (1996) study; (3) The
effective temperature from the result of Poro et al. (2022) analysis of the
relationship between orbital period and the primary (hotter) temperature for
contact systems is $5681_{\rm-(80)}^{+(81)}$ K. We chose Gaia DR2’s
temperature as a fixed parameter for the hotter star. In the continuation, we
performed the initial light curve analysis. Then, based on the following
equations from Kjurkchieva et al. (2018) study, we calculated the temperatures
and performed the final photometric light curve solutions:

$T_{1}=T_{m}+\frac{c\Delta T}{c+1}$ (2)

$T_{2}=T_{1}-\Delta T$ (3)

where $\Delta T$ is the temperature difference between the two stars, and $c$
is $l_{2}$ divided by $l_{1}$.

The $BVRI$ light curves of the AH Mic system are asymmetric in the maxima.
This shows that the O’Connell effect (O’Connell 1951) is present in this
system, as $MaxII$ is brighter than $MaxI$. We employed a stellar spot model
during the light curve solution in all used filters, and we found that
assuming a cool star-spot model on the hotter component leads to the
appropriate solutions for all $BVRI$ light curves (Table 2). Table 3 shows the
characteristic properties of the AH Mic light curves, with the difference in
maxima in each filter displayed in the first row. Therefore, the $B$ filter
had the highest difference between the levels of maximum light and the depths
of the minima, while the $R$ filter had the least. This is expected because of
the induced changes in the magnetic activity.

Table 2 shows the results of the light curve analysis, whereas Figure 3
displays the observational and theoretical light curves of our modeling.
Figure 4 depicts the 3D view and geometrical structure of the AH Mic binary
system.

Table 2: Photometric solutions of AH Mic. Parameter | Result
---|---
$T_{h}$ (K) | $5550_{\rm-(51)}^{+(74)}$
$T_{c}$ (K) | $5274_{\rm-(63)}^{+(40)}$
$q$ | $2.32_{\rm-(4)}^{+(8)}$
$\Omega_{h}=\Omega_{c}$ | $5.57_{\rm-(13)}^{+(48)}$
$i^{\circ}$ | $75.91_{\rm-(35)}^{+(79)}$
$f$ | $0.213_{\rm-(12)}^{+(35)}$
$l_{h}/l_{tot}$ | $0.381(2)$
$l_{c}/l_{tot}$ | $0.619(2)$
$r_{h(mean)}$ | $0.324(2)$
$r_{c(mean)}$ | $0.470(2)$
Phase shift | $-0.005(1)$
$Colatitude(deg)$ | $90(2)$
$Longitude(deg)$ | $274(1)$
$Radius(deg)$ | $20(1)$
$T_{spot}/T_{star}$ | $0.91(2)$
Table 3: Characteristic parameters of the observational light curves in BVRI filters. Light curve | $\Delta B^{(mag.)}$ | $\Delta V^{(mag.)}$ | $\Delta R^{(mag.)}$ | $\Delta I^{(mag.)}$
---|---|---|---|---
$MaxI-MaxII$ | 0.021 | 0.020 | 0.015 |
$MaxI-MinII$ | -0.479 | -0.494 | -0.472 | -0.468
$MaxI-MinI$ | -0.618 | -0.613 | -0.574 | -0.552
$MinI-MinII$ | 0.139 | 0.119 | 0.102 | 0.084
Figure 3: The observed light curves of AH Mic (black dots), and synthetic
light curves obtained from light curve solutions and residuals are plotted,
with respect to orbital phase, shifted arbitrarily in the relative flux.
Figure 4: The positions of the components of AH Mic.

###  4.2 Absolute parameters

The Gaia EDR3 parallax method can be used to estimate the absolute parameters
when only photometric data is available. Several elements from observational
and light curve analysis are used in the computations, despite the fact that
the Gaia EDR3 parallax has high accuracy for this estimation. Recent
investigations have made some corrections to the Gaia EDR3 parallax (e.g.
Lindegren et al. 2021; Ren et al. 2021). The range of these corrections
appears to include the value of observational and light curve solution
uncertainties (Poro et al. 2022).

Therefore, in order to properly estimate the absolute parameters of the stars
of AH Mic binary system, we used the Gaia EDR3 parallax method. This is a
well-known method, which is fully described in the Poro et al. (2022) study.

To perform the estimation, we first obtained the absolute magnitude ($M_{v}$)
of the system using the maximum light of apparent magnitude $V_{max}=12.79(5)$
from our observation, the star’s distance from Gaia EDR3
$d_{(pc)}=470.2(6.5)$, and the extinction coefficient $A_{v}=0.08$ from the
Schlafly & Finkbeiner (2011) study. Then, by using the value of
$l_{h}/l_{tot}$ and $l_{c}/l_{tot}$ from the light curve solutions, we
obtained the $M_{v}$ value of each of the component. In addition, the
bolometric magnitude $M_{bol}$ value of each star were obtained using
$BC_{h}=-0.133$ and $BC_{c}=-0.199$. Furthermore the luminosity ($L$) values
were estimated and the radius values of each star were obtained using the
temperature of the stars from the light curve solutions. The semi-major axis
of the system was also determined by utilizing $a=R/r_{(mean)}$ equation.
According to Kepler’s third law, the value of $M_{h}+M_{c}$ could be estimated
using the semi-major axis and the orbital period of the system. Moreover, the
mass of the hotter and the cooler stars were determined using the value of
mass ratio.

To estimate the uncertainties using mass ratio and temperature from light
curve analysis, we used the average value of uncertainties rather than the
upper and lower ranges. The results of estimating the absolute parameters for
the AH Mic system are shown in Table 4.

For comparison, we calculated absolute magnitude of system using the orbital
period–luminosity calibration of Rucinski (2006):

$M_{V}=(-1.5)-12logP$ (4)

The result of relation 4 shows $M_{V}^{(mag.)}=4.37$, while in our estimation
$M_{V}^{(mag.)}=4.349(10)$ which show relatively close values.

Table 4: Estimation of absolute parameters of AH Mic binary system. Parameter | Hotter star | Cooler star
---|---|---
$M(M_{\odot})$ | 0.702(26) | 1.629(104)
$R(R_{\odot})$ | 0.852(21) | 1.240(28)
$L(L_{\odot})$ | 0.618(3) | 1.067(7)
$M_{bol}(mag.)$ | 5.263(5) | 4.670(7)
$log(g)(cgs)$ | 4.423(5) | 4.463(8)
$a(R_{\odot})$ | 2.634(48) |

##  5 Summery and conclusion

Observations were made on the $BVRI$ filters for the AH Mic binary system in
the Southern Hemisphere. The MCMC approach was used to extract the times of
minima of our light curves, and the times of minima were also collected from
previous observations. Accordingly, we presented a new ephemeris. The first
light curve analysis of the AH Mic system was done in this work, which was
performed using the latest version of the Phoebe Python code and the MCMC
approach. To estimate the system’s absolute parameters, we used the Gaia EDR3
parallax method, and the rest of the parameters were calculated using the
observational parameters and the light curve solutions from this study.

We fixed the Gaia EDR3 temperature for the hotter component at first, and
obtained the cooler component temperature from light curve solutions. This
temperature is from the combination of both stars in the system. So, using
Equation 2, we obtained the final temperature values of each star. The
temperature difference between the two components is 276 K. According to
Allen’s table (Cox 2000), these temperatures indicate that the spectral types
of the hotter and cooler components are G5 and G8 respectively.

The mass ratio is an important parameter for contact binaries and it is
determined through $q$-search when there is just photometric observations. The
mass ratio in this study was determined using a method based on the MCMC
approach. The reliable value of mass ratio through photometric light curve
analysis is still under many investigations to approach the accuracy of this
value using the spectroscopy method. However, investigations have been
conducted, one of which is based on the Machine Learning methods using the
Multi-Layer Perceptron (MLP) regression model and was carried out by the Poro
et al. (2022) study. According to this method, we find the estimated range of
mass ratio for AH Mic to be $0.408$ to $0.572$ with a mean of $0.485$, which
includes our calculated mass ratio of $0.431$ from the light curve solutions.
It’s worth mentioning that the $q$ specified in the Poro et al. (2022) study
is the ratio of the less massive star mass to the more massive star mass,
hence $q$ is never more than unity.

The Orbital Angular Momentum ($J_{0}$) of the AH Mic system was calculated,
and its location is displayed on the $logJ_{0}-logM$ diagram (Figure 5). The
value of $logJ_{0}$ of AH Mic was obtained to be $51.866\pm 0.035$. The
diagram shows that this system is in a contact binary systems region.

Figure 5: The location of AH Mic on the $logJ_{0}-logM$ diagram. Detached and
contact binary systems are separated by a quadratic border line. The quadratic
line is based on a study by Eker et al. (2006).

The H-R diagram depicted the evolutionary state of the components (Figure 6).
The cooler star (which is more massive) has already evolved away from the
terminal-age main sequence (TAMS), and the hotter star stays in the main
sequence belt and above the zero-age main sequence (ZAMS). According to the
value of the mass ratio, the fill-out factor, inclination, and location of
components on the H-R diagram, AH Mic follows the general pattern of the
W-type W UMa systems.

Figure 6: The hotter and cooler stars of the AH Mic system are shown on the
H-R diagram.

## Acknowledgments

This manuscript was prepared based on the Binary Systems of the South and
North Project ($http://bsnp.info$). The goal of the project is to investigate
contact binary systems in the northern and southern hemispheres using
observatories from throughout the world.

We utilized the latest version of Gaia EDR3 (http://www.cosmos.esa.int/gaia)
data from the European Space Agency (ESA) mission Gaia in this study. This
work has made use of the SIMBAD database, operated by CDS, Strasbourg, France.
Also, we used the latest version of the Phoebe related to this project website
($http://phoebe-project.org/$). PHOEBE is funded in part by the National
Science Foundation (NSF 1517474, 1909109) and the National Aeronautics and
Space Administration (NASA 17-ADAP17-68). We thank Fatemeh Davoudi for all her
scientific assistance.

## ORCID iDs

Atila Poro: https://orcid.org/0000-0002-0196-9732
Mark G. Blackford: https://orcid.org/0000-0003-0524-2204
Selda Ranjbar Salehian: https://orcid.org/0000-0002-5223-1332
Esfandiar Jahangiri: https://orcid.org/0000-0002-1576-798X
Meysam Samiei Dastjerdi: https://orcid.org/0000-0002-7502-3907
Mohammadjavad Gozarandi: https://orcid.org/0000-0002-5131-6801
Reihaneh Karimi: https://orcid.org/0000-0003-2295-0499
Tabassom Madayen: https://orcid.org/0000-0002-2653-3281
Elnaz Bakhshi https://orcid.org/0000-0002-8782-6045
Farnad Hedayati: https://orcid.org/0000-0003-0925-6806

## References

  * Arbutina (2007) Arbutina, B. 2007, MNRAS, 377, 1635, doi: 10.1111/j.1365-2966.2007.11723.x
  * Binnendijk (1970) Binnendijk, L. 1970, Vistas in Astronomy, 12, 217, doi: 10.1016/0083-6656(70)90041-3
  * Castelli & Kurucz (2004) Castelli, F., & Kurucz, R. L. 2004, A&A, 419, 725, doi: 10.1051/0004-6361:20040079
  * Collins et al. (2017) Collins, K. A., Kielkopf, J. F., Stassun, K. G., & Hessman, F. V. 2017, AJ, 153, 77, doi: 10.3847/1538-3881/153/2/77
  * Cox (2000) Cox, A. N. 2000, Allen’s astrophysical quantities
  * Diethelm (2014) Diethelm, R. 2014, Information Bulletin on Variable Stars, 6093, 1
  * Eker et al. (2006) Eker, Z., Demircan, O., Bilir, S., & Karataş, Y. 2006, MNRAS, 373, 1483, doi: 10.1111/j.1365-2966.2006.11073.x
  * Flower (1996) Flower, P. J. 1996, ApJ, 469, 355, doi: 10.1086/177785
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * George (2000) George, D. B. 2000, International Amateur-Professional Photoelectric Photometry Communications, 79, 2
  * Hiltner (1962) Hiltner, W. A. 1962, Astronomical techniques.
  * Høg et al. (2000) Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 357, 367
  * Jayasinghe et al. (2019) Jayasinghe, T., Stanek, K. Z., Kochanek, C. S., et al. 2019, MNRAS, 486, 1907, doi: 10.1093/mnras/stz844
  * Kjurkchieva et al. (2018) Kjurkchieva, D. P., Dimitrov, D. P., Ibryamov, S. I., & Vasileva, D. L. 2018, PASA, 35, e008, doi: 10.1017/pasa.2017.68
  * Latković et al. (2021) Latković, O., Čeki, A., & Lazarević, S. 2021, ApJS, 254, 10, doi: 10.3847/1538-4365/abeb23
  * Li & Zhang (2006) Li, L., & Zhang, F. 2006, MNRAS, 369, 2001, doi: 10.1111/j.1365-2966.2006.10462.x
  * Lindegren et al. (2021) Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4, doi: 10.1051/0004-6361/202039653
  * Lucy (1967) Lucy, L. B. 1967, ZAp, 65, 89
  * O’Connell (1951) O’Connell, D. J. K. 1951, Publications of the Riverview College Observatory, 2, 85
  * Poro et al. (2021) Poro, A., Davoudi, F., Alicavus, F., et al. 2021, Astronomy Letters, 47, 402, doi: 10.1134/S1063773721060050
  * Poro et al. (2022) Poro, A., Sarabi, S., Zamanpour, S., et al. 2022, MNRAS, 510, 5315, doi: 10.1093/mnras/stab3775
  * Ren et al. (2021) Ren, F., Chen, X., Zhang, H., et al. 2021, ApJ, 911, L20, doi: 10.3847/2041-8213/abf359
  * Ruciński (1969) Ruciński, S. M. 1969, Acta Astron., 19, 245
  * Rucinski (2006) Rucinski, S. M. 2006, MNRAS, 368, 1319, doi: 10.1111/j.1365-2966.2006.10207.x
  * Salvatier et al. (2016) Salvatier, J., Wieckia, T. V., & Fonnesbeck, C. 2016, PyMC3: Python probabilistic programming framework.  http://ascl.net/1610.016
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Watson et al. (2006) Watson, C. L., Henden, A. A., & Price, A. 2006, Society for Astronomical Sciences Annual Symposium, 25, 47
  * Yang & Qian (2015) Yang, Y.-G., & Qian, S.-B. 2015, AJ, 150, 69, doi: 10.1088/0004-6256/150/3/69
  * Zhang & Qian (2020) Zhang, X.-D., & Qian, S.-B. 2020, MNRAS, 497, 3493, doi: 10.1093/mnras/staa2166