Polarization position angle standard stars: a reassessment of θ and its variability for seventeen stars based on a decade of observations

ABSTRACT

Observations of polarization position angle (⁠|$\theta$|⁠) standards made from 2014 to 2023 with the High Precision Polarimetric Instrument (HIPPI) and other HIPPI-class polarimeters in both hemispheres are used to investigate their variability. Multiband data were first used to thoroughly recalibrate the instrument performance by bench-marking against carefully selected literature data. A novel co-ordinate difference matrix (CDM) approach – which combines pairs of points – was then used to amalgamate monochromatic (⁠|$g^\prime$| band) observations from many observing runs and re-determine |$\theta$| for 17 standard stars. The CDM algorithm was then integrated into a fitting routine and used to establish the impact of stellar variability on the measured position angle scatter. The approach yields variability detections for stars on long time-scales that appear stable over short runs. The best position angle standards are |$\ell$| Car, o Sco, HD 154445, HD 161056, and |$\iota ^1$| Sco, which are stable to |$\le$|0.123^{|$^\circ$|}. Position angle variability of 0.27–0.82^{|$^\circ$|}, significant at the 3|$\sigma$| level, is found for 5 standards, including the Luminous Blue Variable HD 160529 and all but one of the other B/A-type supergiants (HD 80558, HD 111613, HD 183143, and 55 Cyg), most of which also appear likely to be variable in polarization magnitude (p) – there is no preferred orientation for the polarization in these objects, which are all classified as |$\alpha$| Cygni variables. Despite this we make six key recommendations for observers – relating to data acquisition, processing and reporting – that will allow them to use these standards to achieve |$\lt $| 0.1^{|$^\circ$|} precision in the telescope position angle with similar instrumentation, and allow data sets to be combined more accurately.

1 INTRODUCTION

The 21st century has seen the advent of broadband optical polarimeters capable of a precision of 10 parts-per-million or better. Their development was sparked by the hunt for exoplanet signatures (e.g. Hough et al. 2006; Wiktorowicz & Matthews 2008; Piirola, Berdyugin & Berdyugina 2014; Bailey et al. 2015) but instead lead to the discovery of new and predicted stellar polarigenic mechanisms, such as rapid rotation (Cotton et al. 2017; Bailey et al. 2020; Lewis et al. 2022; Howarth et al. 2023), binary photospheric reflection (Bailey et al. 2019; Cotton et al. 2020), linear polarization from global magnetic fields (Cotton et al. 2017, 2019a), and non-radial pulsations (Cotton et al. 2022a). Precise maps of interstellar polarization close to the Sun are now possible (Cotton et al. 2016; Piirola et al. 2020), and inferences have been made about the nature of hot dust (Marshall et al. 2016), debris discs (Marshall et al. 2020, 2023), and even the heliosphere (Frisch et al. 2022). Higher precision studies of known phenomena are also revealing new details about such diverse topics as asteroids (Wiktorowicz & Nofi 2015), gas entrained between binary stars (Berdyugin et al. 2018), the nature of the interstellar medium (Cotton et al. 2019b), and extreme variable stars (Bailey et al. 2024). Alongside this progress, the dream of detecting and characterizing exoplanet atmospheres with polarimetry remains a live ambition (Bailey et al. 2021; Bott et al. 2022; Wiktorowicz 2024). The development of new instruments continues at pace, both for medium to very large sized telescopes (Wiktorowicz & Nofi 2015; Bailey et al. 2020; Piirola et al. 2021) and even amateur-sized telescopes (Bailey, Cotton & Kedziora-Chudczer 2017; Bailey et al. 2023).

Despite the ground-breaking improvements in instrumental precision, polarimetric observations of objects at increasing distance are naturally affected by the interstellar polarization background. The detection of small polarization signals from distant objects is therefore critically dependant upon the accurate calibration of the polarization position angle – a craft that has not progressed at the same rate. We aim to address this issue here.

Linear polarization is defined either in terms of normalized Stokes parameters |$q=Q/I$| and |$u=U/I$| (typically measured in per cent: |$10^{-2}$|⁠, or parts-per-million, ppm: |$10^{-6}$|⁠), or as polarization magnitude

and position angle

measured North over east, relative to the north celestial pole (⁠|$\theta _0$|⁠), i.e. in the equatorial system. Polarimetric data is almost universally reported in either or both of these co-ordinate frames, but collected in an instrument frame, (⁠|$q_i$|⁠, |$u_i$|⁠), and then rotated according to

and

where |$\theta _t$|⁠, usually called the telescope position angle, is the difference between the instrument reference axis and |$\theta _0$| – which is readily accessible in astrometry but not polarimetry (Van De Kamp 1967; Hsu & Breger 1982).¹ Instead polarimetrists often have to determine |$\theta _t$| by reference to high polarization standard stars (Serkowski 1974a, b). For this purpose, |$\theta _t$| must be re-determined for every observing run (and whenever the equipment is disturbed) to reflect the current condition of the instrument and telescope. It is also a difficult task to perform with precision and accuracy, since the available calibration stars vary with observing location and season. Indeed, there can sometimes be no established standards in the sky bright enough for polarimetry on the smallest telescopes (e.g. the |$\lt $|10-inch telescopes used by Bailey et al. 2023, 2024).

Despite some standards apparently having |$\theta$| determined to 0.2^{|$^\circ$|} accuracy (Hsu & Breger 1982), the accuracy is usually considered to be only 1^{|$^\circ$|} (e.g. Wiktorowicz & Nofi 2015; Bailey et al. 2020). With recent advances, 1^{|$^\circ$|} accuracy is not always good enough for the intended science (e.g. Cotton et al. 2020).

A good high polarization standard has two qualities: (i) it is non-variable (especially in |$\theta$|⁠), and (ii) it has a high polarization relative to its brightness, since position angle uncertainty, |$e_\theta$|⁠, is related to polarization magnitude uncertainty, |$e_p$|⁠, (Serkowski 1968; Hsu & Breger 1982):

where |$\theta$| is in degrees, and |$e_p$| is a function of photon count when not limited by instrumentation or seeing.

Most ordinary stars have little intrinsic polarization. Instead the dominant polarizing mechanism is the interstellar medium (ISM; Hall 1949; Hiltner 1949; Serkowski 1968). As light travels from a star to the observer, it interacts with oblate dust grains within the ISM aligned by large-scale magnetic fields; these act like a wire grid polarizer. The interstellar polarization imparted is dependent on the uniformity of the ISM as well as the quantity of dust on the sight line – and hence, indirectly, on distance. Within about 100 pc of the Sun – i.e. within the Local Hot Bubble – interstellar polarization is imparted at a rate of about 0.2–2.0 ppm pc⁻¹ (Bailey, Lucas & Hough 2010; Cotton et al. 2016, 2017), beyond that it is 20 ppm pc⁻¹ (Behr 1959).

The ISM is assumed to be unchanging on relevant astrophysical time-scales, which leads to choosing standards that are relatively distant and bright. Typically, the best small telescope standards have polarizations of several per cent, have |$m_V \lesssim 6$|⁠, and have parallaxes |$\lt \, \sim$|2–4 mas – these are necessarily some of the most extreme stars. The standards used today were mostly chosen in the 1960s and 1970s (Serkowski 1968, 1974a, b; Serkowski, Mathewson & Ford 1975; Clarke 2010), with much of the work establishing wavelength dependence and refining |$\theta$| taking place from the 1970s to 1990s (Serkowski et al. 1975; Whittet & van Breda 1980; Wilking, Lebofsky & Rieke 1982; Whittet et al. 1992; Wolff, Nordsieck & Nook 1996; Martin, Clayton & Wolff 1999). The most comprehensive modern re-examination of the wavelength dependence of interstellar polarization was provided by Bagnulo et al. (2017), but there are scant recent works² looking at the long term stability of the most important stars.

In the earlier literature there was an important debate about which standards might be variable. Hsu & Breger (1982), Dolan & Tapia (1986), Lupie & Nordsieck (1987), Bastien et al. (1988), and Clemens & Tapia (1990) all, often contrastingly, identified standards they considered to be variable. Of these, the most thorough analyses were performed by Hsu & Breger (1982) and Bastien et al. (1988). However, these works have all been criticized as not statistically rigorous by Naghizadeh-Khouei (1991), who pointed out that in most cases only partial data was presented and the data sets were small. The work of Bastien et al. (1988) was the most comprehensive, yet came in for particular criticism by Clarke & Naghizadeh-Khouei (1994), who in reanalysing their data were convinced of the variability of only one star out of the eleven claimed. There, the main objection was that the data were drawn from different sets without this being properly accounted for, and the reanalysis used only a subset of the observations. Some time later Bastien et al. (2007) revisited their work. They applied the cumulative distribution function (CDF) test ‘in a very conservative manner’ that was used and recommended by Clarke & Naghizadeh-Khouei (1994), concluding that 7 of the 11 stars they originally declared variable were, and that the other 4 ‘may be’. This does not seem a particularly satisfactory resolution. Consequently, a pall hangs over the question of which polarization standards are variable on long time-scales, and the caution implied by Bastien et al. (1988)’s findings has gone substantially unheeded by observers.

Putting aside the controversy, more broadly there are three motivations that provoke further study of these stars:

• Interstellar polarization may not be constant on 50-yr time-scales. On any given sight line there will be many different dust clouds, which are in motion with respect to our standard stars. Significant movement of the clouds would cause the observed value of |$\theta$| to vary over time (Bastien et al. 1988; Clarke 2010).

• Extreme stars are the most likely to have large intrinsic polarizations – intrinsic polarization is more common in stars of B-type and earlier³ and K-type and later (Clarke 2010; Cotton et al. 2016), and in more luminous stars (Dyck & Jennings 1971; Clarke 2010; Lewis et al. 2022). Polarization variability could have a very long period, un-captured by prior shorter duration studies, or be episodic as in the case of Be stars (e.g Carciofi et al. 2007) or LBV stars (Gootkin et al. 2020). So, stars seemingly non-variable decades ago may not be so now.

• Modern high precision polarimeters (Wiktorowicz & Nofi 2015; Bailey et al. 2020, 2023; Piirola et al. 2021 ) are up to 100|$\times$| more precise than those used to establish the standards. Consequently, new stellar polarigenic mechanisms are now being detected (Cotton et al. 2017; Bailey et al. 2019; Cotton et al. 2022b). Yet, polarimetric variation associated with these phenomena is usually small, so its study is limited to the nearest stars – without precise |$\theta$| calibration large interstellar polarization overwhelms small intrinsic signals investigated over many observing runs.

Because interstellar polarization increases with distance, the number of objects that can be studied long term at high precision is severely limited and many rarer stellar types are completely unavailable. To enable the discovery of new polarigenic mechanisms this must be remedied. To understand the scale of the problem, consider polarization due to binary reflection: in the Spica system this has an amplitude of 200 ppm (Bailey et al. 2019) – represented by the solid grey horizontal line in Fig. 1. A 1^{|$^\circ$|} error in |$\theta$| can produce errors in the Stokes parameters at that level at a distance of 300 pc (⁠|$p_{\rm ISM} \sim$| 0.55 per cent). The predicted Rayleigh scattering signal from hot-Jupiter exoplanet atmospheres in the combined light of star and planet is, at best, of order 10–20 ppm (Bott et al. 2016, 2018; Bailey et al. 2021). Similarly, the pulsation-driven polarization produced in the |$\beta$| Cep variable |$\beta$| Cru is just 12 ppm (Cotton et al. 2022b) (dashed grey line). For these signal levels a 1^{|$^\circ$|} error can be significant even within 100 pc of the Sun. Improving precision in |$\theta$| to 0.1^{|$^\circ$|} displaces the threshold for hot-Jupiter or |$\beta$| Cru like polarization to 300 pc, and Spica like polarization to 3000 pc.

Our first objective in this paper is to establish mean |$\theta$| offsets between the standards. As it stands, varying the mix of standards changes the calibration. Presumably, zero-point differences between different observers are a source of imprecision. The second objective is to provide an updated assessment of the position angle variability of established polarization standards – especially long-term variability – and in so doing determine which, if any, are suitable for achieving 0.1^{|$^\circ$|} precision.

This paper is structured as follows: Section 2 provides background on each of the high polarization standard stars studied. Section 3 describes our observations; the analysis of which is carried out in Section 4. In Section 5, we discuss the implications of the results. Of particular note, Section 5.4 shows the impact of each correction we made. While, Section 5.5 lists six specific recommendations for observers relating to the acquisition, processing, and reporting of position angle data. The conclusions are presented in Section 6. Appendices  A,  B, and  C detail literature data and calibration details. For easy reference, Appendix  D lists a selection of symbols used through the paper.

2 HIGH POLARIZATION STANDARD STARS

Very bright high polarization standards are rare. The large distances required for significant interstellar polarization mean that only stars with small absolute magnitudes are bright enough. As a result, most standards trace their lineage to the first decades of stellar polarimetric study when the first bright star surveys were being conducted. In particular, the most used standards are drawn from a recommended list first published by Serkowski (1974a). The parameters for those stars were all refined in Serkowski et al. (1975). Other observers have occasionally added to (or subtracted from) this list, according to their needs, but have largely applied the same selection criteria. There are perhaps as many as two dozen standards in irregular use, depending on what brightness criteria are applied. These stars are far from evenly distributed across the sky. Overwhelmingly, the standards are located in dusty regions fairly close to the Sun, such as the Sco-Cen association; the few that are not can be very important. For instance, Matsumura, Seki & Kawabata (1998) described reports of variability in HD 43 384 as a ‘serious problem’, stressing that there was no bright alternative within |$\sim$|6 h right ascension in the Northern hemisphere.

We have largely worked from southern mid-latitudes, and so most stars we report on here are accessible primarily from there, but the transportation of an instrument to the Monterey Institute for Research in Astronomy (MIRA), has allowed us to add a number of northern stars. The standard stars in this study all appear in the catalogues of Serkowski et al. (1975), Hsu & Breger (1982), and/or Bagnulo et al. (2017); their properties are summarized in Table 1. They are all either well-established standards or have been used as such in making observations with the High Precision Polarimetric Instrument (HIPPI) and other HIPPI-class polarimeters. Appendices  A and  B provide references and describe, in meticulous detail, how we came to favour the tabulated polarization and reddening properties. The co-ordinates and magnitudes for each standard given here – that define which telescopes they are accessible to – are taken directly from SIMBAD. Below is an account of other pertinent details, including variability found by other methods that might portent polarimetric variability, as well as a detailed account of claims and counter-claims of polarimetric variability for each star.

2.1 HD 7927

HD 7927 (⁠|$\phi$| Cas) is a bright yellow supergiant star of spectral type F0 Ia (Gray, Napier & Winkler 2001) that is likely, though not conclusively, a member of the NGC 457 moving group (Eggen 1982; Rosenzweig & Anderson 1993). It has two notable visual companions, the brightest companion (⁠|$\phi ^2$| Cas) is |$m_V=$| 7.04, 132.8 arcsec away, and the closest companion is a |$m_V=$| 12.3 at 48.4 arcsec separation (Mason et al. 2001). Small amplitude variations with no defined period have been found in RV (Adams, Joy & Sanford 1924; Arellano Ferro, Parrao & Giridhar 1988)⁴ and in photometry by Percy (1989), who note that the photometric variations are too small compared to RV to indicate Cepheid-like behaviour.

First measurements of HD 7927’s polarization were made by Hiltner (1951). The star was not found to be variable by Coyne (1972) but he did note its |$p(\lambda)$| as anomalous. No variability was found by Hsu & Breger (1982), whose claimed detection thresholds are 0.01 per cent in p and 0.2^{|$^\circ$|} in |$\theta$|⁠. Wavelength dependence of |$\theta$| in HD 7927 has been observed on multiple occasions (Gehrels & Silvester 1965; Coyne & Gehrels 1966; Hsu & Breger 1982) but only Dolan & Tapia (1986) found that the wavelength dependence varied from night to night; they emphasize this as critically problematic for a position angle standard. Dolan & Tapia (1986) also found |$\theta$| variable. Furthermore, Bastien et al. (1988) found HD 7927 to be variable in both p and |$\theta$|⁠, although the results of this paper are heavily criticized and this result refuted by Clarke & Naghizadeh-Khouei (1994). Earlier Naghizadeh-Khouei (1991) had described this star as displaying ‘definite polarization variability’ both in |$\theta$|⁠, and in p in R band (but not in p in B band) based on his own observations.

2.2 HD 23512

HD 23 512 (BD|$+$|23 524) is an A0 V type star (Fitzpatrick & Massa 2007) and is a member of the Pleiades cluster (Abt & Levato 1978). The star has a companion, discovered by lunar occultation, with a brightness difference of 2 mag and a separation of |$0.1^{\prime \prime }$| (Mason et al. 2001) or |$0.05^{\prime \prime }$| (Torres, Latham & Quinn 2021). It has been a candidate for having a variable RV (Smith & Struve 1944) but this was not confirmed by Abt et al. (1965). The star has also been a double line candidate (Liu, Janes & Bania 1991) but this was not corroborated by Torres (2020). The polarization of HD 23512 was found not to be variable by Hsu & Breger (1982). It is described as ‘clearly’ variable in both p and |$\theta$| by Bastien et al. (1988), which was refuted by Clarke & Naghizadeh-Khouei (1994).

2.3 HD 43384

HD 43384 (9 Gem) is of spectral type B3 Ib (Rachford et al. 2009) classified as an |$\alpha$| Cyg variable star (ESA 1997). Hsu & Breger (1982) found that the star’s polarization angle is variable at a level of 0.8 |$\pm$| 0.2^{|$^\circ$|} on the short term, with larger long term variations apparent (⁠|$\Delta \theta \sim 2^\circ$|⁠; |$\Delta p= 0.25$| per cent over a decade). Coyne (1972) had previously described variability around thrice as much in both p and |$\theta$|⁠. Matsumura et al. (1998) found that the polarization variability (⁠|$\Delta \theta \sim 1^\circ$|⁠; |$\Delta p=0.2~{{\ \rm per\ cent}}$|⁠) was phase locked with the 13.70 d period observed in Hipparcos photometry (ESA 1997). In contrast Dolan & Tapia (1986), though noting an extreme |$\Delta \theta /\Delta \lambda$| found neither that parameter to be complex nor |$\theta$| to be variable.

2.4 HD 80558

HD 80 558 (LR Vel) is a B6 Ia supergiant (Houk 1978) with prominent photometric variability (van Genderen et al. 1989). The polarization of HD 80 558 was first studied by Serkowski & Robertson (1969) and it has been used as a high polarization standard since then. Dolan & Tapia (1986), in comparing their data to Serkowski (1974a)’s, found no significant difference in p or |$\theta$|⁠. Hsu & Breger (1982) also reported no variability. Bastien et al. (1988) found HD 80 558 to have variable polarization over 35 nights of observation. This result was refuted by Clarke & Naghizadeh-Khouei (1994)’s reanalysis of Bastien et al. (1988)’s data.

2.5 HD 84810

HD 84810 (⁠|$\ell$| Car) is a classical Cepheid variable with a spectral type that ranges from F8–G9 (Albrecht 1921) and a period of |$\approx$|35.5 d (Trahin et al. 2021). Owing to its brightness and proximity it has been extensively observed from ultraviolet (UV) to infrared (IR) wavelengths for more than a century (e.g. Bohm-Vitense & Love 1994; Kervella et al. 2006). In principle, the purely radial pulsations of a Cepheid variable should produce no polarization change (Odell 1979), and HD 84810 has been found to be invariable in p and |$\theta$| by Hsu & Breger (1982), Bastien et al. (1988), and Clarke & Naghizadeh-Khouei (1994). Sensitive measurements by Bailey et al. (2023) show only small variations in p of 0.023 |$\pm$| 0.005 per cent from 48 observations over about a year.

2.6 HD 111613

HD 111613 (DS Cru) is a supergiant of spectral type A2 Iab (Ebenbichler et al. 2022) and a member of NGC 4755 (Humphreys 1978). Hsu & Breger (1982) find no variability for HD 111613 in p or |$\theta$|⁠. Dolan & Tapia (1986) saw no change in p over a 4-yr period (1980–84), but found |$\theta$| and its wavelength dependence to be inconsistent between observing runs. Bastien et al. (1988) observed for 41 nights and saw significant variations in both |$\theta$| and p (⁠|$\Delta \theta =2.4^\circ$|⁠, |$\Delta p=0.105$| per cent) on a time-scale of |$\approx$|32 d, a result confirmed by Clarke & Naghizadeh-Khouei (1994)’s reanalysis. Bastien et al. (1988) noted that the polarization was seen to vary slowly, and supposing a binary system, derived an inclination of |$84 \pm 1 ^{\circ }$| based on an assumed a 64-d period.

2.7 HD 147084

HD 147084 (o Sco) is an A4 II bright giant (Martin et al. 1999) in Upper-Scorpius (de Geus, de Zeeuw & Lub 1989). Small amplitude RV variations were measured by Levato et al. (1987) who state that this range may be due to intrinsic motions in the atmosphere. HD 147084 is noteworthy for being a standard for circular as well as linear polarization. It has a maximum fractional circular polarization of approximately 0.04 per cent at |$2.32 ~\mu m$| (Kemp 1972), indicating that the light passes through at least two regions of the interstellar medium with differently aligned dust particles.

HD 147084 has substantial coverage in polarization data spanning ultraviolet to infrared wavelengths, owing to its large |$\lambda _{\rm max}$|⁠, making it particularly useful as a standard (Kemp 1972; Martin et al. 1999). No variability in p or |$\theta$| was found by Hsu & Breger (1982) or Dolan & Tapia (1986). In contrast, Bastien et al. (1988) find it to be variable in both p and |$\theta$|⁠, a result refuted by Clarke & Naghizadeh-Khouei (1994). A small potential variability in p of 0.028 |$\pm$| 0.008 per cent has been found by Bailey et al. (2023) from 108 observations over more than a year.

2.8 HD 149757

HD 149 757 (⁠|$\zeta$| Oph) is a well-studied single star with an O9.5 V spectral type (Hubrig et al. 2011). Its rapid rotation velocity of 400 km s⁻¹ causes it to lose mass through a strong wind (Hubrig et al. 2011), and gives rise to a variable surface temperature through its oblateness (Balona & Dziembowski 1999). Periodic variability for this star has been noted in both photometry and spectroscopy (helium line profiles) consistent with a |$\beta$| Cephei-type classification (Hubrig et al. 2011). The spectral variability is likely the result of non-radial pulsations (Balona & Kambe 1999), where these modes may be excited periodically by lower order modes (Walker et al. 2005). HD 149 757 was one of ten O-type stars included in a study of polarimetric variability by Hayes (1975), from 12 observations over many weeks, he did not find it to be variable. Lupie & Nordsieck (1987) describe the star as non-variable but caution they have few observations. McDavid (2000) carried out a study of nine O-type stars with variable winds, including HD 149757, using agglomerated data from 1949 to 1997; none exhibited statistically significant variability, but small amplitude, short term variability amongst the targets was hinted at by a multi-technique campaign.

2.9 HD 154445

HD 154 445 (HR 6353) is a B1 V spectral type star (Reed 2003); it has no identified companions (Eggleton & Tokovinin 2008). The first reported polarimetric observations of HD 154 445 were made by Hiltner (1951). Repeated observations at optical (e.g. van Panhuys Smith 1956; Serkowski et al. 1975; Cikota et al. 2018) and near-infrared wavelengths (Dyck, Forbes & Shawl 1971; Dyck & Jones 1978) have demonstrated consistency in p and |$\theta$|⁠. Hsu & Breger (1982) find the HD 154445 to be invariable in p and |$\theta$|⁠. The star was claimed as variable in p (but not |$\theta$|⁠) at the 2|$\sigma$| level by Bastien et al. (1988), but this was not borne out by reanalysis (Clarke & Naghizadeh-Khouei 1994). Recently reported observations by Wiktorowicz et al. (2023) present little evidence for variability.

2.10 HD 160529

HD 160 529 (V905 Sco) is a luminous blue variable (LBV) of spectral type A2Ia (Stahl et al. 2003). It has a prolific history as a photometric and spectroscopic variable star (e.g. Wolf, Campusano & Sterken 1974; Sterken et al. 1991). Decades of photometry from Sterken et al. (1991) show that the star’s magnitude dimmed by |$\sim$|0.5 mag over 18 yr. More recent AAVSO data spanning the last 20 yr shows similar time-scales of variability, with as much as a magnitude in brightness changes. A spectroscopic study by Wolf et al. (1974) highlighted many signatures that could likely be attributed to strong mass-loss including line profile variations, line splitting, P-Cygni, and inverse P-Cygni profiles. This large photometric and spectroscopic variability has likely led to the difficulties in classifying the spectral type; the presence of strong, sharp emission lines, and H |$\alpha$| excess likely complicated it as well. Early classifications of HD 160 259 included, A4 se|$\alpha$| (Merrill & Burwell 1933), A2-3 Ia (Wallerstein 1970), and A9 Ia (Houk 1982).

Polarization measurements of HD 160 529 reach back as far as the early 1950’s (Hiltner 1951; Markowitz 1951). No polarimetric variability has been ascribed to the star (e.g. Treanor 1963, Hsu & Breger 1982, Dolan & Tapia 1986), but Dolan & Tapia (1986) do note a complex |$\Delta \theta /\Delta \lambda$|⁠.

2.11 HD 161056

HD 161056 (HR 6601) is a B1/2V star (O’Donnell 1994). Telting et al. (2006) included it in a study looking for line profile variations associated with pulsation; none were indicated, albeit from a single observation. HD 161056 was first observed polarimetrically by van Panhuys Smith (1956) as part of her survey of interstellar polarization in the Southern Milky Way and it is often included in polarimetric studies of the interstellar medium (e.g. Piccone & Kobulnicky 2022). Bastien et al. (1988) only suspected variability in p but reported |$\theta$| variability of 0.5^{|$^\circ$|}; however, later reanalysis calls into question this conclusion (Clarke & Naghizadeh-Khouei 1994). Berdyugin, Snare & Teerikorpi (1995) constrain any variability to |$\lt $| 1^{|$^\circ$|}.

2.12 HD 161471

HD 161 471 (⁠|$\iota ^1$| Sco) is a luminous red supergiant star of spectral type F2Ia (Houk 1978; Gray & Garrison 1989). It is a spectroscopic binary (Pourbaix et al. 2004) and has a 13th mag companion at 37 arcsec separation. It’s H |$\alpha$| line width probably indicates a weak stellar wind (Danks & Dennefeld 1994). It is not a widely used position angle standard, but has been so utilized by Bailey et al. (2023) to calibrate the position angle of polarization in 20-cm PICSARR observations. They find potential variability in p of 0.020 |$\pm$| 0.004 per cent from 18 observations spanning more than a year.

2.13 HD 183143

HD 183143 (HT Sge) is an extremely luminous hypergiant star of spectral type B7Iae (Chentsov 2004). It was first found to have a high broadband polarization by Hall & Mikesell (1950), who described it as a ‘star of special interest’. Serkowski (1974a) later named it as a standard. Clemens & Tapia (1990) found their observations of it to be consistent with those of Serkowski et al. (1975) and Schulz & Lenzen (1983) claim its polarization as a function of wavelength is consistent with an interstellar origin. Spectropolarimetric data from Lupie & Nordsieck (1987) marks it as their least variable standard, but shows |$\sigma _\theta =$| 0.40^{|$^\circ$|}. However, Hsu & Breger (1982) convincingly showed HD 183143 exhibited polarimetric variability (⁠|$\Delta \theta \sim 1^\circ$|⁠, |$\Delta p = 0.2$| per cent) on a time-scale of days to weeks, and Dolan & Tapia (1986) found that its |$\Delta \theta /\Delta \lambda$| character varied from night-to-night, along with |$\theta$| itself.

2.14 HD 187929

HD 187929 (⁠|$\eta$| Aql) is a classical Cepheid with spectral type F6 Ib-G4 Ib and a pulsational period of 7.18 d (Benedict et al. 2022). As the first Cepheid discovered (Pigott 1785), it has been well-studied over the years, and perhaps most particularly during the era of space observations (Evans 1991; Benedict et al. 2007; Evans et al. 2013). The current understanding is that it is a triple system containing, in addition to the Cepheid, a late-B close-in companion as well as an F-type companion lying some 0.66 arcsec from the primary. Polarimetric measurements of the star have focused on attempts to detect a magnetic field using spectropolarimetry to varying degrees of success (Borra, Fletcher & Poeckert 1981; Plachinda 2000; Wade et al. 2002; Grunhut et al. 2010). In linear polarization, Hsu & Breger (1982) report a particularly large |$\Delta \theta /\Delta \lambda$| of |$-7.3 \pm .3 ^\circ /\mu m$|⁠, but no variability. Dolan & Tapia (1986) report measurements that differ by |$1.5^\circ$| from Hsu & Breger’s, along with a complex |$\Delta \theta /\Delta \lambda$| behaviour not explainable by a two-cloud model, and they name intrinsic polarization as a possibility. Bastien et al. (1988) categorized HD 187929 as a suspected variable, but later retract this assessment (Bastien et al. 2007).

2.15 HD 198478

HD 198478 (55 Cyg) is a blue supergiant star of type B3 Ia and a prominent |$\alpha$| Cyg variable with asymmetric contraction varying over hours to days (Wilson 1953). Periods of variability (in pressure, gravity, and modes) appear to correlate with – and are well-modelled by – mass-loss episodes (Kraus et al. 2015; Yadav & Glatzel 2016). HD 198478 may also experience macroturbulence from convection significant enough to contribute to measurable line broadening beyond that from rotation (Jurkić et al. 2011), which may further influence the consistency of some parameters like surface gravity.

Although it is widely used as a standard polarization star (e.g. Cox et al. 2007), large changes in the polarization of HD 198478 have been observed previously by Hsu & Breger (1982) and Wiktorowicz et al. (2023). In particular, Hsu & Breger (1982) saw changes in |$\theta$| and p of 1^{|$^\circ$|} and |$\lesssim$|0.06 per cent, respectively, within a short run – several days. They associated this variability with emission variability seen to occur on the same short time-scale as reported by Underhill (1960) and Granes & Herman (1972). Dolan & Tapia (1986) also suspected variability, partly on the basis that their |$\theta$| determination differed substantially from earlier literature but also because Treanor (1963) noted p was unusually high for its location. Naghizadeh-Khouei (1991) describes its variability in p and |$\theta$| as ‘very obvious’, in particular reporting |$\Delta \theta \sim$| 4.8|$^\circ$|⁠, but observing that |$\Delta \theta /\Delta \lambda$| is consistent from night-to-night. Naghizadeh-Khouei (1991) was critical of the statistical approach of some of the early polarimetric studies, and advocated use of the CDF to aid in matching the polarimetric mechanism to the observed variability. In this specific case he noted the similarity of the variability of this star to that of other supergiants, ascribing it to mass loss and the presence of a stellar wind.

2.16 HD 203532

HD 203532 (HR 8176) is a B3 IV subgiant in the constellation Octans. It is the southernmost standard in the current study with a declination of |$-82.683^{\circ }$|⁠. With this latitude, it is placed close to the molecular clouds south of the Chamaeleon complex, which are associated with the Galactic plane (Larson et al. 2000). Due to coordinate precession being larger for coordinates close to the celestial poles, the position angle changes more over time than for the lower latitude stars (see Table A1). HD 203532 has no known companions nor is it a known variable star (Samus’ et al. 2017).

The first polarization measurement was made by Mathewson & Ford (1970). Later Serkowski et al. (1975) made four measurements yielding position angles between |$126.2 ^{\circ }$| and |$127.6^{\circ }$| in the V band. It has not been reported to be variable, but measurements made by Bagnulo et al. (2017) and Bailey et al. (2020) with modern equipment disagree in |$\theta$| by 2.6 |$\pm$| 0.9|$^\circ$|⁠.

2.17 HD 210121

HD 210121 (HIP 109265) is a B-type star, sharing a line of sight with a single, high latitude cloud sitting |$\sim 150$| pc from the Galactic plane (Welty & Fowler 1992). The star is of uncertain spectral type, with several incongruent classifications having been assigned at different points in the literature. For example, its spectral type has been listed as B3 V (Welty & Fowler 1992; Larson, Whittet & Hough 1996), B3.5 V (Siebenmorgen et al. 2020; Krełowski et al. 2021), B5-6 V (de Vries & van Dishoeck 1988), B7 II (Valencic, Clayton & Gordon 2004; Fitzpatrick & Massa 2007; Bagnulo et al. 2017; Piccone & Kobulnicky 2022), and B9 V (Voshchinnikov et al. 2012). On the whole, a critical reading of the literature suggests to us that an earlier spectral type is more likely; however, we have opted to use the B7 classification in this work, since this is what was used for determining the Serkowski parameters. With the foreground cloud characterized by a high abundance of small grains, HD 210121 is often cited with reference to its extremely low |$R_V$| and high UV extinction.

Initial polarimetric measurements were made by Larson et al. (1996), which show no significant rotation with wavelength. Bagnulo et al. (2017), with greater sensitivity, identified a gradient of |$\Delta \theta /\Delta \lambda = 0.86 \pm 0.07 ^{\circ }/100$| nm, but no variability has been implied.

3 OBSERVATIONS

The data for this work comes from 88⁵ observing runs (or sub-runs⁶) on six different telescopes using four different HIPPI-class polarimeters of three different designs, spanning 10 yr of operation. It includes every observation we have made of the 17 different standards listed in Table 1 during these runs. For 5 standards our data spans the full (almost) 10 yr, two stars (HD 7927 and HD 198478) were observed from the northern site for only a year, the HD 43384 data span is only slightly longer than a year, every other star has a multiyear data set, most of which are at least 5 yr. Many of the observations were made solely for the purpose of |$\theta$| calibration, and some multiband observations were made to check the modulator efficiency (see Appendix  C). However, from June 2020 a number of observing runs were made specifically for this work.

The six telescopes observed with were the 3.9-m Anglo-Australian Telescope (AAT) at Siding Spring Observatory, the Penrith 60-cm (24-in) telescope at Western Sydney University (WSU), a 14-inch Celestron C14 at UNSW Observatory (UNSW), a Celestron 9|$\frac {1}{4}$|-inch telescope at Pindari Observatory in suburban Sydney (PIN), the 36-in telescope at MIRA’s Oliver Observing Station (OOS), and the 8.1-m Gemini North Telescope (GN).

The four HIPPI-class polarimeters include the original HIPPI (Bailey et al. 2015), Mini-HIPPI (Bailey et al. 2017), and two different HIPPI-2 units (differentiated by the hemispheres they operated in, Bailey et al. 2020; Cotton et al. 2022a). They are dual-beam⁷ aperture photo-polarimeters that share common design elements, namely a ferro-electric crystal modulator operating at 500 Hz and optimised for blue wavelengths; a Wollaston (or Glan–Taylor) prism analyser; and modern, compact photomultiplier tube detectors (PMTs). The PMTs are manufactured by Hamamatsu; we mostly used blue (B) sensitive units (Nakamura et al. 2010), but a few observations were also made with PMTs with a redder response curve (designated R). Most observations were made in the SDSS |$g^\prime$| filter or unfiltered (Clear), but a range of other filters were used as described in Appendix  C.

Each of these instruments measure only a single Stokes parameter at once. To measure the other Stokes parameter, the instrument is rotated through 45^{|$^\circ$|}. With HIPPI this was accomplished with the AAT’s Cassegrain rotator. The other instruments instead used their own instrument rotator. In practice an observation was actually made up of four measurements: at angles of 0, 45, 90, and 135^{|$^\circ$|}, where the perpendicular measurements are combined in a way that minimizes any instrumental contribution to the polarization.

A database containing a summary of the instrument and telescope set-up for each run, and the details of every high polarization standard observation in machine readable format is available at CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/MNRAS or on MIRA’s website at http://www.mira.org/research/polarimetry/PA. The catalogue contains both the |$g^\prime$| and Clear (unfiltered) observations analysed in the main part of the paper, as well as those made in other bands used for calibration in Appendix  C (583 observations in total). Note that only the errors derived from the internal statistics of the observation are reported in the catalogue, for an assessment of accuracy see initially, Appendix C4, then Sections 4.3 and 4.4.

3.1 Reduction

Standard reduction procedures for HIPPI-class polarimeters are described by Bailey et al. (2020). These include a bandpass model that takes account of the spectral type of the target. For this work, we have re-evaluated and updated the modulator calibrations to improve the accuracy of this procedure (see Appendix  C). Here, we improve the bandpass model by including the expected polarization and reddening of the standards (given in Appendix  B), the computer code is adapted from the PICSARR reduction (Bailey et al. 2023).

The telescope polarization (TP) is measured by many observations of low polarization standard stars assumed to be unpolarized. The full list of standards used is: |$\alpha$| Aql, |$\alpha$| Boo, |$\gamma$| Boo, |$\eta$| Boo, |$\tau$| Cet, |$\alpha$| Cen, |$\alpha$| CMa, |$\alpha$| CMi, |$\beta$| Hyi, |$\alpha$| Lac, |$\beta$| Leo, |$\alpha$| Ser, |$\beta$| Tau, and |$\beta$| Vir. Aside from |$\alpha$| CMi, all of the additions from those listed in Bailey et al. (2020) have been restricted to Northern hemisphere use. The standard bandpass model is applied prior to the straight average of the standard measurements in each band being taken, then these values are vector subtracted from all the other observations. Aside from run N2018JUN (see Bailey et al. 2020), the TP never exceeded 170 ppm in the |$g^\prime$| band, and varied by typically 10s of ppm or less between bands. Typical TP levels for the telescopes we observed with can be seen in Bailey et al. (2017, 2020) and Cotton et al. (2022b). The nominal errors of this process are mostly less than 10 ppm with larger uncertainties occurring only for some of the smaller telescope runs. These errors are wrapped into the observational errors on an RMS basis – and for the most part are negligible.

The telescope position angle is usually calibrated as the straight average difference between the expected and observed position angles of standards observed in a |$g^\prime$| or Clear filter. The nominal position angles, |$\theta _{g^\prime }$| of the standards have been redetermined from the literature in Appendix  A. The expected position angles include some corrections to these. The two most important are related to the change in position angle with wavelength, |$\Delta \theta /\Delta \lambda$|⁠, and the precession of the co-ordinate system – these are also described in Appendix  A.

The precession turns out to be particularly important. Stars in the North precess positively in |$\theta$|⁠, and vice versa for those in the South – with typical values for |$\Delta \theta /\Delta t$| being several tenths of a degree per century. The literature values for |$\theta$| were first established for some standards more than 50 years ago. Failing to account for precession can artificially induce a degree or more difference between some pairs of standards.

Here, we also apply a correction for the Faraday rotation of polarization (Faraday 1846) by the atmosphere in the presence of the Earth’s magnetic field. This step is recommended by Clarke (2010) when |$p/e_p\gtrsim 1000$|⁠, which many of our measurements nominally surpass. Though never observed, the possible impacts of Faraday rotation on ground-based astronomical polarimetry have been discussed for nearly a century (Lyot 1929; Serkowski 1974b; Hsu & Breger 1982; Clarke 2010). Here, we make the same simplifying assumptions as Hsu & Breger (1982), namely that the Earth’s magnetic field is aligned to geographic poles, and employ,

where |$B_{||}$| is the component of the magnetic field (⁠|$\approx$|0.5 Gauss) parallel to the line of sight, h the scale height of the atmosphere equal to 80 000 cm, X the airmass, and |$V(\lambda)$| the Verdet constant, which is dependant on wavelength. West & Carpenter (1963) give V as 6.27|$\times 10^{-6}$| arcmin Gauss⁻¹ cm⁻¹ for the air at standard temperature and pressure at 567 nm, and we derive the values at other wavelengths by scaling and extrapolating from fig. 2 in Finkel (1964) – at 470 nm this gives 9.53|$\times 10^{-6}$| arcmin Gauss⁻¹ cm⁻¹. At the airmasses of our observations, the correction never comes to much more than 0.01^{|$^\circ$|}.

4 ANALYSIS

4.1 Position angle precision by instrument

Clarke & Naghizadeh-Khouei (1994)’s primary criticism of Bastien et al. (1988) was that measurements from different set-ups were combined in an unweighted way. To conduct a long-term analysis, we need to combine data from multiple (albeit similar) instruments across many sub-runs, where each corresponds to a new mounting of the instrument on the telescope. This task requires some care. So, before we attempt it, we first seek an understanding of the |$\theta$| variability attributable to different telescope/instrumental set-ups.

In Table 2, the standard deviation of |$\theta$|⁠, |$\sigma _\theta$|⁠, of repeat observations in the same filter (limited to |$g^\prime$| or Clear) of the same star within a sub-run is determined; such observations are directly comparable with each other. The table reports the median standard deviation, |$\eta (\sigma _\theta)$|⁠, of each set of observations, where the number of observations, |$N_o\ge 2$|⁠; and the median of |$\epsilon _\theta$|⁠, a metric designed to determine the scatter independent of known noise sources

where |$e_{\rm m\, \theta }$| is the measured uncertainty derived from the internal statistics of our measurements, it is largely photon-shot noise but also incorporates other noise sources associated with seeing and the detector.⁸

It should be noted that this approach is strictly only valid for a Gaussian distribution, and |$\theta$| is not Gaussian. However, it approaches Gaussian at high signal-to-noise, i.e |$p/e_{p} \gtrapprox 5$|⁠, and in our case |$p/e_{\rm m\, p} \gt \gt 100$|⁠, typically.

The different instruments had different rotation mechanisms, which is likely to contribute to scatter in |$\theta$|⁠. HIPPI used the Cassegrain rotator of the AAT to move position angle, while HIPPI-2 and Mini-HIPPI have their own instrument rotator. HIPPI-2’s rotator is a heavy-duty precision component, and on the AAT the value for |$\epsilon _\theta$| is small. Mini-HIPPI’s rotator is not as well made, so the lower precision for Mini-HIPPI was expected. Additionally, the Pindari telescope, unlike the others, is not on a fixed mount so is susceptible to external forces, e.g. wind, incidental operator intervention.

The smaller telescopes do not have as good a value for |$\epsilon _\theta$| as the AAT. Potential mechanical explanations would be a slight polar misalignment or play in the mounting. However, the number of observation sets is small, and so the difference could just be a result of which stars were used in the measurements. For the MIRA-36 this is likely to be a factor since, proportionally more standards were observed that are robustly reported as variable in the literature (and we have excluded HD 198478 from this part of the analysis because it especially biases the results). More reliable standards make up a larger majority of the sets for the other telescope/instrument combinations. To get a more robust measure for the two smaller HIPPI-2 compatible telescopes, we combine their data below the midline of Table 2.

We do not have repeat observations of any standards during the Gemini North run (N2018JUN). Hence, to provide a representative figure we make observation sets for four stars by combining |$g^\prime$|⁠, 500SP and one Clear observation, wavelength corrected using our bandpass model according to the relations in Appendix  A. This is less than ideal because the observations were taken sequentially, but they do at least probe a small range of paralactic angles for the Alt-Az mount. This is important because the telescope polarization during the run was large – probably owing to an inhomogenously aluminized secondary (Wiktorowicz et al. 2023) – and the position angle of its wavelength dependence was not well determined (Bailey et al. 2020).

4.2 Relative position angles using the co-ordinate difference matrix

Informed by the relative performance of the instruments, we now seek to re-establish |$\theta$| values for our standards using our data to facilitate literature comparisons. We use a co-ordinate difference matrix (CDM; Baechler et al. 2020) approach to combine |$\theta$| data from sub-runs where at least two standards were observed in |$g^\prime$| (with the B PMT). Using a single filter reduces the data available, but ensures uncertainties in wavelength effects are minimized. For reasons explained in Section 4.3, we have also removed the Pindari observations from this calculation.

A relatively new tool, CDMs are similar to a Euclidian Distance Matrices (EDM; Dattorro 2005) and are finding applications in a number of different fields (e.g. Mozaffari, Akbarzadeh & Vogel 2019; Krekovic 2020; Liu et al. 2023; Chen et al. 2024). The CDM algorithm optimally calculates the relative mean differences between objects by combining such information from pairs of points. The algorithm is applicable even when the matrix is incomplete, i.e. when measurements do not exist for every object pair. Our application requires a weighted 1D CDM. From Baechler et al. (2020), algebraically, the CDM, |$\boldsymbol{C}$| is made up of elements |$C_{jk} = x_j-x_k$|⁠, i.e. the differences between points |$x_j$| and |$x_k$|⁠. We have a noisy CDM

where |$\boldsymbol{Z}$| is a noise matrix and |$\boldsymbol{W}$| the weight matrix. To optimally recover the set of points |$\lbrace x_{j}\rbrace _{j=1}^{N}$| that generated it, the solution is

where |$\boldsymbol{A} = \boldsymbol{\Lambda }-\boldsymbol{W}$| and |$\widetilde {\boldsymbol{\nu }} = (\boldsymbol{C}\circ \boldsymbol{W})\boldsymbol{\mathbb {1}}$|⁠, where |$\boldsymbol{\mathbb {1}}$| is the all-one vector and

The first point is then fixed to zero in the algorithm by removing the first point in |$\boldsymbol{x}$| and |$\widetilde {\boldsymbol{\nu }}$| along with the first row and column of |$\boldsymbol{A}$| and likewise for all the corresponding matrices. From this process, we take the nominal errors on the recovered |$x_j$| points to be

Here, our x terms are |$\theta$| values, and we calculate the weights using the root-mean-square (RMS) sum of |$e_{\rm m\, \theta }$| and |$\eta (\epsilon _{\rm i\, \theta })$| for each measurement. Where multiple observations of a standard have been made during a sub-run, we take the error-weighted-average, and thus also the resulting RMS-error for the weight calculation. The runs are combined in the same way. The result therefore makes no account of stellar variability – only mean position angles are calculated.

Fig. 2(a) shows the number of object pairs used in the calculation and 2(b) the associated standard errors (where the weights are the inverse of the error squared). It should be noted that this does not match the number of observations because only a single agglomerated measurement is made per sub-run for a given standard, single-standard sub-runs are discarded, and because more standards were observed in some runs than others.

The results of the CDM procedure are presented in columns 6–10 of Table 3 where the standard error, |$\theta _e$|⁠, given in column 7 is the error in the mean |$\theta$| given in column 6; this does not account for the error in the zero-point of the co-ordinate system or does it describe the distribution of |$\theta$| values. Columns 9 and 10 are the weighted standard deviation of |$\theta _{\rm obs}-\theta _{\rm pred}$|⁠, |$\sigma _\theta$|⁠, and the average error, |$\bar{e}_\theta$| (which is equal to |$\theta _e\sqrt{N_o}$|⁠) after calculation of |$\theta _t$| using the newly determined values of |$\theta$| for each standard. Column 9 may be compared to column 5, which is the unweighted standard deviation after calibration of |$\theta _t$| instead using |$\theta _{\rm lit}$| – the way calibration of |$\theta _t$| has previously been done for HIPPI-class instruments. It can be seen that the CDM-derived |$\sigma _\theta$| is reduced for most stars compared to the previous method. Our |$\sigma _\theta$| and |$\bar{e}_\theta$| terms are the equivalent of Bastien et al. (1988)’s |$\sigma _2(\theta)$| and |$\sigma _1(\theta)$|⁠, respectively.

The CDM only calculates the relative position angles of the standards (the value of the first listed object in the matrix is arbitrarily set to zero, so that |$\theta _{\rm diff}$| represents the difference from it for each standard). For the absolute value, we have calculated and applied an offset, |$\zeta$|⁠, based on the literature values (as given in Table A1). Here, |$\zeta$| is calculated by finding the median difference between the determined, |$\theta _{\rm diff}$|⁠, and the literature, |$\theta _{\rm lit}$| values

We estimate the error in |$\theta _0$| by this method to be 0.177^{|$^\circ$|}. In this case, we calculated errors as the RMS sum of the values in Table 3 and 0.5^{|$^\circ$|}, which we consider an appropriate typical uncertainty in the original literature measurements.⁹ HD 203532 is removed from the calculation as an obvious outlier.

An alternative method for determining |$\zeta$| would be to use the error-weighted average difference

where the |$w_j$| is the weight of the jth element, equal to the inverse of the error squared. Hsu & Breger (1982) used this method to tie their measurements to those of Serkowski (1974a).

Calculating |$\zeta$| using equation (13) rather than equation (12) results in all the values in Table 3 being shifted by |$-0.210^\circ$|⁠, which is not so different from the zero-point uncertainty calculated above. We prefer the median approach to mitigate sparse sampling of potentially variable objects. These values describe the possible offset of our entire network of standards but bare no relevance for their relative position angles.

4.3 Estimating stellar variability in |$\theta$|

In Section 4.1, we estimated the uncertainty associated with each instrument/telescope combination, we now wish to do the same for each standard star, |$e_{\star \theta }$|⁠. A first step is to determine the (error weighted) standard deviation of position angles from their expected values for each standard – after re-calibrating |$\theta _t$| for each run – these are the values in column 9 of Table 3. The general formula for weighted standard deviation is (Heckert & Filliben 2003),

where |$\bar{x}$| is the weighted mean, |$w_j$| the weights on each element |$x_j$|⁠, |$N^\prime$| is the number of non-zero weights¹⁰, and the term |$(N^\prime -1)/N^\prime$| comes about from Bessel’s correction to the variance, applied when the population mean is being estimated from the sample mean.

However, estimating |$e_{\star \theta }$| in this way is very crude because |$\theta _t$| depends on the weighted means of each observation for each run. Typically, each run contains only a handful of measurements where we have considered each standard to be as good as any other. If just one star is variable and significantly different to its assumed value of |$\theta$|⁠, that will shift the measurements of every other star in the run. The very parameters we aim to determine are corrupted by assuming them to be zero in the first instance.

To overcome this problem, we employ a scheme that iteratively fits for |$e_{\star \theta }$|⁠: For this we use SciPy’s optimize function, employing the Sequential Least SQuares Programming Optimizer (SLSQP), to minimize a function f,

where |$\chi _r^2$| is given by

where d is the degrees of freedom, equal to |$N_o$| less the number of fit parameters, |$(O_j-E_j)$| is the difference between the observed and expected |$\theta$| value and |$\sigma _j$| the uncertainty associated with each measurement after re-calibration, which has four contributions,

|$e_{\rm m\, \theta }$| is the measured uncertainty, as in equation (7); |$e_{i\, \theta }$| and |$e_{\star \theta }$| are the standard deviations of the variability (errors) associated with the telescope/instrument set-up and the star being observed, respectively, and |$e_{t\, \theta }$| is the RMS error in the determined telescope zero-point for the run (which, for each run, is derived from the other three terms).

However, |$(O_j-E_j)$| depends on the values we assign |$\theta$| for each star. The assignments are made using the CDM algorithm, which in turn depends upon the uncertainty assigned to each observation. Our assignments in Section 4.2 are made, essentially, with |$e_{\star \theta }=0$|⁠. Therefore our fit function incorporates a recalculation of the CDM-derived assignments based on the current values of the fit parameters.

We found that the fit function was prone to getting stuck in local minima, leading to some values of |$e_{\star \theta }$| being inconsistent with the subsequently calculated values of |$\sigma _\theta$|⁠. This occurs because the CDM calculates only one |$\theta$| value per star per run, and the fit function just demands a model where the sample variance is appropriately described, it does not care where the uncertainty is placed. To overcome this problem, we used the |$\sigma _\theta$| values determined after the fit as initial parameters for |$e_{\star \theta }$| in subsequent fits and iterated until the |$\sigma _\theta$| values converged (i.e. did not change by more than 5|$\times 10^{-5\, \circ }$| between subsequent iterations) though the changes are small after about the third iteration.

Initially we performed the analysis of Section 4.2 using the same set of observations as in Section 4.1, but found that the three standards observed repeatedly at Pindari gave elevated values of |$e_{\star \theta }$| compared to otherwise. We presume this to be due to the complications associated with not using a fixed mount, and so we excluded those observations from further analysis. We also removed the Clear observations; stellar variability is likely to have a wavelength dependence and we wanted an unbiased comparison between stars. We could not apply the same restriction in Section 4.1 since this would have left too few sets for some set-ups.

It follows that the uncertainty for each set-up could use some fine tuning. We tried fitting |$e_{\rm i\, \theta }$| along with |$e_{\star \theta }$| but this just transferred uncertainty from the set-ups to the stars, where the proportion was overly sensitive to the initial conditions. Though, if boundary conditions were enforced, the uncertainty transferred was mostly small. However, the ratios of |$e_{\rm i\, \theta }$| to each other were typically only slightly changed, suggesting that the results of Section 4.1 are not biasing observations made with one set-up over another. It also means that if most stars are variable, our estimates of |$e_{\star \theta }$| will underestimate the short term variability of the most stable stars. Ultimately we decided against fitting |$e_{\rm i\, \theta }$|⁠, as this represented the most conservative approach, and retain the values assigned in Section 4.1.

The results are given in columns 11–18 of Table 3. There are only minor changes in |$\theta$| (column 11) compared to the pure CDM results (column 6). There are marked decreases in |$\sigma _\theta$| for the least variable stars, that are only partly offset by increases in the more variable stars. This result should be obvious because we are down-weighting the contributions of the most variable stars to achieve a more accurate determination of |$\theta _t$|⁠. The errors are increased because they now include a component associated with |$e_{\star \theta }$|⁠.

The significance of the determined |$\sigma _\theta$| we calculate (in column 16 of Table 3) by scaling each observation to the corresponding error, and then taking the error-weighted standard deviation. This is better than dividing |$\sigma _\theta$| by |$\bar{e}_\theta$| because it takes proper account of the different uncertainties associated with each observation. In column 17 the fit values of |$e_{\star \theta }$| are given. Due to the nature of the fitting function |$e_{\star \theta }^2 \approx \sigma _\theta ^2-\bar{e}_\theta ^{\, 2}$|⁠. Five stars are found to be variable with 3|$\sigma$| confidence – all of these are A or B type supergiants. Another four stars are seemingly variable only at the 2|$\sigma$| level.

The final column of Table 3 is a calculation of polarization, following equation (2), that would be required in |$u^\prime$| to rotate by |$e_{\star \theta }$| assuming |$q^\prime = p_{\rm lit}$|⁠. This allows the variability of the objects to be compared directly, independent of the interstellar polarization imparted by the ISM.

The |$u^\prime$| and |$q^\prime$| values are depicted graphically as the vertical and horizontal axes, respectively, of QU diagrams in Fig. 3 (full explanatory details are in the caption). The error bars plotted for |$q^\prime$| are derived below in Section 4.4, but we present this figure here to make an important point: the improvements made in |$\theta$| calibration have resulted in this being more accurate than our p calibration. We have used three different modulators, where the performance of at least one of them has evolved over time, this affects modulation efficiency, including as a function of wavelength. Our ability to correct for efficiency changes is limited by the calibration data obtained (see Appendix  C). But any change in modulator efficiency applies equally to each Stokes parameter so that the effects on |$\theta$| are largely negated – especially in the bands most closely corresponding to the optimum operating wavelength of the modulator (i.e. the |$g^\prime$| band, see Table C1) – the instrumental errors in |$\theta$| come from other sources (largely mechanical), as already discussed in Section 4.1, and these are independently parameterised and accounted for. The consequences are easy to see in Fig. 3 – where the co-ordinate system is rotated, for each object, so that |$q^\prime$| corresponds to changes in p and |$u^\prime$| to changes in |$\theta$| – where the horizontal scatter is almost always greater than the vertical scatter. This means we cannot employ a CDF test reliably. However, where we see variability in |$\theta$| we can be confident it is real, but if there is only scatter in p this probably represents only instrumental variability.

4.4 Variability in p

The p statistics for all of the available |$g^\prime$| observations are given in Table 4, except for those acquired at the Pindari observatory, which are given in Table 5. The Pindari observations have larger nominal (photometric) errors, but otherwise should be comparable.

Some standards do not agree as well with the predictions as others. The mean disagreement is quantified in Table 4, in terms of the unweighted mean difference, |$\bar{\Delta p}$|⁠. HD 210121 is almost 6 per cent under prediction as a ratio. Other stars that differ by more than two per cent are HD 23512 and HD 183143 which are also underpolarized; and HD 203532 and HD 154445, which are overpolarized.

Conservatively, a long-term change in p from the literature value would be indicated for any star with |$\bar{\Delta p} \gt 3\sigma _{\rm \Delta p}$|⁠. None of the standards listed in Table 4 meet this condition. The nearest is HD 210121 which is significant only at |$2\sigma$|⁠. Some correction is, however, probably still justified, but because the measurements are monochromatic it is not possible to say if it is |$p_{\rm max}$| or |$\lambda _{\rm max}$| that is different. Indeed, we have calibrated the modulator polarization efficiency to all the standard observations collectively, so all that one can say is that several stars deviate in p compared to others based on the source literature, we can’t say which are inaccurate.

Table 4 also reports the unweighted standard deviation in |$\Delta p$| compared to the mean, and the mean measured error, |$\bar{e}_{\rm m\, {\it p}}$|⁠. If there are no further contributions to the error in polarization, |$e_p$|⁠, then the scatter in p attributed to unaccounted for noise sources, |$\epsilon _{\Delta p}=\sqrt{\sigma _{\Delta p}^2-\bar{e}_{\rm m\, {\it p}}^2}$| represents the stellar variability in p, i.e. |$e_{\star \, p}$| or equivalently |$e_{\star \, q^\prime }$| if we rotate the reference frame in the same way as first done in Section 4.2 (i.e. this is the counter-part to column 18 in Table 3). In Fig. 4, we have plotted |$\epsilon _{\Delta p}$| against |$e_{\star \, u^\prime }$|⁠. Owing to the large interstellar polarizations of these stars, whether |$e_{\star \, u^\prime }$| or |$e_{\star \, q^\prime }$| is larger, should be purely a matter of chance. We would therefore expect an equal number of points to fall either side of the blue line if |$p_e$| is dominated by |$e_{\rm m\, {\it p}}$|⁠. This is clearly not the case, so there is an additional instrumental error, |$e_{i\, p}$|⁠, that contributes to |$\epsilon _{\Delta p}$| that needs to quantified, i.e.

Many more points are above the blue line in Fig. 4 at low values of |$e_{\star \, u^\prime }$|⁠, so a fixed error is more appropriate than one that scales with p. In re-calibrating the modulators (Appendix C4), we found a typical disagreement between |$p_{\rm obs}$| and |$p_{\rm pred}$| of 460 ppm. Using this figure overestimates the error (magenta line in Fig. 4) because it incorporates inaccuracies in the literature values of p into the metric. If instead we take the median difference to the mean value of |$\Delta p$| for each star, the result is |$\approx$|290 ppm for the observations in Table 4 and |$\approx$|155 ppm for the Pindari observations. The fact the precision is better in p for the Pindari observations, which used only a small telescope but a single instrument/telescope set-up, is evidence that variation in instrument and set-up is the dominant source of the extra error. The purple line in Fig. 4 represents 280 ppm of added error; this appears to be an underestimate – too many points fall above the line, probably as a consequence of taking the mean from small data sets with uneven temporal sampling. The red line in the same figure is the intermediate value of 375 ppm, which more evenly divides the points, and which we adopt as |$e_{i\, p}$|¹¹

Having established a value for |$e_{i\, p}$|⁠, we can now calculate error-weighted figures for |$\sigma _p$| and |$\bar{e}_p$| to determine the significance of the result, this is given in the final column of Table 4. Two stars are found to be variable at 3|$\sigma$| significance in p; both are also 3|$\sigma$| significant variables in |$\theta$|⁠. Another four stars are 2|$\sigma$| significant; two of these are not similarly significant in |$\theta$|⁠: HD 154445 and HD 161056; they both have early B spectral types. HD 111613 is a 3|$\sigma$| significant variable in |$\theta$|⁠, but does not approach this level of significance in p.

4.5 Variability time series

Here, we present the |$\theta$| and p data for each star as a time series in Figs 5 and 6, respectively. The latter includes data from the Pindari observatory, whereas the former does not, but otherwise excludes any point not on the first plot. Both figures include not just the |$g^\prime$| data from runs with multiple stars observed, as used in Sections 4.2 and 4.3, but also Clear data (open symbols), and |$g^\prime$| and Clear data which were re-calibrated based on Clear observations (where there are insufficient |$g^\prime$| standard observation for the run). The additional data may be less accurate but is useful to fill in gaps in the time series (e.g. the 2016 HD 160529 datum). We have colour coded the observations by run as described in the caption. This is important because despite our precautions it is still a more precise matter to compare observations made within the same run. The error bars incorporate the error in |$\theta _t$|⁠, but if comparing observations intrarun the uncertainty is less than this.

Long-term trends, or regular periodicity – where sufficiently sampled – would be revealed in Figs 5 and 6, but neither behaviour is obvious. Considering both Figs 5 and 6, the time-scales of variability in |$\theta$| and p mostly appear correlated. For instance, the fast variability of HD 198478 is very easy to see within a run. Whereas the slower, but no less pronounced, variability seen in HD 160529 or HD 80558 occurs on longer time-scales.

A noteworthy discrepancy is seen in Fig. 6 when comparing data from two early 2018 runs for HD 154445 – this is the only abrupt change in polarization observed for this star and is most likely not real. The Gemini North run (2018JUN) has the least reliable calibration, and this represents the first of these data points. Closer inspection reveals that all four observations from this run appear over-polarized. If we remove these points from the analysis in Section 4.4 (see below the mid-rule in Table 4) then the significance of variability in p for both HD 154445 and HD 161056 falls below 2|$\sigma$|⁠, and so we regard neither as a variability candidate.

5 DISCUSSION

5.1 Changes in |$\theta$| over decades?

Four stars have refined |$\theta$| values (column 11 in Table 3) significantly different to the literature using the criteria |$|\Delta \theta | \gt 3\sigma _{\theta }$|⁠: HD 149757, HD 161056, HD 203532, and HD 210121 (other differences might be ascribed to stellar variability). This could indicate a slow change in |$\theta$| over many decades. However, inspection of Fig. 5 does not support this notion. The data for these four stars do not hint at a long term trend despite spanning 4 years or more. If there is one it must only be apparent on longer time scales.

We sourced |$\theta _{\rm lit}$| from only Bagnulo et al. (2017) for two stars: HD 161056 and HD 210121. HD 203532, includes data from Bagnulo et al. (2017) and Serkowski et al. (1975). Together, these are three of the five stars where |$\theta$| was sourced from Bagnulo et al. (2017) (another star sourced from Bagnulo et al. (2017) is HD 80558, it has a large |$\Delta \theta$| value, but is apparently also more variable). This is curious, because Bagnulo et al. (2017) is our most recent reference source, leading us to discount a slow drift in |$\theta$| over time. Of the stars mentioned, three were observed by Bagnulo et al. in (2015) and have large negative |$\Delta \theta$| values, HD 161056 was observed later in 2017 and has a large positive |$\Delta \theta$| value. We conclude the difference is probably due to inaccurate |$\theta _t$| calibration with FORS2 on the VLT.

This leaves HD 149757 as the only star of interest; for it we sourced |$\theta _{\rm lit}$| from Serkowski et al. (1975), modified for |$\Delta \theta /\Delta \lambda$| with data from Wolff et al. (1996) to give 127.2^{|$^\circ$|}. Our determination is almost 1^{|$^\circ$|} below this. The spectropolarimetry of Wolff et al. (1996) produces a figure of |$\approx 125^\circ$| for the |$g^\prime$| band, while that of Wolstencroft & Smith (1984) gives |$\approx 126^\circ$|⁠. B band observations described by McDavid (2000) range from |$129^\circ \pm 0.4$| to |$126.2^\circ \pm 0.2$|⁠; this seems like a lot of variation but it is not too different to that reported for other standards in McDavid (2000)’s agglomerated tables. Fig. 5 shows that few, if any, of these stars are really variable in |$\theta$| over such a range.

More tellingly, perhaps, is the difference in |$\Delta \theta /\Delta \lambda$| behaviour between Wolff et al. (1996) and Wolstencroft & Smith (1984) – the slope is completely opposite in the two cases. This may point to intrinsic polarization that is episodic in nature. This would be in keeping with the irregular variability of the star’s photometry and spectroscopy (see Section 2.8). If so, it is not captured by the |$g^\prime$| observations we analysed, we can assign only 0.211^{|$^\circ$|} to stellar variability. However, there is some evidence for a more active era early in 2017 in Clear observations made with Mini-HIPPI at the UNSW observatory – these can be seen to range over |$\sim 2^\circ$| in Fig. 5.

5.2 Variability in |$\theta$| compared to the literature

5.2.1 Literature data

In Fig. 7, we compare our determinations of intrinsic stellar variability to those available in the literature. As noted by Naghizadeh-Khouei (1991) only partial data is presented in some of these sources, and others make only a handful of measurements. Since there are so few determinations available, we employ very relaxed criteria to include as many comparisons as possible.

The core of the literature determinations presented come from Bastien et al. (1988) who derived variability in |$\theta$| from the standard deviation of q and u (⁠|$\sigma _2(\theta)$| comparable to our |$\sigma _{\theta }$|⁠) and compared that to the average error, |$\sigma _1(\theta)$| (equivalent to our |$\bar{e}_\theta$|⁠). Their data includes typically tens of observations per star agglomerated from half a dozen different sources that made use of one of two of the better instruments available between 1983 and 1986.

Hsu & Breger (1982) present only representative data. They describe the data in detail, but do so in inconsistent ways that are not conducive to analysis. They actually observed each star tens of times between 1979 and 1981 (Breger & Hsu 1982), but this is lost to us. All we have been able to do is digitise the data they plotted in their Figs 1 and 2. We take their typical reported error of 0.1^{|$^\circ$|} as representative of |$\bar{e}_\theta$| and calculate the standard deviation of the presented data to derive |$\sigma _{\theta }$|⁠. The data in their Fig. 1 represents a single 6 night observing run from August of 1981, Fig. 2 adds observations from 4 nights in September of the same year. We neglect the R band data available for one target and use only the V band observations.

Dolan & Tapia (1986) observed three of the same stars we did. They were mostly concerned with |$\Delta \theta /\Delta \lambda$| but in their Table 1 they present an average |$\theta$| from all bands for 5 observations of HD 7927 from 1984, 3 of HD 111613 from 1980 and 1984, and 8 of HD 183143 from 1980 and 1984. We calculate standard deviations from this data. They report that their standard error is 0.09^{|$^\circ$|} based on a test from injecting a polarized signal into their system.

From Wiktorowicz et al. (2023) we take the B band data for 3 stars, including HD 154445 which has only 3 observations; the others between 6 and 13. All the data sets span many years. Two data sets, from two different telescopes, are presented for HD 149757, we have plotted both. They give what is effectively |$\bar{e}_\theta$| directly in their Table 17, and we calculate |$\sigma _{\theta }$| from their |$q_{\rm var}$| and |$u_{\rm var}$| values in their Table 18 using Bastien et al. (1988)’s formula.

The variability of HD 43384 has been much discussed (see Section 2.3), we have two sources of data for HD 43384. Coyne (1972) presented data from 1969 to 1971 in two bandpasses, we take a high |$S/N$| (⁠|$e_p\le 0.2~{\%}$|⁠) subset of the more extensive dataset (G band) which then totals 9 observations; |$\theta$| is only reported to the nearest degree, and so we add an additional 0.5^{|$^\circ$|} RMS error. Matsumura et al. (1998) made 17 observations relative to HD 21291¹² with a single instrument with a typical error of 0.4|$^\circ$| between 1991 and 1996. In the paper these data are plotted against phase assuming a 13.70-d period, with a sine curve fit. We digitised the points and calculated the standard deviation.

Finally, the values for one point for HD 149757 comes from McDavid (2000), who reports an insignificant variability based on agglomerated data from 1949 to 1996; the figures we have used come from the ‘Grand Average’ of multiband data in their Table 11 (though all the bands are pretty similar).

5.2.2 Comparison

The picture that emerges from Fig. 7 is one of surprisingly little disagreement. In particular Bastien et al. (1988)’s determinations are quite a good match for ours; with the only notably discrepant point being HD 7927. There are actually four literature determinations of HD 7927’s variability, and they disagree wildly. Our data on this target spans only a year. So, if the different estimates are to be reconciled it can only be through an appeal to episodic behaviour.

More broadly, it can be seen that variability determinations made from very short observing runs (open symbols in Fig. 7) yield low values (the exception is HD 198478 which exhibits obvious variability within a few nights). Amongst such determinations, we find good agreement with HD 147084 and HD 154445 which are two of the least variable stars in |$\theta$|⁠. In every other case there is disagreement with our data that spans years.

We do find agreement, largely within 1|$\sigma$|⁠, with almost everything that has been observed over a longer period. The biggest disagreements, still within 2|$\sigma$|⁠, are for HD 43384 which has some of the largest literature uncertainties, and HD 183143 for which we have only 7 observations over a year’s duration.

Despite the misgivings expressed by Clarke & Naghizadeh-Khouei (1994) and Naghizadeh-Khouei (1991), our analysis is consistent with the findings of earlier authors that many standards are in fact variable. A casual perusal of stellar polarimetry papers from the last 40 years will reveal that observers have largely ignored these reports of variability. This is partly because bright high polarization stars are not that common and alternative calibration methods not readily available. Therefore, the point is not to be able to declare a star variable or not based on statistical criteria, since this will not help observers much. There is enough evidence to suggest variability is common that we should consider it so. Yet, without quantification, this is an equally unhelpful statement. Any departure from spherical symmetry will produce a polarized signal, so it is a truism to write that everything is a polarimetric variable. The important question then is not which standards are variable? but rather what uncertainty should we ascribe to our measurements when using these standards? Our contention is that a careful analysis will assume the best estimate available, even if it is smaller than the formal detection threshold.

5.3 Astrophysical inferences

The data presented in this paper represents a rare opportunity to study the long-term small-scale polarimetric variability of high polarization stars with precision. In Table 6 – broken down as early-type supergiants, late-type supergiants, and other stars; and ordered by spectral type – the values of |$e_{\star u^\prime }$| and |$e_{\star q^\prime }$| are tabulated; the final column combines the two figures in a variability metric,

The groupings reveal B/A-type supergiants to be more variable as a category than either the F/G-type supergiants or the stars of other classes. In Figs 5 and 6 time series data is plotted, by star, in |$\theta$| and p, which we use along with the table and the QU diagrams in Fig. 3 to elucidate the nature of the variability below.

A thorough review of the QU patterns associated with different polarigenic mechanisms is given by Clarke (2010). In brief, periodic mechanisms, relating to binarity or persistent surface features – like magnetic spots – will be revealed by loops or figure-of-eight-type patterns. Polarization associated with a Be mechanism – material ejected from the equator followed by slow decretion from a disc – which has a preferred orientation, will fall predominantly along a straight ‘intrinsic line’. Polarization generated by randomly distributed clumps of gas within a (perhaps irregularly driven) stellar wind will manifest as a scatter plot. Whereas, an interstellar polarization drift would probably look like a random walk, but predominantly in p.

5.3.1 |$\alpha$| Cyg variables

Coyne (1972) was the first to (statistically) study the polarimetric variability of supergiant stars; which he found twice as variable as his control group. Together Fig. 3 and Table 6 are quite revealing in showing that it is the B/A-type supergiants – which are all |$\alpha$| Cyg variables – that are clearly more variable as a class, with all but one being formally 3|$\sigma$| detected polarimetric variables. By and large their QU diagrams show no favoured direction nor clear pattern, something generally associated with a clumpy stellar wind mechanism (see refs. within Bailey et al. 2023).

In the era of space-photometry many early-type supergiants have been classified as |$\alpha$| Cyg variables on evidence of what is sometimes called stochastic low-frequency variability. In the GCVS catalog |$\alpha$| Cyg variables were originally classified as A- or B-type supergiants that displayed semi-regular radial velocity (RV) variability, said to be associated with the beating of many long period non-radial pulsations (Samus’ et al. 2017). Often, it is only assumed that the photometric and RV variability are associated. Nevertheless, many such stars are also known polarimetric variables, including for instance |$\lambda$| Cep (Hayes 1975), |$\chi ^2$| Ori (Vitrichenko & Efimov 1965), Rigel (Hayes 1986), and Deneb itself (Cotton et al. 2024).

HD 80558, HD 111613, HD 183143, HD 183143: All four of these stars have past claims of variability. HD 183143 and HD 198478 were found to be variable in position angle by Hsu & Breger (1982). HD 183143 and HD 198478 are among the most clearly variable stars in our data; they have values of |$e_{\star \theta }$| larger than 0.5|$^\circ$|⁠. Dolan & Tapia (1986) later reported that HD 111613 and HD 183143 were variable in both p and |$\theta$|⁠. HD 183143 shows perhaps the most dramatic shift in polarization we observed, and we also find it variable in p. It is not surprising that it should be a polarimetric variable given the noteworthy emission in its spectral lines. It exhibits P Cygni profiles in H|$_{\alpha -\delta }$| indicative of a strong wind, that if clumped would be an asymmetric scattering medium. Bastien et al. (1988) reaffirmed HD 111613 and added HD 80558 as |$\alpha$| Cyg stars that were polarimetric variables. Clarke & Naghizadeh-Khouei (1994) criticised Bastien et al. (1988) but nevertheless agreed that HD 111613 was a polarimetric variable. We find HD 80558 to be variable at 3|$\sigma$| significance in |$\theta$| – easily seen in Fig. 3 – corresponding to an intrinsic variability of almost 0.4^{|$^\circ$|} – and almost 2|$\sigma$| in p. Thus, it seems the variability of these four stars is persistent.

HD 43384 has broad past evidence for variability (see Section 2.3). In particular, it was found to be variable in |$\theta$| by Hsu & Breger (1982). However, it has the lowest value of |$e_\star$| of any of the six early-type supergiants we observed, and is the only one for which there is no formal variability detection. Matsumura et al. (1998) claim the star’s variability is periodic, but our data represent poor phase coverage for testing the reported period. So this is a potential cause of the disagreement. However, since we have just over a year’s data on this star, it is also possible this is simply insufficient for capturing episodic variability.

HD 160529 is the only B/A-supergiant without a prior variability claim; it is invariable in the short-run data presented by Hsu & Breger (1982). Repeat observations we made of the star in the same run also agree well. However, we see very large changes in |$\theta$| and p over years. Deneb, which is of the same spectral class and type, also exhibits polarimetric variability only on longer time-scales than other |$\alpha$| Cyg variables (Cotton et al. 2024). As an LBV star, HD 160529 is undergoing mass-loss that varies on a range of time-scales associated with winds, pulsation, and rotation (Stahl et al. 2003). As such, it is rather surprising that it has not previously been noted as a polarimetric variable. Another LBV star, P Cygni, has recently been the subject of a decade of spectropolarimetric study by Gootkin et al. (2020), who find that it displays a variety of polarimetric behaviour, in part associated with free-electron scattering off of clumps (mostly) uniformly distributed around the star. We deduce that the large polarimetric variability we see here in HD 160529 – amounting to |$e_\star =$| 0.21 per cent; |$e_{\star \theta } =$| 0.64|$^\circ$| – is associated with similar phenomena. Mass-loss from these stars is greatest during outburst. For HD 160529 this corresponds to an increase in temperature of |$\sim$|4000 K and a size increase of 180 |$R_\odot$|⁠. This process can take a number of years, consistent with the slow variation seen in Figs 5 and  6.

Hayes (1984) made studies of the two early type |$\alpha$| Cyg stars |$\alpha$| Cam (O9.5 Ia) and |$\kappa$| Cas (B1 Ia). These two stars show both slow long term variability and faster short term variability in their polarization, in much the same way as do those early-type supergiants observed here with sufficient regularity. For instance, slow longer term variability is evident in HD 80558, faster short-term variability can be seen in both HD 80558 and HD 198478. We have not sampled all the B-type supergiants sufficiently on both time-scales to test that these behaviours are simultaneously displayed, none the less it is a fair hypothesis. Hsu & Breger (1982) found early-type supergiants more variable on short time-scales than other standards and recommended against their use as such on that basis. Hayes (1984) concluded that the polarization variations he saw were the result of varying mass-loss in the stellar extended atmospheres and that he was witnessing the ‘waxing and waning of non-periodic clumpiness in the envelope’, such as to explain the variability at both time-scales.

That similar polarimetric behaviour might be universal in early-type supergiants is suggested by Hayes (1986) similar findings regarding polarimetric variability in the later type (B8 Iae) Rigel. There, the variability was slower and the conclusion was that it resulted from localized disturbances in the stellar envelope originating at or below the photosphere. Cotton et al. (2024) recently observed analogous behaviour for (A2 Ia) Deneb and came to similar conclusions. The variability of (A2 Ia) HD 111613 seen by Bastien et al. (1988) is clearly slow in the same way, as is that we see in (A2 Ia) HD 160529. Together these data present a picture where the time-scale of polarimetric variability in |$\alpha$| Cyg stars varies with spectral type. Given the contention that non-radial pulsations might be responsible for RV and photometric variability in similar stars (Bowman et al. 2019) and the associated debate about the plausibility of that hypothesis (Lenoir-Craig et al. 2022), we suggest that the time is right to examine the polarimetric variability of early-type supergiants as a class with a view to informing this question.

5.3.2 Suspected variables

Our observations of little intrinsic polarization amongst the F/G-supergiants is in keeping with expectations based on their cooler photospheres, which are less ionized, and produce less prominent winds. The standards of the other classes are of earlier types, but overall their variability is similar. There are no formal 3|$\sigma$| detections of polarization variability amongst any of these other stars. However, the levels of variability we record are generally a match for Bastien et al. (1988)’s findings, and four stars show 2|$\sigma$| detections; these might therefore be regarded as suspected variables. We briefly comment on each of these below:

HD 7927: We have already briefly discussed HD 7927 in Section 5.2. Bastien et al. (1988) found a greater |$\theta$| variability for this star than we record, Hsu & Breger (1982)’s presented data gives a lower value. Bastien et al. (1988) also found the star to be significantly variable in p, however we do not—our average error is about double theirs, though this does not wholly account for the discrepancy. It seems likely the star’s behaviour is episodic. Being of spectral type F0 Ia it does not meet the formal definition for an |$\alpha$| Cyg classification, though it is adjacent. As noted in Section 2.1 it displays the characteristics associated with the class, namely irregular RV and stochastic photometric variability, albeit at a low level. If the star does belong to the same class of objects then its polarimetric variability profile looks more like Rigel’s in Hayes (1986) than the earlier type supergiants studied in Hayes (1978, 1984).

HD 23512: Variability in |$\theta$| for HD 23512 is 2|$\sigma$| significant, corresponding to |$e_{\star \theta }=$| 0.256|$^\circ$| – with which Bastien et al. (1988)’s results are consistent. This is surprising since as an ordinary A0 V star it should be one of the least intrinsically polarized (Cotton et al. 2016), a fact already alluded to by Hall (1951).¹³ This suggests that its companion may be a polarimetric variable. The companion is responsible for 16percnt of the combined flux, but being of small separation and discovered by Lunar occultation¹⁴, little else is known about it. Breger (1984) confirmed Hall (1951)’s initial finding that HD 23512 is much more polarized than any of the other bright stars in the cluster, but showed its wavelength dependence was typical. This leads us to speculate that the companion might instead be a bright background star. As discussed by Guthrie (1987), HD 23512 is unusual in other ways: it displays anomalously weak absorption at 2200 Å and in several other bands for a star of its spectral type, and is depleted in calcium. These factors may provide some clues as to the nature of the other star. If the ‘B’ component is predominantly responsible for the polarization variability, the scale of it would be |$\sim$|2000 ppm, which is not unreasonable for a supergiant or a Be star.

HD 187929: Bastien et al. (1988) describe |$\theta$| variability consistent with our measurements. The star has a magnetic field detection, but it is too small to produce significant linear broadband polarization (Barron et al. 2022). Classically, |$\delta$| Cep variables (Cepheids) should not be polarimetric variables as radial pulsations produce no net polarization change (Odell 1979). However, Polyakova & Sudakov (1981) and Polyakova (1984, 1987, 1990) claim to have measured phase-dependant polarization that is bound by a three-lobed ‘rosette’ in the QU diagrams of many such stars. Polyakova & Sudakov (1981) ascribe this to an artefact of calibration owing to broadband measurements of colour-variable objects, but Polyakova (1990) later suggested an intrinsic polarization mechanism due to a circumstellar envelope 20–30 per cent larger at the equator than the pole – Kervella et al. (2006) thought this worthy of investigation.

We plotted the p and |$\theta$| against phase for both |$\delta$| Cep variables in Fig. 8, since nothing was obvious from the QU diagrams in Fig. 3. Surprisingly, there is a hint of phase dependant behaviour in the binned data of HD 84810 consistent with that described in the literature, albeit well below our formal detection thresholds—investigating this is beyond the scope of the paper. There is no clear phase-dependant behaviour for HD 187929; therefore, we discard the two aforementioned hypotheses as variable polarization mechanisms for it. Another possibility involves the star’s binarity. The unseen B9.8 V companion (Evans 1991) has an unknown orbit. Both a long period orbit and one as short as 55 d that is face-on are viable based on current data (Benedict et al. 2022). A truly face-on orbit would result in a rotation of |$\theta$| with no change in p. Detailed calculations are required to determine if either photospheric reflection or entrained gas are good hypotheses in this case, though even a 55 d period seems like it would not produce a lot of polarization. If it does, further observations could reveal the geometry of the orbit (Brown, McLean & Emslie 1978; Cotton et al. 2020).

HD 210121: The final star displaying a |$\theta$| variability more significant than 2|$\sigma$| is HD 210121 – |$e_{\star \theta }=$| 0.264^{|$^\circ$|}. This star has one of the smallest p values, so this corresponds to an |$e_{\star \, u^\prime }$| of only 126 ppm. If the star really is a bright giant, winds might be responsible. However, this polarization is small enough that many mechanisms could explain such variability. Any of non-radial pulsations, a binary mechanism, emission behaviour, or interstellar drift could be responsible, but no clear leads present themselves in the literature.

5.4 Relative impact of corrections

In Sections 3.1 and 4, we have carried out a series of corrections aimed at more accurately calibrating the telescope position angle for each run. In Table 7, we rank each correction by the median difference between our final determination, and what we’d get neglecting the specified corrections. For this, we include all except the Pindari runs.

Assuming the new values of |$\theta$| reported are the real values, the typical improvement in |$\theta _t$| determination is 0.3|$^\circ$|⁠. The most significant corrections are due to precession and the re-determined |$\theta$| values. Here, we have calculated |$\theta$| for each star in a 2020 equinox. Our observations span a decade, from 2014 to 2023. It can be seen that making precession corrections is now essential to achieve 0.1^{|$^\circ$|} accuracy if working from older literature. However, using contemporary sources will be suitable for most applications.

The largest |$\Delta \theta _{t}$| values for |$\Delta \lambda /\Delta \theta$| (airmass) correspond to runs with Clear rather than |$g^\prime$| observations. For instruments where |$\theta$| rotates with wavelength, narrower bands are better for calibration. If one wants to simplify their reduction routine, this analysis suggests they need not worry about wavelength correction post-fact to account for airmass, but it is a good idea to make band corrections before beginning if the literature data is in a different filter.

Applying the weightings associated with |$e_{i\, \theta }$|⁠, |$e_{\star \theta }$|⁠, and |$e_{\rm m\, \theta }$| is only important in a few cases – those where only a few standard measurements were obtained and one of them was one of the most variable stars. If one observes only the least variable standards, then weighting will generally not be needed. However, if forced to use less reliable standards, one should use the best weightings available, and observe as many different standards as possible – since a standard variable on long time-scales may be invariable on short ones and because |$e_{t\, \theta }$| scales |$\propto \theta _e/\sqrt{N_o}$|⁠.

5.5 Recommendations for observers

We reiterate the best determinations for position angle parameters in Table 8, as they are important for our recommendations to observers seeking better position angle calibration. Both |$\theta$| and its associated variability |$e_{\star \theta }$| are derived in Section 4.3, the other parameters come from the literature as described in Appendix  A. Using this data and applying the advice below will allow for better than 0.1^{|$^\circ$|} error in |$\theta _t$| from a handful of observations in most circumstances. In order of importance our recommendations are:

• Correct |$\theta$| for the equinox of observations using the final column in Table 8 or equation (A1). Report the equinox.

• Use the best available determination of |$\theta$| for your standards, report the value assigned and provenance of the standards employed. Ideally all determinations should be made with the same |$\theta _0$|⁠. For the stars studied here, we recommend the values in Table 8. If other standards are needed, we prefer the values given in Hsu & Breger (1982), and then, if those are unavailable, determinations tied to Serkowski et al. (1975)’s |$\theta _0$| calibration and determined from long run data.

• Use the available |$\Delta \theta /\Delta \lambda$| information in determining |$\theta _{\rm lit}$| for your bandpass. In many cases a full bandpass model will not be required, but if carrying out one for this purpose, use the spectral types and values of |$E_{(B-V)}$|⁠, |$R_{\rm V}$|⁠, |$p_{\rm max}$|⁠, |$\lambda _{\rm max}$|⁠, and K given in Table 1 for the stars studied here, or follow the methodologies laid out in Appendix  B to determine values for others.

• Characterize the position angle error associated with your instrument independent of photon shot noise. To achieve this, one should observe at least one, preferably a few, of the least variable standards multiple times within a few nights, calculate the standard deviation of each, subtract the RMS shot-noise contributions, and take the median. Such a determination should be made and reported as part of commissioning an instrument on a new telescope. Without this information combining data between groups or even just runs to a high precision in |$\theta$| is fraught.

• Choose the most reliable standards available based on reported |$\theta$| variability. For the stars studied here, we recommend the values in Table 8. For other standards, taking |$\sqrt{\sigma _2(\theta)^2-\sigma _1(\theta)^2}$| from Bastien et al. (1988) is recommended. Where |$e_{\star \theta }$| of a potential standard is unknown, the following equation might be used as a first approximation:

  $$\begin{eqnarray} e_{\star \theta } = \frac{28.65}{\sqrt{2}}\frac{e_\star }{p_{\rm ISM}}, \end{eqnarray}$$

  (20)

  where |$e_{\star \theta }$| is in degrees, |$p_{\rm ISM}$| is the interstellar component of the polarization (equal to |$p_{\rm obs}$| for an ideal standard), and |$e_\star$| should be taken as 1130 ppm for an |$\alpha$| Cyg variable, and 300 ppm otherwise, on the basis of the means in Table 6.

• Especially where the available standards might be less reliable, observe multiple standards (ideally |$\ge 3$|⁠) for calibration and weight the result according to the RMS sum of measurement uncertainty and the best estimate of variability.

6 CONCLUSIONS

Long-duration studies of high-polarization stars are rare, leading us to explore the astrophysical implications of our findings. Extreme stars are extremely polarized. The most polarized object in our study, the LBV star HD 160529, is found to also be one of the largest amplitude polarimetric variables. This star has not previously been noted as a polarimetric variable; seemingly it is variable only on longer time-scales.

Large amplitude polarization variability is common and shows no preferred orientation in |$\alpha$| Cygni variables (HD 80558, HD 111613, HD 183143, and HD 198478), a behaviour that may extend to the F0 Ia HD 7927 at lower levels. By analogy with historical studies, clumpy winds are the most likely polarigenic mechanism, but pulsation might also play a role and a dedicated observational program and investigation are recommended.

Later type supergiants and other ordinary stars that happen to be extremely reddened make more reliable standards. The best standards are HD 84810, HD 147084, HD 154445, HD 161056, and HD 161471, for which we attribute variability |$\le 0.123^\circ$|⁠. However, polarization variability may also be present, at lower levels, where we make 2|$\sigma$| detections. There is also evidence for variability in the literature or extended data for other stars in our sample (i.e. HD 43384 and HD 149757). The |$g^\prime$| data analysed might be too thin to capture episodic behaviour in these cases; Hayes (1975) came to a similar conclusion after finding no variability for HD 149757.

The companion to HD 23512 may be a large-amplitude polarimetric variable. HD 187929 displays position angle variability not correlated with the Cepheid phase, an observation that favours a companion in a face-on orbit.

Our results largely vindicate the findings of Bastien et al. (1988), who were criticized for a lack of statistical rigor by Clarke & Naghizadeh-Khouei (1994) in combining data from multiple sources. The methodology we apply alleviates the problems associated with combining data from disparate observing runs. Key to this is the application of a co-ordinate difference matrix, which works by amalgamating difference measurements of pairs of points.

By combining a decade’s worth of data, we have been able to make more accurate and precise determinations of |$\theta$|⁠, as well as estimates of its variability, |$e_{\star \theta }$| which are tabulated in Table 8. The other standard properties, as given in Tables 1 and 8 , have been derived from literature we have curated to minimize the impacts of disagreement between different researchers; these should also be adopted. Along with these improvements, which are quantitatively assessed in Section 5.4, we provide specific recommendations for observers in Section 5.5 that will allow for better than 0.1^{|$^\circ$|} calibration in telescope position angle from a handful of observations in most circumstances.

ACKNOWLEDGEMENTS

This research has made use of the SIMBAD data base, operated at CDS, Strasbourg, France; NASA’s Astrophysics Data System; and the VizieR catalogue access tool, CDS, Strasbourg, France (doi:10.26093/cds/vizier). Based on spectral data retrieved from the ELODIE archive at Observatoire de Haute-Provence (OHP). We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research.

JPM acknowledges support by the National Science and Technology Council of Taiwan under grant NSTC 112-2112-M-001-032-MY3. Funding for the construction of HIPPI-2 was provided by UNSW through the Science Faculty Research Grants Program. We thank the Friends of MIRA for their support.

We thank the former Director of the Australian Astronomical Observatory, Prof. Warrick Couch, the current Director of Siding Spring Observatory, Prof. Chris Lidman, and all of the staff at the AAT for their support of the HIPPI and HIPPI-2 projects. We thank Prof. Miroslav Filipovic for providing access to the Penrith Observatory. We thank all of the additional student volunteers at UNSW, MIRA, WSU, and elsewhere, who assisted with observations at the various observatories. We thank Dr. Wm. Bruce Weaver for useful comments on the manuscript. We thank Dr. Sarbani Basu for help acquiring a reference.

Based in part on data obtained at Siding Spring Observatory. We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present.

Based in part on data obtained at Western Sydney University, Penrith Observatory. We acknowledge the traditional owners of the land on which the WSU Penrith Observatory stands, the Dharug people, and pay our respects to elders past and present.

DATA AVAILABILITY

All of the data for this work are available in a catalogue described in Section 3.

Footnotes

REFERENCES

APPENDIX A: LITERATURE POSITION ANGLE DATA

In this work, a comparison to literature measurements is used to determine the position angle zero-point. We predominantly work in the SDSS |$g^\prime$| band, which for a typically reddened standard observation has |$\lambda _{\rm eff}\approx 470$| nm. Most literature determinations of |$\theta$| have been made in Johnson filters. The closest is B band, which has |$\lambda _{\rm eff}\approx 440$| nm. Here, we predominantly rely on |$\theta$| determinations made in/for that band by three sources (Serkowski et al. 1975; Hsu & Breger 1982; Bagnulo et al. 2017). Combining the work from these sources gives the values in Table A1. HD 161471 is not listed in these works, instead we use Serkowski & Robertson (1969).

The change in position angle with time due to precession, |$\Delta \theta /\Delta t$|⁠, is described by Van De Kamp (1967) as

where |$\Delta t$| is in years and |$(\alpha _0, \delta _0)$| are the initial co-ordinates in degrees. In Table A1, we adjust the value to a 2020 equinox.

Hsu & Breger (1982) and Bagnulo et al. (2017) determine the change in position angle with wavelength, |$\Delta \theta /\Delta \lambda$|⁠; where this is available, or estimates are possible from other sources, we have made an adjustment corresponding to 30 nm. More rigour than this is not justified by the precision of the data, which is nominally 0.2^{|$^\circ$|} for Hsu & Breger (1982)’s position angle determinations, and more than this for |$\Delta \theta /\Delta \lambda$| (in |$^{\circ }/\mu \rm m$|⁠). We round to 0.1^{|$^\circ$|} in Table A1.

APPENDIX B: LITERATURE POLARIZATION DATA

FLC modulator calibration in HIPPI-class polarimeters requires data be checked against the polarization of standard stars. For this we use a bandpass model that requires polarization and extinction parameters for each star. So far this data has been assembled in an ad hoc way (Bailey et al. 2020). Here, we develop a more systematic approach that seeks to favour the most reliable polarization (Serkowski Law) data and homogenizes the available extinction and reddening data on each standard star.

B1 Serkowski law parameters

The intrinsic polarization of high polarization standard stars is assumed to be negligible, and their polarization accomplished entirely by the interstellar medium.¹⁵ Interstellar polarization is described by the empirically determined Serkowski Law (Serkowski 1968)

where |$\lambda$| is the wavelength in |$\mu$|m, p the polarization, |$p_{\rm max}$| the maximum polarization, |$\lambda _{\rm max}$| the wavelength corresponding to |$p_{\rm max}$|⁠, and the dimensionless constant K is the inverse half-width of the curve.

We can determine the most reliable values of |$p_{\rm max}$| and |$\lambda _{\rm max}$| independently of K, since according to Hsu & Breger (1982) and Bagnulo et al. (2017) that parameter does not significantly impact the determination of the other two. If properly calibrated, a spectropolarimeter offers the best precision for this, so long as |$\lambda _{\rm max}$| is within its operating range. For this reason, we prefer determinations made by Martin et al. (1999), Bagnulo et al. (2017), and Wolff et al. (1996) (in that order) where available. Bagnulo et al. (2017) use the FORS2 instrument to take data over the range 0.375–0.940 |$\mu$|m, whereas Wolff et al. (1996) obtain data from 0.400 to 0.700 |$\mu$|m and then supplement this with broadband infrared data from Bailey & Hough (1982) and Nagata (1990). Similarly, Martin et al. (1999) made fits to data from the FOS instrument on the Hubble Space Telescope supplemented with infrared ground based data (of which only Wilking et al. 1980 is relevant here).

Where a Serkowski law fit is instead made to broadband measurements the accuracy of the determinations depends to a greater extent on the quality of the bandpass model used, the reddening law adopted and the bands in which data is obtained. Because |$p \propto 1/\lambda$| in the Serkowski law, infrared data is particularly important to getting a good fit. In acknowledgement of these factors, we next prefer the work of Wilking et al. (1980, 1982), Hsu & Breger (1982), and Serkowski et al. (1975) in that order.

The work of Serkowski et al. (1975) lists the most extensive set of observations. They apply a universal reddening law and are careful to report the details of their bandpass model. However, the longest passband used is R, and their final reported values are actually a straight average of their own work and the earlier determinations of others, who were not always as thorough. Hsu & Breger (1982) obtained multiple observations of each star in |$UBVR$| and a longer 0.75 |$\mu$|m band. They conducted a very careful study – their bandpass model took account of airmass and they also applied a more modern universal reddening law. Wilking et al. (1980, 1982) obtained |$JHK$| band data and added it to optical and UV data of a variety of earlier workers to achieve the widest wavelength range of the broadband works considered here.

The third Serkowski parameter, K, has in some reported cases been fit, and in others assumed. Serkowski et al. (1975) assumed 1.15 – which was the best mean fit to their data and that of Coyne, Gehrels & Serkowski (1974) – this value was also assumed by Hsu & Breger (1982). Wolff et al. (1996) and Martin et al. (1999) applied the relation found by Whittet et al. (1992)

This was the mean determination from fitting |$UBVRIJHK$| data obtained for 105 stars. Earlier Wilking et al. (1980), in fitting K for each individual object, identified a very similar mean relationship

which they revised in later work (Wilking et al. 1982) to

Given the |$\lambda _{\rm max}$| values for our standards, and the typical errors, there is no practical difference between equations (B2), (B3) and (B4).

Bagnulo et al. (2017) also fit K individually for each object; in doing so they obtained quite different values to earlier workers. This could be due to calibration discrepancies or differences in the regions of the ISM probed. However, their favoured explanation is a lack of infrared data, in line with the findings of Clarke & Al-Roubaie (1983), who showed that the relationship recovered by Whittet et al. (1992) depended on the wavelength range probed. This may point to a real, if subtle, feature of interstellar polarization that is masked by fitting a Serkowski curve to data that includes infrared bands.

This presents a difficulty for us in selecting appropriate values of K to use in our calibration. Not all the standards have individually fit values available. Of those that do, the wavelength range of the fit data varies. Re-fitting all of the data in a consistent way is desirable but beyond the scope of this work. Another complicating factor is that different reddening laws, and reddening parameters have been adopted by different workers, which will also impact K.

Our goal is consistency. For this reason, we have decided to adopt the relationship of Whittet et al. (1992) as given in equation (B2), except where K has been fit using a similar wavelength range – in practice, for the standards we have, this means only the Wilking et al. (1980) determinations are applicable. The difference between these values is negligible in half the cases, and very little data has been taken at infrared wavelengths with HIPPI-class instruments anyway, so the difference between the possible approaches here may be academic.

The final Serkowski parameters we adopt for the standards used in this paper, and more broadly are given in Table B1.

B2 Extinction parameters

As described in Bailey et al. (2020), our bandpass model applies a Castelli & Kurucz (2003) atmosphere model based on spectral type.¹⁶ For distant stars, the spectral energy distribution can be reddened using the relationship described by Cardelli et al. (1989). For stars as distant as our standards, this is necessary for a thorough treatment.

Extinction and reddening are related

where |$R_{\rm V}$| is the normalized extinction, |$A_{\rm V}$| is the extinction in the V band, and |$E_{(B-V)}$| the selective extinction. Our algorithm requires |$R_{\rm V}$| and |$E_{(B-V)}$|⁠.

Sometimes the extinction parameters have been reported with Serkowski parameter determinations. However, we have found these to be inconsistent, sometimes superseded by better measurements, and in others simply assumed. We therefore conducted a thorough, though non-exhaustive, search of the literature with the aim of determining appropriate mean values from a reliable sample. We give this data in Table B2. In some cases only |$A_{\rm V}$| has been reported, in such cases we have calculated |$R_{\rm V}$| using the mean |$E_{(B-V)}$| we determine; these are marked with 'c' in the table. The spectral types given in Table B2 are the most common indicated in the cited extinction literature.

Table B2 gives the standard deviation of |$E_{(B-V)}$| and |$R_{\rm V}$| as the 1|$\sigma$| value. With few exceptions, |$E_{(B-V)}$| is well defined. The range of |$R_{\rm V}$| values reported in the literature is much wider; there are also fewer determinations and we have been forced to use whatever we can find. In particular there are few determinations of |$R_{\rm V}$| for the Cepheid variables (HD 84810 and HD 187929). Nevertheless, many of the adopted means match expectations well. The galactic mean is known to be 3.10 (Cardelli et al. 1989) and 12 of 17 stars fall in the range 2.89–3.29. We adopt 3.10 for HD 187929 for which we found no |$R_{\rm V}$| value in the literature.

Of the |$R_{\rm V}$| outliers, HD 147084 and HD 111613 have high values. Whittet & van Breda (1980) report a mean |$R_{\rm V} =$| 3.9 for the Sco-Cen association, of which both these stars are members (Upper Sco and Lower Cen-Cru, respectively; Whittet & van Breda’s measurement for HD 147084 is 3.82). So, the adopted values are reasonable. The other two outliers, HD 210121 and HD 161471, have lower values. Both are amongst the closer high polarization standards at about 450 and 600 pc, respectively (parallax values, as currently listed in SIMBAD, are given in Table B2 for interest’s sake; these values are not used in our calculations). We might expect nearer stars to typically have less extinction since the diffuse Local Hot Bubble will make up more of the sight line to them, but many of the other standards are just as close. The determination for HD 210121 is based on two recent measurements from reliable sources (that are both low), so is likely to be robust.

On the other hand, the only |$A_{\rm V}$| value for HD 161471 comes from Luck (2014), whose |$E_{(B-V)}$| value of 0.47 we had to exclude as an outlier. We therefore examined the extinction maps of Lallement et al. (2019), which indicate lower than typical extinction on HD 161471’s sight line, graphically integrating the extinction density along that sight line, while imprecise, gives answers consistent with Luck (2014)’s |$A_{\rm V}$| determination.

APPENDIX C: MODULATOR RE-CALIBRATION

C1 Original calibration efficiency check

In this work, we aim to assess the long-term variability of established polarization standard stars as observed by HIPPI-class instruments. A number of different modulators, with different characteristics, have been used with these instrument. Consequently, we carried out a preliminary analysis (Table C1) to check the predicted polarization against our observations. The predictions depend on the literature adopted parameters of the stars. The determination of the observed polarization depends on our bandpass model, particularly the modulator efficiency.

The data presented in Table C1 is derived from modulator calibrations presented in Bailey et al. (2015, 2020) and Cotton et al. (2022a) and from standard star data indicated in our past works (Bailey et al. 2015, 2017, 2020, 2023; Cotton et al. 2022b).¹⁷ For this purpose only data taken with the ‘B’ modulator in either the SDSS |$g^\prime$| band or with no filter (Clear), which has a similar effective wavelength, has been used – this is the data relevant to this study. Additionally, because these bands are used for |$\theta$| calibration, there are many more such observations, and limiting this analysis to them is more robust for our purpose here.

In the Table C1, we have grouped together the different performance eras of the boulder nonlinear systems (BNS) modulator, which each have a separate calibration but in some cases represent only a handful of |$g^\prime$| or clear observations. Runs with the Micron Technologies (MT) modulator have been divided by instrument, since these show a clear difference. Table C1 reveals a mismatch between the different modulators. The final line in the table represents data taken with the MT modulator with HIPPI-2 at MIRA, as descibed in Cotton et al. (2022b) and subsequently, where our calibration method varied, as described in the next section C2.

The efficiency differences represent minor sources of error in high precision observations involving small polarizations. Yet, if combining data on high polarization targets as in this work, the impact is significant. It motivated us to homogenize the standard data (Appendix  B), and perform new modulator calibrations.

C2 Recalibration methodology

Ferro-electric Liquid Crystal (FLC) modulators, as used in HIPPI-class polarimeters, have an efficiency as a function of wavelength, |$e(\lambda)$|⁠, that can be described by the terms: |$\lambda _{\rm 0}$|⁠, the wavelength of peak efficiency; |$e_{\rm max}$|⁠, the maximum efficiency; and |$Cd$|⁠, the product of terms describing the crystal birefringence (C), and the plate thickness (d). These terms are related by one of two expressions¹⁸; near total polarization by

and at smaller polarizations (⁠|$\lessapprox$|10  per cent) by

where in both equations

In past work (Bailey et al. 2015, 2020), we have used lab based data derived from injecting a polarized source beam through broad and narrow-band filters and applying equation (C1) to derive |$e_{\rm max}$|⁠; subsequently equation (C2) was used to fit the other parameters, based on multiband observations of high polarization standard stars.

In Cotton et al. (2022b), we switched to fitting |$e_{\rm max}$| along with the other parameters to the observed data. This was possible because more multiband standard data was taken specifically for this purpose. Whereas shorter runs prior to calibration did not allow such complete data sets in the other instances. We now have a greater amount of standard data and are in a position to re-evaluate the earlier calibrations based on a similar approach.

The filters used are fully described in Bailey et al. (2015, 2020) and Cotton et al. (2022a); they include SDSS |$g^\prime$| and |$r^\prime$|⁠, Johnson–Cousins U, V (or a similar substitute), along with a 425 and 500 nm short pass filters (425SP and 500SP, respectively), and a 650 nm long pass filter (650LP). Each was paired with either the B or R PMTs described in Section 3. Note also that HIPPI used SDSS |$g^\prime$| and |$r^\prime$| filters manufactured by Omega Optics, whereas the other instruments used Astrodon Generation 2 equivalents. The two HIPPI-2 instruments used different V band filters, with the southern instrument using a filter with a typical Johnson-like profile, and the northern instrument using a square profiled filter approximating the band, labelled |$V_p$| (see Cotton et al. 2022a).

We use the Python scipy routine curve_fit to match the observed polarizations to predictions using our bandpass model, with |$\lambda _{\rm 0}$|⁠, |$Cd$|⁠, and |$e_{\rm max}$| as the fit parameters. Previously we have error weighted the data, but here equal weighting (nominal 200 ppm error) is given to each observation in recognition of stellar variability likely being a greater source of error than shot noise.

C3 Analysis and adopted parameters

C3.1 Meadowlark

The original Meadowlark calibration (EA¹⁹) was carried out with data from the Gemini run (N2018JUN) and two runs from early 2019 (2019FEB, 2019MAR). To check for any evolution of the modulator, we first fit both the original data set (EA) and that collected since, broken down into two periods (EB and EC) for which we had sufficient multiband data, with the results presented in Table C2. The EA determination of |$Cd$| is different to that previously found from the lab based method (Bailey et al. 2020), but fitting |$e_{\rm max}$| makes the most difference; our lab-based method previously gave us a value close to 1, now we find only 0.936.

Table C2 appears to show evolution of parameters from era EA to EB to EC. However, when we removed two U band observations from the 2019FEB run and refit N2018JUN as E0 and everything else as (E1), a satisfying fit was obtained for E1. The U band data has a significant position angle rotation compared to all the other bands, which almost certainly indicates rotation within the band, suppressing the observed polarization. With those two observations removed, that data is no longer such a good match for the N2018JUN run, which has a redder and quite different fit. The telescope polarization at Gemini North was very large at blue wavelengths and not well fit in position angle, which probably accounts for this. Alternatively, the reflectance of the mirrors with wavelength might not be so well determined, and could be wrapped into the modulator parameters.

C3.2 Boulder nonlinear systems

As previously reported (Bailey et al. 2020), the performance of this modulator evolved over time, particularly rapidly during the 2018AUG run, after which it was removed from service.

We initially fit the data as per the previously delineated performance eras, as shown in Table C3. From this it was noted that only the |$\lambda _{\rm 0}$| values were more than 1.5|$\sigma$| from the error-weighted mean. Only a small number of stars were available for some of the Eras, reducing the robustness of the determinations, so it was decided to refit |$\lambda _{\rm 0}$| for each, but keeping |$Cd$| and |$e_{\rm max}$| fixed to the error-weighted means. In general the reduced |$\chi ^2$| statistic was improved by this, which validates the approach.

The adopted values result in an |$e_{\rm max}$| that is greater than previously established for the BNS modulator. The opposite is true for the ML unit, as expected from Table C1.

C3.3 Micron technologies

The MT modulator is the oldest unit and has seen service on many hundreds of nights. It has operated reliably throughout that time. Yet Table C1 foreshadows the issue we investigate here – divergent performance between the earliest (2014MAY) run and later Mini-HIPPI runs (both at UNSW and Pindari Observatory).

Our analysis is presented in Table C4. The difference between the prior fit to run N2022AUG (E3) and the initial fit results from a combination of the newly adopted standard star parameters and the addition of data from the N2023FEB and N2023MAY runs – which are wholly in the |$g^\prime$| band. We also refitted, but did not tabulate, just the N2022AUG run (⁠|$N_o=15$|⁠) to gauge the effect of the updated stellar parameters in isolation: |$\lambda _{\rm 0} = 503.0 \pm 1.7$|⁠, |$Cd = 1.949 \pm 0.113 \times 10^{7} nm^{3}$|⁠, |$e_{max} = 0.981 \pm 0.005$| – within error of those reported in Cotton et al. (2022b).

Our analysis of other eras is complicated by the data available. Only |$g^\prime$| standard observations were made during the Pindari observatory runs (E2). Only a single |$r^\prime$| standard observation was made during the 2014MAY HIPPI run (E0), the other half dozen observations were either in |$g^\prime$| or Clear. In general there is a dearth of multiband data amongst these runs, and the situation is further complicated by the restricted number of standards available to the two small telescopes – the smallest used at Pindari Observatory in particular. This means we cannot reliably fit |$Cd$| nor |$\lambda _{\rm 0}$| for era E0 nor E2 in isolation. Consequently, we initially fixed those values at those derived from the E1 (UNSW Mini-HIPPI runs) and E3 (MIRA HIPPI-2 runs) to gauge |$e_{\rm max}$|⁠, as shown in Table C4. We fit E0 with both E3- (shown) and E1-seeds (not shown); |$e_{\rm max}$| was similar but the E3 seed resulted in a |$\chi ^2$| half the value.

Mini-HIPPI used a Glan–Taylor prism (Thorlabs GT5-A; Bailey et al. 2017), rather than the Wollaston prism employed by HIPPI/-2, the second beam of which is much less efficient²⁰; this is not otherwise accounted for and is likely the reason for lower |$e_{\rm max}$| for E1 and E2. The similarly high values of |$e_{\rm max}$| for E0 and E3 indicate the modulator itself is unchanged.

Subsequently, we fit eras E1 and E2 together, and eras E0 and E3 together. These secondary fits with larger data sets showed no significant difference between them in terms of |$\lambda _{\rm 0}$| and |$Cd$|⁠. Therefore, we determined their error-weighted averages and fixed those parameters to determine two values to adopt for |$e_{\rm max}$| for the MT modulator, when used with either Mini-HIPPI or HIPPI/-2.

Assuming the primary Glan–Taylor beam to have the same polarization efficiency as each beam in the Wollaston prism therefore implies |$e_{\rm max}$| for the secondary beam is 0.857. Or, put another way, the polarizing efficiency of the secondary beam is 87.7 per cent that of the primary beam – in line with expectations.

C4 Accuracy assessment

Some measures of the correspondence between the standard observations and bandpass predictions after the modulator recalibration are given in Table C5. The agreement is best for |$g^\prime$| observations and the bands nearest to it, owing to the much greater number of observations made in |$g^\prime$|⁠. In aggregate the deviation from the predictions is no more than a few per cent²¹ in any of the regularly used bands, and typically less than one percent for HIPPI and HIPPI-2. These figures compare favourably to those presented for HIPPI-2 in Bailey et al. (2020).

The greatest deviation occurs for U band, which we excluded from our fitting routines because a position angle rotation indicated the polarization was likely suppressed – the tabulated data supports this. The three observations made with the 425SP and 500SP filters (which both have a red leak we have characterized) with the R PMT are also in poor agreement. These were all made when the BNS modulator was evolving quickly, which probably explains this. Deficiencies in the components of the bandpass model will also show up here. For instance, the efficiency curves of the PMTs are based on manufacturer data and not well defined at the extremes of the wavelength range. Aside from the obvious (U, 425SP, 650LP) we also expect the |$r^\prime$| band paired with the B PMT to be slightly less accurate, since it includes a large contribution from the edge of the B PMT’s sensitivity range.

The standard deviation in |$p_{\rm obs}/p_{\rm pred}$| is typically a few percent; this is better than our previous calibration (Bailey et al. 2020) by a factor of two. Contributors to the scatter include inaccuracy and imprecision in the literature determinations – which have only been reported to the nearest 100 ppm (0.01 per cent); see Table B1 – as well as instrumental effects and any intrinsic variability of the stars.

Some of our V band observations are clearly discrepant, and this manifests in a larger standard deviation. Inspection of the filter revealed a small mark near its centre which probably explains this.

The standard deviation is noticeably larger for the 425SP band. For hot and/or luminous stars with significant electron scattering in their atmospheres a greater magnitude of intrinsic polarization is to be expected from some mechanisms. However, bluer wavelengths are also most impacted by seeing noise and airmass changes during an observation. In the case of the instruments used here, the comparative steepness of the modulator efficiency curves in the blue will also be a factor. For most standards |$\lambda _{\rm max}$| is around 0.5–0.6 |$\mu$|m, which means inaccuracy in its determination will be magnified in the shortest wavebands.

Mean and median values for the difference between |$p_{\rm obs}$| and |$p_{\rm pred}$| are also given for comparison to the nominal measurement errors in Table C5. These are only intended as representative, since this metric will be very sensitive to stars with larger polarizations, and the mix of stars is heterogeneous. However, it does show the scale of unaccounted for scatter is similar for HIPPI/-2 and Mini-HIPPI, despite utilizing completely different telescopes and sites and, gives a rough indication of what level of improvement in standard determinations is desirable and/or what intrinsic variability might be present. In particular, it is noteworthy that the median disagreement in the most reliable central bands is about 460 ppm (the number weighted average of the two medians).

APPENDIX D: SELECTED SYMBOLS