The Orbit and Dynamical Mass of Polaris: Observations with the CHARA Array

Polaris (α UMi, HR 424, HD 8890, V = 2.02 mag) is the nearest and brightest classical Cepheid. It is a member of a triple system with a resolved eighth magnitude physical companion at a separation of 18''. The Cepheid has been known for many years to be a single-lined spectroscopic binary with a period of about 30 yr (Roemer 1965; Kamper 1996), with components designated Aa (the Cepheid) and Ab. A thorough compilation of RV data has recently been produced by Torres (2023), resulting in a definitive orbit.

While Polaris Aa is a Cepheid, albeit with a small amplitude, it has several characteristics that are unusual.

• 1.  
  It pulsates in the first overtone (Feast & Catchpole 1997).
• 2.  
  It has an unusually rapid period change (Neilson et al. 2012). For fundamental mode pulsators, the period change agrees with expectations of evolution through the instability strip. For overtone pulsators, some additional factor seems to be involved (Evans et al. 2018).
• 3.  
  The pulsation amplitude can vary. Specifically, it decreased from about 1960 to 1990 (Arellano Ferro 1983). However, it did not die out completely (cease pulsation), but in about 2000 began to increase (e.g., Bruntt et al. 2008; Anderson 2019). Anderson, indeed, found the amplitude to be very stable over the course of 7 yr.
• 4.  
  In addition, there was a "glitch" (phase jump) in the O–C period residuals about 1960 (Arellano Ferro 1983; Turner et al. 2005).
• 5.  
  The instability strip crossing in which Polaris is located has been extensively discussed. Based on the most recent distance and the sign of the period change, Evans et al. (2018) concluded it is on the third crossing, although it has been argued that it is on the first crossing (Anderson 2018).
• 6.  
  The possibility that the Cepheid itself may be the product of a merger is discussed by Bond et al. (2018), Evans et al. (2018), and Anderson (2018), based on the fact that the distant component Polaris B is too cool to match the isochrone of 100 Myr that fits the Cepheid, Polaris Aa.
• 7.  
  The single observation of Polaris to detect polarization (Barron et al. 2022) found a much more complicated Stokes V profile than for other Cepheids implying a complex magnetic field. However, the profile is similar to the nonvariable supergiant α Per (Grunhut et al. 2010), which perhaps reflects the very low pulsation amplitude.
• 8.  
  Several studies have searched for additional periodicities in the RVs (Lee et al. 2008; Anderson 2019, and other studies discussed therein). A periodicity of approximately 120 days was found by Lee et al., for instance, which they identify as due to rotation in a star with spots. Two periodicities (40.2 and 60 days) were found by Anderson; 120 days is an integral multiple of these.

These topics are fully discussed and references to all the velocity data sources are provided by Torres (2023). In this discussion these apparently exceptional properties are placed in the context of Cepheids. The times of periastron passage of Polaris coincide with the amplitude changes and also the phase jump at about 1960. In the extensive study of period changes in Cepheids, Csornyei et al. (2022) find a number with phase jumps, as was found by Szabados previously (Szabados 1991, 1989, and references therein). They are often found in short-period overtone pulsators, and it was suggested by Szabados (1992) that they are likely to be found in binaries. This would link a number of the peculiarities related to Polaris to phenomena seen in other Cepheids, and provide a physical explanation. On one hand, this is surprising because the Polaris Aa and Ab system is comparatively wide. Torres estimates that the periastron separation between the two stars is 6.2 au, or approximately 29 times the Cepheid radius. On the other hand, the orbit is eccentric (e = 0.635), providing a variation of conditions around the orbit. The wide orbit would make it surprising for a characteristic such as a phase jump to be produced by simple tidal interaction. However, overtone pulsators are more erratic in their periods than fundamental mode pulsators, presumably for reasons we do not completely understand. Of course, these peculiarities may exist in other stars but have not been identified because they do not have such long and well-covered data sets.

An important parameter in determining properties of Polaris is its distance. Specifically, the determination of the total mass in the system depends on the distance to convert the angular separation to au. The distance has been controversial, partly because of difficulties because Polaris is very bright. Recent distance determinations are summarized by Engle et al. (2018). Polaris B is the wide companion at a distance of 18'' from the Polaris A system. It can be used for distance determination. Bond et al. (2018) argued that Polaris A and B form a gravitationally bound system, which allows us to use Polaris B for distance determination. The distance to Polaris B using Gaia DR3 (Gaia Collaboration et al. 2023) including the Lindegren parallax offset is 136.90 ± 0.34 pc (Lindegren et al. 2021; Gaia Collaboration et al. 2023). The distance will be reviewed when the Gaia DR4 is available, including corrections such as discussed in Khan et al. (2023).

Several previous steps have led up to the measurement of the mass of the Cepheid. The binary system discussed here is Polaris A made up of the Cepheid Polaris Aa and the companion Polaris Ab. Kamper (1996) published an orbit (with additional velocities by Kamper & Fernie 1998). In addition, Wielen et al. (2000) determined the inclination and the position angle of the line of nodes by comparing the proper motion from the Hipparcos satellite with the average long-term ground-based proper motion. The final parameter to determine the mass was the separation of Aa and Ab. This system was resolved with HST (Evans et al. 2008) using the Advanced Camera for Surveys (ACS). This provided the mass of the Cepheid of ${4.5}_{-1.4}^{+2.2}{M}_{\odot }$. This was followed up with three observations between 2007 and 2014 using WFPC2 and WFC3 since the ACS high-resolution channel (HRC) was no longer available (Evans et al. 2018). The astrometry from these observations was less accurate. However, since they cover a quarter of the orbit, the inclination could be determined from them as well as the separation between the components. The present paper discusses observations around the orbit continued with CHARA and a speckle observation from APO.

The sections below contain the CHARA interferometric observations including data reduction and analysis, the diameter and surface imaging of Polaris, APO speckle observations, the orbit fitting, and discussion of the results.

We began our analysis of the interferometric data by fitting the squared visibility amplitudes for the angular diameter of Polaris Aa. The limb darkening was parameterized with the linear function I_λ (μ) = I_λ (1)[1 − u_λ (1 − μ)] (Hanbury Brown et al. 1974), with $\mu =\cos \zeta$, where ζ is the angle between the line of sight and a surface element of the star. We solved for the limb-darkened diameter (θ_LD), limb-darkening coefficient (u_λ ), and visibility scaling factor (V₀) to account for visibility miscalibrations or the presence of extended flux. The squared visibilities of the limb-darkened model are multiplied by V₀ before comparing with the measured values. We first performed a grid search through the three parameters to find the best fit and then followed a Monte Carlo bootstrap approach with 1000 iterations to determine uncertainties. For each bootstrap iteration, we randomly selected measurements with repetition to construct a new set of visibilities with some measurements repeated and some left out. We then applied Gaussian uncertainties to the resampled measurements and fit a limb-darkened model using the IDL MPFIT package (Markwardt 2009) to optimize the fit for each iteration. We adopted the standard deviation of the bootstrap distributions as the uncertainties on each parameter. The angular diameters had been corrected using the wavelength calibration factors discussed in Section 2. To minimize scatter in the observations, the angular diameters were fit to the squared visibilities averaged over the 10 minute observing sets. The results from the diameter fits are presented in Table 2.

We measured a mean diameter of θ_LD = 3.143 ± 0.027 mas and a mean limb-darkening coefficient of u_λ = 0.120 ± 0.043. At the Gaia distance of 136.90 ± 0.34 pc this corresponds to a mean radius of 46.27 ± 0.42 R_⊙. On most nights, the scaling factor V₀ was greater than 1.0, indicating that the V² measurements are higher than the limb-darkened model. This implies a "negative" background flux, indicating that the scaling term likely arises because of a miscalibration of the visibilities, rather than a physical cause like extended flux outside the field of view. The calibration issues discussed in Section 2, as well as uncertainties in the partially resolved calibrator diameters, could have contributed to the miscalibration. The mean diameter is consistent within the uncertainties with the interferometric diameter measured with the FLUOR instrument at the CHARA Array (Mérand et al. 2006). However, the calibration of the MIRC/MIRC-X visibilities was not accurate enough to confirm the presence of the extended circumstellar envelope reported by Mérand et al. (2006).

To investigate whether we see changes in the diameter over the 4 day pulsation cycle, we fit the diameters separately for each configuration observed on each night as shown in Table 2. We adopted the ephemeris of Berdnikov & Pastukhova (1995) for the times of maximum light:

$$\begin{eqnarray*}&&{\mathrm{HJD}}_{\max }=2,428,260.727+3.969251E.\end{eqnarray*}$$

To account for the changing pulsation period, we used the O–C diagram in Figure 6 of Torres (2023) to calculate the times of maximum light nearest to each angular diameter measurement. The columns in Table 2 are the date of the observation, the Heliocentric Julian Date (HJD) –2,400,000, the time of maximum light used to compute the phase, the angular diameter and the limb-darkening coefficient, the scaling factor, the angular diameter from a fixed limb-darkening coefficient u_λ = 0.120, and the number of visibility measurements. We plot the limb-darkened diameters with the fixed limb-darkening coefficient against the pulsation phase in Figure 3. The dashed line represents the mean angular diameter, and the dotted lines show 0.4% variation in diameter expected from pulsation (Moskalik & Gorynya 2005). The standard deviation of the mean diameter is 2.1 times larger than the variation expected from pulsation. Cepheids are expected to reach a minimum angular diameter around phase 0.8, which is roughly consistent with our measurements. However, the variations in the measured diameters could also be impacted by the visibility miscalibrations discussed in Section 2.

Zoom In Zoom Out Reset image size

The limb-darkened diameter model does a reasonable job fitting the squared visibility amplitudes. However, the strong nonzero and non-180^° closure phases (see Figure 1) suggest the presence of asymmetries on the surface of Polaris Aa. Initial modeling of the data indicated that these strong deviations in the closure phases are not produced by the faint binary companion Ab. Modeling with PMOIRED ¹⁸ (Mérand 2022) and a custom-made routine in IDL, indicates that a few starspots on the surface of Polaris Aa can account for most of the closure phase signal. However, the size, location, and contrast of the spots are not well constrained due to the limited (u, v) coverage during the observations.

We also reconstructed images of the stellar surface of Polaris Aa using SQUEEZE ¹⁹ (Baron et al. 2010, 2012), SURFING (e.g., Roettenbacher et al. 2016), and ROTIR ²⁰ (Baron 2020). Both SURFING and ROTIR perform the image reconstructions on the surface of a spheroid, constraining the size and shape of the star, which is important given the limited (u, v) coverage of the observations. ROTIR is written in Julia and uses the OITOOLS libraries, ²¹ a package for visualizing, modeling, and imaging interferometric data. We used the total variation regularizer with a weight of 0.01–0.05 for ROTIR. To prepare the OIFITS files for imaging, we corrected the visibilities for each individual data set using V₀ scaling factors computed from limb-darkened diameter fits; we used the data averaged over the 10 minute observing sets. We merged together the data sets from all configurations on a given night to improve the (u, v) coverage for the image reconstructions. For the data from 2021 April, we merged the data together from the three consecutive nights.

Figure 4 shows the image of Polaris Aa reconstructed with SURFING and ROTIR for the epoch with the most complete (u, v) coverage from UT 2021 April 2, 3, and 4. Both images show a large bright spot to the north of center. We show the SURFING and ROTIR images from four different nights and the corresponding (u, v) coverage in Figure 5. All of the images show spots on the surface; however, some of the surface features could be created by artifacts in the reconstruction process due to gaps in the (u, v) coverage or systematic errors in the visibility calibration. The bottom row in Figure 5 shows reconstructed images from simulated data of a simple limb-darkened model without spots and with the same (u, v) sampling of Polaris to show which features are likely to be artifacts and which are likely to be spots on the surface. The contrast of the surface spots in the real data is larger than the contrast in the simulated limb-darkened disk data. In addition to having much lower contrast, the spots in the simulated data are distributed in point-symmetric patterns around the center of the star. The spots in the real data that do not follow these point-symmetric patterns are more likely to be true features on the stellar surface. Changing the simulated noise in the limb-darkened models slightly changes the shape and orientation of the artifacts in the simulated images; however, the large-scale structure of the patterns remains consistent between iterations with different simulated noise. The mirroring of bright and dark spots about the origin during some of the Polaris epochs is likely an artifact of the reconstruction process; however, these features are indicative of an asymmetry in brightness distribution. The interferometric data and observables extracted from the reconstructed images are shown in Appendix A.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We used the data sets with 30 s integration times to search for the faint companion (Polaris Ab). The shorter integration time reduces time smearing of the companion signal given the separations of 40–85 mas expected from the previous orbit fit. To account for the asymmetries on the surface of Polaris Aa, we first reconstructed images with SQUEEZE using all of the data from a given night. We then used a modified version of the IDL Binary Grid Search Procedure ²² (Schaefer et al. 2016) where we extract visibilities and closure phases directly from the reconstructed image of Polaris Aa and analytically add in the companion. We started the grid search at the expected location of the binary based on the orbit computed by Evans et al. (2018) and searched within ΔR.A. = ± 20 mas and Δdecl. = ±20 mas from the initial values. We applied a correction for bandwidth smearing to the binary fit. The results from the binary grid search are shown in Appendix B.

We detected the Ab companion on five nights; these nights are indicated by a superscript next to the UT date in Table 1. On the nights with nondetections, the data quality was poor (e.g., 2016 November 18) or Polaris Ab was located outside of the interferometric field of view given by λ²/Δλ. The field of view corresponds to 50 mas in the low spectral resolution mode of MIRC and MIRC-X, so with an expected separation of ∼70 mas, we suspect that bandwidth smearing compromised the binary signal on UT 2018 August 27 and 2019 April 9. We were able to recover the detection of the companion on subsequent dates when we used the higher spectral resolution Prism102 and Grism190 modes with a field of view of 100 and 190 mas, respectively.

The average contrast of the companion was ΔH = 7.8 ± 0.4 mag, beyond the typical detection limits of the MIRC and MIRC-X instruments (Gallenne et al. 2015). To improve our confidence in the solutions, we also performed the binary grid search separately for each configuration observed on a given night. For the nights with high confidence solutions, the separations measured from each configuration separately were consistent with the results based on the full night of data. Additionally, the three successive nights in 2021 April produced consistent binary solutions.

We used a bootstrap approach with 100 iterations to estimate the uncertainties based on the 67.5% confidence ellipses for two parameters from the bootstrap distributions. These errors are added in quadrature with the uncertainties from the wavelength calibration. The best-fit binary positions (corrected for wavelength calibration) and flux ratios (f_Ab/f_Aa) are listed in Table 3. The positions of Polaris Ab relative to Aa are given as the separation (ρ) and position angle (PA) measured east of north. We also list the separations projected into R.A. (Δα) and decl. (Δδ). The error ellipses are given by the semimajor axis (σ_maj), semiminor axis (${\sigma }_{\min }$), and the PA of the major axis (σ_PA). The dates are given as HJD and Julian year (JY). In Table 3, we also collected the HST astrometric measurements reported by Evans et al. (2008, 2018) in both polar and Cartesian coordinates. These papers reported the dates in Besselian years, which we have converted to JY.

Table 3. Relative Astrometric Position of the Polaris Ab Companion

HJD	JY	ρ	PA	Δα	Δδ	σmaj	${\sigma }_{\min }$	σPA	fAb/fAa
−2,400,000		(mas)	(deg)	(mas)	(mas)	(mas)	(mas)	(deg)
CHARA
57643.735	2016.6974	33.008	71.696	31.338	10.366	5.757	1.159	18.49	0.00089 ± 0.05838
58728.831	2019.6683	75.537	340.239	−25.538	71.088	0.833	0.547	40.73	0.00072 ± 0.00744
59306.890	2021.2509	84.507	315.900	−58.810	60.686	1.200	0.663	94.49	0.00110 ± 0.00789
59307.982	2021.2539	84.447	315.787	−58.888	60.527	10.736	3.954	123.45	0.00070 ± 0.02214
59308.785	2021.2561	84.232	315.669	−58.861	60.252	0.880	0.804	135.12	0.00045 ± 0.00363
Speckle
60237.913	2023.800	98.8	286.7	−94.6	28.4	5.3	4.2	0.0	...
HST
53585.488	2005.5866	172.4	231.407	−134.75	−107.54	2.1	2.1	0.0	...
53961.416	2006.6158	169.6	226.385	−122.79	−116.99	3.1	3.1	0.0	...
54298.514	2007.5387	180.0	223.000	−122.76	−131.64	20.0	20.0	0.0	...
55153.597	2009.8798	150.0	216.000	−88.17	−121.35	20.0	20.0	0.0	...
56835.302	2014.4841	85.0	175.000	7.41	−84.68	20.0	20.0	0.0	...

Note. Positions in Table 3 are for J2000.

Download table as:  ASCIITypeset image

3.3.1. Confirmation of Binary Positions with CANDID

The orbital analysis depends on the combination of two data sets, the RVs and relative astrometry. The RV data set is very large and made up of measurements from many sources to cover multiple cycles of the 30 yr orbit of Polaris Aa–Ab. We used the complete set of 3659 RVs with pulsation removed that are tabulated in Table 3 of Torres (2023). ²⁴ To those we added an additional 16 Hermes velocities from the program VELOCE I (Anderson et al. 2024). The data were prepared using the template approach detailed in Anderson (2019). Initially, we applied the RV jitter and offset terms in Table 5 of Torres (2023) and the additional Hermes data (jitter is as defined in Torres 2023; the uncertainty for each data set from the orbit fitting). After experimentation and preliminary fits, the offsets and jitter were redetermined as listed in Table 4. The sources and references are as listed in Table 5 of Torres (2023). Successive pairs of columns to the right provide the RV offset and jitter for the original orbit in Torres (2023) and the solution for RVs and CHARA and HST astrometry here.

The astrometry data set has only 11 two-dimensional data points that cover only 16 yr, or half the orbit.

Table 5 lists the orbital parameters from the spectroscopic RV orbit computed by Torres (2023) and Anderson et al. (2024) and a fit to only the visual orbit computed through a Newton Raphson technique using the IDL orbit fitting library ²⁵ described by Schaefer et al. (2016). The orbital parameters include the orbital period P_orb, time of periastron passage T_p, eccentricity e, argument of periastron with reference to the primary ω, primary RV amplitude K₁, systemic velocity γ, longitude of ascending node Ω for J2000, inclination i, and angular semimajor axis a. The eccentricity and ω_Aa differ by 1.4σ and 2.3σ, respectively, between the spectroscopic and visual orbits.

Our methods for fitting the RVs and astrometry simultaneously are described in the following subsections. It was realized early on that the analysis approach is particularly important since the orbit has a substantial eccentricity. In a preliminary solution weighted solely by the errors on the data points in the two data sets, the HST astrometry was systematically offset from the orbit. The eccentricity and ω_Aa are set essentially by the large number of RVs. We explored several approaches to combine the different data sets. The various weighting methods are described in more detail below with the orbital parameters from each solution listed in Table 5.

To determine the orbital parameters, we simultaneously fitted the RVs and astrometric position using a Markov Chain Monte Carlo (MCMC) routine, ²⁶ whose log-likelihood function is given as (Gallenne et al. 2019a)

$$\begin{eqnarray*}&&\mathrm{log}({ \mathcal L })=-\frac{1}{2}\,{\chi }^{2},\mathrm{with}\,{\chi }^{2}={\chi }_{\mathrm{RV}}^{2}+{\chi }_{\mathrm{ast}}^{2}.\end{eqnarray*}$$

In our previous works on binary Cepheid systems with astrometric data (see, e.g., Gallenne et al. 2019b, 2018), both the pulsation and orbital motions of the Cepheids were fitted together with the astrometry. In this work, we decided to ignore the pulsation motion of the Cepheid by directly fitting the pulsation-corrected RVs from Torres (2023). In this case, our RV model is simply defined as

$$\begin{eqnarray*}&&{\chi }_{\mathrm{RV}}^{2}=\sum \frac{{\left({V}_{1}-{V}_{1{\rm{m}}}\right)}^{2}}{{\sigma }_{{{\rm{V}}}_{1}}^{2}},\end{eqnarray*}$$

in which V₁ and σ₁ denote the measured RV and uncertainties for Polaris, and V_1m is the Keplerian velocity model characterized by Heintz (1978):

$$\begin{eqnarray*}\begin{array}{rcl}{V}_{1{\rm{m}}} & = & \gamma +{K}_{1}\,[\cos (\omega +\nu )+e\cos \omega ],\\ \tan \displaystyle \frac{\nu }{2} & = & \sqrt{\displaystyle \frac{1\,+\,e}{1-e}}\tan \displaystyle \frac{E}{2},\\ E-e\sin E & = & \displaystyle \frac{2\pi (t-{T}_{{\rm{p}}})}{{P}_{\mathrm{orb}}},\end{array}\end{eqnarray*}$$

where γ is the systemic velocity, e is the eccentricity, ω is the argument of periastron, ν is the true anomaly, E is the eccentric anomaly, t is the observing date, P_orb is the orbital period, T_p is the time of periastron passage, and K₁ is the primary RV amplitude. The ${\chi }_{\mathrm{ast}}^{2}$ measure for astrometry depends on the astrometric measurements as

$$\begin{eqnarray*}&&\begin{array}{l}{\chi }_{\mathrm{ast}}^{2}={\chi }_{{\rm{a}}}^{2}+{\chi }_{{\rm{b}}}^{2},\\ {\chi }_{{\rm{a}}}^{2}=\sum \frac{{[({\rm{\Delta }}\alpha -{\rm{\Delta }}{\alpha }_{{\rm{m}}})\sin {\sigma }_{\mathrm{PA}}+({\rm{\Delta }}\delta -{\rm{\Delta }}{\delta }_{{\rm{m}}})\cos {\sigma }_{\mathrm{PA}}]}^{2}}{{\sigma }_{\mathrm{maj}}^{2}},\\ {\chi }_{{\rm{b}}}^{2}=\sum \frac{{[-({\rm{\Delta }}\alpha -{\rm{\Delta }}{\alpha }_{{\rm{m}}})\cos {\sigma }_{\mathrm{PA}}+({\rm{\Delta }}\delta -{\rm{\Delta }}{\delta }_{{\rm{m}}})\sin {\sigma }_{\mathrm{PA}}]}^{2}}{{\sigma }_{\min }^{2}},\end{array}\end{eqnarray*}$$

in which $({\rm{\Delta }}\alpha ,{\rm{\Delta }}\delta ,{\sigma }_{\mathrm{PA}},{\sigma }_{\mathrm{maj}},{\sigma }_{\min })$ denote the relative astrometric measurements with the corresponding error ellipses, and (Δα_m, Δδ_m) represent the astrometric model defined as

$$\begin{eqnarray*}&&\begin{array}{l}\,{\rm{\Delta }}{\alpha }_{{\rm{m}}}=r\,[\sin {\rm{\Omega }}\cos (\omega +\nu )+\cos i\cos {\rm{\Omega }}\sin (\omega +\nu )],\\ \,{\rm{\Delta }}{\delta }_{{\rm{m}}}=r\,[\cos {\rm{\Omega }}\cos (\omega +\nu )-\cos i\sin {\rm{\Omega }}\sin (\omega +\nu )],\\ \,r=\displaystyle \frac{a(1-{e}^{2})}{1\,+\,e\cos \nu },\end{array}\end{eqnarray*}$$

where Ω is the longitude of ascending node, i is the orbital inclination, and a is the angular semimajor axis.

As a starting point for our 100 MCMC walkers, we performed a least-squares fit using orbital values from Torres (2023) as first guesses. We then ran 100 initialization steps to well explore the parameter space and get settled into a stationary distribution. For all cases, the chain converged before 50 steps. Finally, we used the last position of the walkers to generate our full production run of 1000 steps, discarding the initial 50 steps. All the orbital elements are estimated from the distribution taking the median value and the maximum value between the 16th and 84th percentiles as uncertainty (although the distributions were roughly symmetrical).

Polaris is a single-lined spectroscopic binary that does not allow the individual masses to be determined independently without additional information. To this end, we adopt the Gaia distance for Polaris B (above). In our MCMC procedure, we included the parallax measurement using a normal distribution centered on 7.3045 mas with a standard deviation of 0.0178 mas. Masses are then derived with the following equations:

$$\begin{eqnarray*}\begin{array}{rcl}{M}_{{\rm{T}}} & = & \displaystyle \frac{{a}^{3}}{{P}_{{\rm{orb}}}^{2}\,{\varpi }^{3}},\\ {M}_{2} & = & \displaystyle \frac{0.03357\,{a}^{2}\,{K}_{1}\,\sqrt{1-{e}^{2}}}{{\varpi }^{2}\,{P}_{{\rm{orb}}}\sin i},\\ {M}_{1} & = & {M}_{{\rm{T}}}-{M}_{2},\end{array}\end{eqnarray*}$$

with M_T = M₁ + M₂ the total mass in M_⊙, a and ϖ in mas, P_orb in yr, and K₁ in km s⁻¹.

Initially, when fitting the orbit simultaneously to the RVs and astrometry, we weighted each measurement by their respective uncertainties. We performed a similar fit weighted by measurement uncertainties using the IDL orbit fitting library, but with uncertainties computed through a bootstrap technique. The orbital parameters from this fit are listed in the "σ Weights" column of Table 5.

To more evenly weight the two data sets, we can compute an average χ² statistic multiplying the χ² from each data set by 1/N, where N is the number of data points (Mérand et al. 2015, Equation (5)). The total χ² for the simultaneous fit then becomes

$$\begin{eqnarray*}&&{\chi }^{2}=\frac{1}{{N}_{\mathrm{RV}}}\,{\chi }_{\mathrm{RV}}^{2}+\frac{1}{{N}_{\mathrm{ast}}}\,{\chi }_{\mathrm{ast}}^{2},\end{eqnarray*}$$

where N_RV is the number of RV measurements and N_ast is the number of astrometry data points. The orbital parameters are listed in the "N σ Weight" column of Table 5. Uncertainties in the orbital parameters were determined from MCMC distributions.

Another weighting scheme is motivated by the Synthetic Minority Over-sampling Technique (SMOTE) used in the machine learning community (Chawla et al. 2002). In order to weight the two sets of measurements more evenly, we replicated each position measurement 167 times to give a roughly even number of observables (N = 167 × 11 × 2 = 3674). We applied random Gaussian uncertainties to the replicated measurements, creating a swarm of simulated observables for each data point. We then fit a simultaneous orbit to the RVs and the replicated positions using the IDL orbit fitting library. ²⁷ We followed a Monte Carlo bootstrap approach to compute uncertainties. For each bootstrap iteration, we randomly selected from the 11 astrometric measurements, with repetition. We then applied Gaussian uncertainties to the resampled measurements and replicated each data point as described above. We repeated this process 1000 times and computed uncertainties in the orbital parameters from the standard deviation of the bootstrap distributions. The orbital parameters from this fit are listed in the "Replicated" column of Table 5. The spectroscopic and astrometric orbits are shown in Figure 6. The orbital parameters and uncertainties are consistent with the 1/N weighting method. The averaged χ² in Section 5.4 underweights the RVs, while the replicated points overweight the astrometry; but both methods provide similar results in the end.

Zoom In Zoom Out Reset image size

In summary, all three weighting schemes produced similar masses within the errors. As stated in Section 5.2, in the solution based on measurement errors (Section 5.3) of the astrometric points from HST ACS/HRC are systematically displaced from the orbit solution. Thus, the astrometry is poorly represented. This is possibly due to systematic errors in the historical RVs, which are difficult to quantify, and which may on their own be distorting the shape of the orbit and imposing that shape on the much less numerous astrometric measurements. For this reason, a solution giving higher weight to the astrometry (Section 5.5) is preferred with a mass of the Cepheid of 5.13 ± 0.28 (5%) M_⊙.

Approximately three-quarters of the 30 yr orbit of Polaris has now been observed with HST, CHARA, and speckle measurements. The resulting mass of the Cepheid is 5.13 ± 0.28 (5%) M_⊙, which is larger than before the CHARA and speckle observations were included, but with smaller errors. The previous mass was 3.45 ± 0.75 M_⊙ (Evans et al. 2018; based on a slightly different distance).

Cepheids provide a quantitative test of evolutionary calculations. In the case of Polaris, the mass derived here is combined with a luminosity to compare with theoretical predictions. For MW Cepheids in general, definitive luminosities will be produced by Gaia in the release that includes orbital fits in the analysis. For Polaris, we use the distance to Polaris B as discussed in Section 1.1. Since Polaris B is not a close binary, an orbit is not required in the Gaia solution.

Figure 7 shows the Cepheid mass–luminosity relation with the location of Polaris in comparison with the predictions from evolutionary calculations of several groups, which predict the luminosity for a given mass.

Zoom In Zoom Out Reset image size

Figure 7 includes V1334 Cyg, which has the most accurate Cepheid mass at present (Gallenne et al. 2018). The solution also provides an accurate (3%) distance to the system. It is a good comparison for Polaris since it has similar characteristics. It pulsates in the first overtone and has a period (luminosity) similar to Polaris. Also shown in Figure 7 are masses for LMC Cepheids in eclipsing binaries from Pilecki et al. (2021). Luminosities for the MW Cepheids are from the distance to Polaris (above) and V1334 Cyg. Luminosities for the LMC Cepheids are from Pilecki et al. (2021). The observed mass–luminosity combinations are compared with evolutionary predictions from Bono et al. (2016) and Anderson et al. (2014). The Bono tracks cover MW and LMC metallicities. For the MW the relations are shown for no main-sequence convective core overshoot and moderate overshoot. The Anderson relation is for main-sequence rotation and moderate convective overshoot. The lower metallicity of the LMC stars results in a higher luminosity than for solar abundance stars, as indicated by the models.

The determination of the mass of Polaris in this study is part of a series of studies to measure that of Cepheids. Additional MW Cepheids are to be added to Figure 7. However, the first indication is that in the examples so far (Polaris, V1334 Cyg, and the LMC Cepheids) the luminosity is larger than predicted by current evolutionary calculations.

To explore the details of the calculation results, we provide Figure 8 updated from Figure 7 in Evans et al. (2018). Temperatures are from (B−V)₀, as discussed in Evans et al. (2018). Figure 8 illustrates the location of Polaris and also V1334 Cyg with respect to evolutionary calculations for 4, 5, and 7 M_⊙ from Georgy et al. (2013). The tracks contrast evolution for stars with zero rotational velocity on the main sequence and those with 0.95 breakup velocity. Both Polaris and V1334 Cyg are more luminous than the 5 M_⊙ track in the blue loop, even for large main-sequence rotation. For a comparison between codes, Figure 8 in Evans et al. (2018) shows the location of Polaris compared with three sets of tracks for a 5 M_⊙ model: Geneva (Georgy et al. 2013), MIST (Choi et al. 2016), and PARSEC ²⁸ (based on the Bressan et al. 2012 group). All three sets of tracks are shown for zero initial main-sequence rotation; for the Geneva and MIST calculations, a track is also shown for substantial rotation. In this figure, the MIST tracks (both with and without substantial main-sequence rotation) approximately coincide with the Georgy tracks with rotation.

Zoom In Zoom Out Reset image size

One further criterion for the comparison between observed masses and evolutionary tracks is that Cepheids are found in a region of the luminosity–temperature diagram (Figure 8) where the blue loops penetrate the instability strip. This was a motivation of the exploration of the effect of rotation in main-sequence progenitors (Anderson et al. 2014), since increasing the main-sequence core convective overshoot increases luminosity at the Cepheid phase, but truncates the blue loops. Evolutionary tracks for additional masses probe this further.

This study provides astrometry from CHARA interferometry and APO speckle observations for the Polaris Aa and Ab system. When combined with previous HST observations and RVs (Torres 2023), a mass of 5.13 ± 0.28 M_⊙ is derived. The observed luminosity for Polaris is brighter than predicted by current evolutionary tracks. The predicted luminosity depends on the rotation of the main-sequence progenitor, but Figure 8 shows that Polaris is at least 0.4 mag brighter than the predicted tracks including a small correction from the 5 M_⊙ tracks to the mass of Polaris (based on a 0.5 critical breakup velocity; Anderson et al. 2014). This study is part of a series of studies of Cepheid masses using interferometry as well as Gaia astrometry to examine these questions.

In the case of Polaris several of its properties are unusual in the context of Cepheids. It pulsates in the first overtone, which may be linked with the rapid period change. Overtone pulsators have more instability in their pulsation periods, which can sometimes be interpreted as rapid period change. The variation in pulsation amplitude may also be related to its pulsation mode. The recent discussion of velocities (Torres 2023) postulates that the "glitch" in pulsation may occur at periastron passage. This would add a new factor to the interpretation of the observations. It has been suggested for other Cepheids with a "phase jump" that they are likely to be binary systems (Csornyei et al. 2022). In summary, while these properties are unusual for Cepheids, they exist in other stars. Thus, Polaris fits in the framework of pulsating supergiants, particularly if orbital motion is included. While these characteristics complicate the interpretation of observations, all are found in stars without abnormal evolution.

The identification of starspots is consistent with several properties of Polaris. It has a very low pulsation amplitude, which sets it apart from full amplitude Cepheids. This may mean that the atmosphere is like that of a nonvariable supergiant, which have often indicators of activity. It is not clear how full amplitude pulsation affects the atmosphere and magnetic field in pulsators, so Polaris is an interesting test case.

Polaris has a couple of other characteristics that would be consistent with magnetically related spots. Polarization has been found (Barron et al. 2022). However, the polarization measurements look more like those of the nonvariable supergiant α Per than other Cepheids (Grunhut et al. 2010), but this may just reflect the very low pulsation amplitude. Polaris has also been detected in X-rays (Evans et al. 2022).

Our identification of starspots opens the prospect of determining a rotation period. It also can explain why additional periodicities have been difficult to find, since the distribution of spots is variable. The long-period (≃120 days) variation found by Lee et al. (2008), for instance, might be a rotation period. Periods of 40 and 60 days found by Anderson (2019) may be related.

Another possible interpretation of surface features on Polaris is convective supergranules as discussed by Schwarzschild (1975).

A final characteristic of the orbit of Polaris is the high eccentricity. This is often found in systems that have undergone three-body interaction, which would be consistent with the Cepheid being a merger product from a former triple system, as suggested by Evans et al. (2018). However, high eccentricity is frequently found in long-period systems (Shetye et al. 2024), so it does not require special conditions.