# 13 Years of P Cygni Spectropolarimetry: Investigating Mass-loss Through
H$\alpha$, Periodicity, and Ellipticity

Keyan Gootkin Department of Astronomy, University of Washington, Seattle, WA,
USA  Trevor Dorn-Wallenstein Department of Astronomy, University of
Washington, Seattle, WA, USA  Jamie R. Lomax Physics Department, United States
Naval Academy, 572C Holloway Rd, Annapolis, MD, 21402, USA  Gwendolyn Eadie
Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON,
Canada  Department of Statistical Sciences, University of Toronto, Toronto,
ON, Canada  Emily M. Levesque Department of Astronomy, University of
Washington, Seattle, WA, USA  Brian Babler  Department of Astronomy,
University of Wisconsin-Madison, 475 N. Charter St., Madison, WI 53706, USA
Jennifer L. Hoffman Department of Physics and Astronomy, University of Denver,
2112 East Wesley Ave., Denver, CO 80208, USA  Marilyn R. Meade  Department of
Astronomy, University of Wisconsin-Madison, 475 N. Charter St., Madison, WI
53706, USA  Kenneth Nordsieck  Department of Astronomy, University of
Wisconsin-Madison, 475 N. Charter St., Madison, WI 53706, USA  John P.
Wisniewski  Homer L. Dodge Department of Physics and Astronomy, University of
Oklahoma, 440 W Brooks Street, Norman, OK 73019, USA Keyan Gootkin
<EMAIL_ADDRESS>
###### Abstract

We report on over 13 years of optical and near-ultraviolet spectropolarimetric
observations of the famous Luminous Blue Variable (LBV), P Cygni. LBVs are a
critical transitional phase in the lives of the most massive stars, and
achieve the largest mass-loss rates of any group of stars. Using
spectropolarimetry, we are able to learn about the geometry of the near
circumstellar environment surrounding P Cygni and gain insights into LBV mass-
loss. Using data from the HPOL and WUPPE spectropolarimeters, we estimate the
interstellar polarization contribution to P Cygni’s spectropolarimetric
signal, analyze the variability of the polarization across the H$\alpha$
emission line, search for periodic signals in the data, and introduce a
statistical method to search for preferred position angles in deviations from
spherical symmetry which is novel to astronomy. Our data are consistent with
previous findings, showing free-electron scattering off of clumps uniformly
distributed around the star. This is complicated, however, by structure in the
percent-polarization of the H$\alpha$ line and a series of previously
undetected periodicities.

††software: Python 3.7.4, Astropy v4.0.1 (Robitaille et al., 2013; Price-
Whelan et al., 2018), Scipy v.1.18.1 (Jones et al., 2001), Numpy v1.17.2 (van
der Walt et al., 2011),Matplotlib v3.1.1 (Hunter, 2007), REDUCE (Nook et al.,
1990; Wolff et al., 1996; Davidson et al., 2014)

′´

##  1 Introduction

Luminous Blue Variables (LBVs)—also called S Doradus or S Dor variables—are
incredibly rare, with less than 20 confirmed members of this class in the
Milky Way (Richardson & Mehner, 2018). Despite the short lifetime of the LBV
phase ($\sim 10^{5}yr$), these stars are an important set of objects for
understanding the post-main sequence evolution of stars with initial masses
$\geq 20M_{\sun}$ (Groh et al., 2014). This eclectic group represents a
critical transitional phase in the lives of the most massive stars.

LBVs reside at the top of the Hertzsprung-Russell diagram. This unique
position—near the Humphreys-Davidson limit (Humphreys & Davidson, 1979), and
on the S Dor instability strip—produces strange behavior in these stars. LBVs
are highly unstable, losing large amounts of mass, and dramatically varying
both photometrically and spectroscopically. In their quiescent, hot state
(12,000 to 30,000 K), LBVs resemble blue supergiants. However, during an
outburst they move horizontally across the HR diagram to temperatures between
7500 K and 9000 K—maintaining a consistent bolometric luminosity even while
$V$-band magnitudes fluctuate (Humphreys & Davidson, 1994).

Other than Wolf-Rayet stars, LBVs achieve the highest mass-loss rates of any
group of stars, $10^{-5}-10^{-4}M_{\sun}$ $yr^{-1}$, and beyond the “typical”
S Dor outbursts they undergo catastrophic (but non-terminal) eruptions in
which significant mass is lost (Humphreys & Davidson, 1994). The only observed
eruptions in our galaxy were from P Cygni (also known as P Cyg, 34 Cyg, or
Nova Cyg 1600), in which 0.1 $M_{\sun}$ was lost (Smith & Hartigan, 2006), and
$\eta$ Car which is estimated to have lost over 10$M_{\sun}$ of material in
its massive 19th century eruptions (Smith et al., 2003). However, the physical
mechanisms behind the eruption, outburst, and quiescent modes of mass loss are
poorly understood.

The subject of this paper, P Cygni, was discovered on August $18^{th}$, 1600.
The Dutch cartographer, globe-maker, and former student of Tycho Brahe, Willem
Janszoon Blaeu, observed a Nova Stella in the heart of Cygnus (Blaeu, 1602).
He chronicled this discovery in an inscription on a celestial globe, made in
his Amsterdam workshop in 1602111The first edition of this globe, as of the
writing of this paper, is in the collections of Skokloster Castle, a castle
and museum north of Stockholm and can be viewed online here.. This new star
was P Cygni. Although it was called a nova at the time, it was later
recognized as a Luminous Blue Variable (LBV).

During the 1600 eruption (the Nova Stella), P Cygni brightened to the point of
visibility (up to $3^{rd}$ magnitude) for the first time, and remained visible
for 26 years before fading to $\sim 6^{th}$ magnitude. It reached $\sim
3^{rd}$ magnitude again in 1654, experienced variability for a number of
decades, but has been in a stage of quiescence since the late $18^{th}$
century (de Groot & Mart, 1988).

Even in a rare, enigmatic class of objects, P Cygni is unique. It does not
exhibit all of the same photometric behavior of other S Doradus variables. It
has stayed in its quiescent state, at roughly 18,700 K and 61,000 $L_{\sun}$,
for several hundred years. The nebula surrounding P Cyg is faint and
morphologically different from those of most LBVs (Nota et al., 1995). It has
also been studied for longer than any other LBV; an excellent summary of work
done in the $20^{th}$ century can be found in Israelian & de Groot (1999).
Being such a singular and well studied object, P Cygni is critically important
to our understanding of the evolution of the most massive stars in the
universe.

Additionally, since it has been so stable for the past 300 years, ejecta from
the outbursts have had the time to move far away from the star, and
polarimetry can now open a window into the less dramatic, quiescent mode of
mass-loss. There have been extensive studies of the shells of material which
were ejected in the 17th century eruptions (Barlow et al., 1994; Nota et al.,
1995; Meaburn et al., 1996), first resolved by Leitherer & Zickgraf (1987),
but much less is known about the nature of the ambient stellar wind.

In order to study the mass-loss in P Cygni’s current state we turn to
polarimetry. In the study of massive stars, polarimetry is an invaluable tool.
Due to these stars’ enormous luminosities and optically thick nebulae, it is
difficult to observe the innermost regions of their stellar winds. Polarimetry
circumvents this problem, in fact, in certain geometries polarimetric signals
can be enhanced in the presence of an optically thick nebula (Wood et al.,
1996a, b). While P Cygni is in a quiescent (hot) state, we expect electron
scattering to be the dominant mechanism causing any net polarization. This
scattering is only efficient very near to the star (Nordsieck et al., 2001;
Taylor et al., 1991a; Davies et al., 2005, 2006) and, as is discussed in the
following section, is only sensitive to asymmetries when unresolved.
Therefore, polarimetry is a good way to probe asymmetries at the base of P
Cygni’s stellar wind.

While P Cygni had long been known to be polarized (Coyne & Gehrels, 1967;
Coyne et al., 1974; Serkowski et al., 1975), the first extensive study of its
polarization was conducted by Hayes (1985). Taylor et al. (1991a) and
Nordsieck et al. (2001) (hereafter T91 and N01) followed the results of Hayes
(1985) using optical spectropolarimetric observations. Spectropolarimetric
observations measure the Stokes parameters q and u as a function of
wavelength, just as spectroscopic observations measure flux as a function of
wavelength. The following conclusions were made by all three studies
(Nordsieck et al., 2001; Taylor et al., 1991a; Hayes, 1985):

  * •

P Cygni is intrinsically polarized, indicating the presence of asymmetries in
the polarizing region.

  * •

The polarization of P Cygni does not appear to prefer a position angle. This
implies that the asymmetries which produce net polarization are equally likely
to emerge in any direction on the plane of the sky.

  * •

With previous data sets the polarization of P Cygni did not appear to be
periodic, although there was evidence for a 10-15 day characteristic timescale
for polarimetric changes.

As presented in the discussion of T91, these results seem to favor a roughly
axisymmetric wind with regions of enhanced electron density, referred to as
inhomogeneities or clumps, which are ejected from the star and travel outwards
through the polarizing region within its wind. We will use these terms in this
work as well.

In this paper, we combine the previously published data from T91 and N01 with
unpublished observations from the early 2000s, which were taken with the same
instrument as these archival data. In section 2, we discuss the details of
these observations. We discuss our attempts to remove the polarization due to
the interstellar medium in Section 3. We analyze the remaining polarization
signature, which is intrinsic to the P Cyg system, with novel methods in
Section 4. In Section 5, we discuss our results in the larger context of what
is already known about the P Cyg system. Finally, we summarize our conclusions
in Section 6.

##  2 Observations

###  2.1 HPOL

Our data are comprised of spectropolarimetric observations taken on 80 nights,
spanning 13 years, using the University of Wisconsin’s Half-wave
Spectropolarimeter (HPOL) while it was mounted on the 36” telescope at Pine
Bluff Observatory. In 1995 HPOL switched from using a Reticon dual channel
photo-diode array to a $400\times 1200$ pixel CCD. As a result our data are
split into the 46 observations taken before this switch (hereafter Reticon
data) and the 34 taken after (CCD data). The Reticon detector covered a
wavelength range between 3200 and 7600 Å at a 15 Å resolution, while the CCD
recorded data using blue (3200-6000 Å) and red (6000-10,500 Å) gratings in
combination to obtain a full observation. The blue grating had a spectral
resolution of 7.5 Å, while the red grating was 10 Å. More complete reviews of
this instrument are in Nordsieck & Harris (1996), Wolff et al. (1996), and
Davidson et al. (2014). Basic reduction steps were performed on all of the
HPOL observations using the custom built Fortran software package REDUCE (Nook
et al., 1990; Wolff et al., 1996; Davidson et al., 2014).

We downloaded these data from the Mikulski Archive for Space Telescopes (MAST)
archive222The website for this archive is https://archive.stsci.edu and
omitted 4 nights of data from our final data set. These four nights of CCD
observations only used the red grating (27 September 1997, 27 August 1998, 30
May 2000, and 13 November 2002), which does not cover the full wavelength
range we were interested in analyzing. Therefore, we omitted them from our
final data set. Table 1 lists the civil and modified Julian dates of our
remaining HPOL observations (46 Reticon observations and 30 full CCD
observations). Of the full 80 nights of data, T91 used 20 of the Reticon
observations and N01 used 15 of the CCD observations.

From the normalized Stokes parameters, q and u, percent polarization,
$P_{\%}$, and position angle, $\Theta$, are defined as

$P_{\%}=\sqrt{q^{2}+u^{2}}$ (1)

$\Theta=\frac{1}{2}\arctan{\frac{u}{q}}.$ (2)

####  2.1.1 Synthetic V Band Data

To characterize the optical broadband behavior of P Cyg, we convolved each of
our HPOL observations with a synthetic Johnson V-band filter (Bessell, 1990).
We report these values in Table 1 as our “observed” V band values. Figure 1
shows a graphical representation of how this convolution process works
mathematically.

HPOL’s instrumental polarization was monitored on an approximately monthly
basis by observing polarized and unpolarized standard stars. Davidson et al.
(2014) reported the instrument’s systematic uncertainties in the Johnson UBVRI
bands for the CCD detector, which were calculated using those observations. In
order to match HPOL’s historically reported values, we calculated one overall
root-mean-square systematic uncertainty for each CCD observation in the V band
using the individual q and u systematic uncertainties compiled by Davidson et
al. (2014). Two observations (10 September 1998 and 16 September 1998) fell
outside of all of the date ranges for which Davidson et al. (2014) reported
systematic uncertainties. Therefore, we used the values reported for the 14
June 1998-16 August 1998 range because it was closest in time to these
observations. Systematic uncertainties for the Reticon detector are less well
known, but our previous experience with HPOL data suggests they are less than
0.02%. In order to characterize the uncertainty in the V band for a given
observation, we use the larger of its systematic and internal uncertainties,
which are from photon statistics, in our analysis. We report both these
uncertainties in Table 1.

Figure 1: Graphic illustrating how the $V$-band (§2.1.1) and H$\alpha$ $pfew$ (§2.1.2) values tabulated in Table 1 are calculated. These data were taken from a representative observation after ISP correction as described in §3. For the pfew side of the figure (right) the dashed lines represent regions contributing to the continuum polarization and the color represents regions contributing to the line polarization. Table 1: HPOL Observations  |  |  | V Band | H$\alpha$ Continuum | H$\alpha$ Line
---|---|---|---|---|---
Detector | MJD | Date | Observed | Intrinsic |  |  | Intrinsic | Observed | Intrinsic | Observed
|  |  | % q | % u | % q | % u | Int Err | Sys Err | % q | % u | Error | % q | % u | Error | % q | % u | Error | % q | % u | Error
Reticon | 47668 | 1989-05-22 | 0.305 | 0.997 | -0.134 | -0.135 | 0.020 | 0.004 | -0.150 | -0.171 | 0.065 | 0.240 | 0.877 | 0.065 | 0.046 | 0.027 | 0.172 | 0.435 | 1.075 | 0.172
— | 47672 | 1989-05-26 | 0.611 | 1.057 | 0.172 | -0.076 | 0.020 | 0.004 | 0.214 | -0.042 | 0.048 | 0.603 | 1.006 | 0.048 | -0.110 | 0.115 | 0.108 | 0.279 | 1.162 | 0.108
— | 47697 | 1989-06-20 | 0.289 | 1.128 | -0.150 | -0.005 | 0.020 | 0.004 | -0.058 | 0.042 | 0.053 | 0.332 | 1.090 | 0.053 | -0.227 | -0.143 | 0.101 | 0.163 | 0.905 | 0.101
— | 47714 | 1989-07-07 | 0.437 | 0.863 | -0.002 | -0.269 | 0.020 | 0.005 | -0.077 | -0.100 | 0.052 | 0.313 | 0.948 | 0.052 | -0.138 | -0.099 | 0.131 | 0.251 | 0.949 | 0.131
— | 47808 | 1989-10-09 | 0.527 | 1.106 | 0.088 | -0.026 | 0.020 | 0.004 | 0.069 | -0.028 | 0.055 | 0.458 | 1.019 | 0.055 | -0.098 | -0.099 | 0.099 | 0.292 | 0.949 | 0.099
— | 47822 | 1989-10-23 | 0.174 | 0.826 | -0.265 | -0.306 | 0.020 | 0.003 | -0.193 | -0.269 | 0.037 | 0.197 | 0.779 | 0.037 | 0.028 | -0.168 | 0.064 | 0.417 | 0.879 | 0.064
— | 47913 | 1990-01-22 | 0.716 | 1.317 | 0.276 | 0.184 | 0.020 | 0.007 | 0.198 | 0.169 | 0.096 | 0.587 | 1.217 | 0.096 | -0.064 | -0.292 | 0.159 | 0.326 | 0.756 | 0.159
— | 47918 | 1990-01-27 | 0.547 | 1.056 | 0.108 | -0.076 | 0.020 | 0.003 | 0.154 | -0.148 | 0.035 | 0.544 | 0.901 | 0.035 | 0.020 | -0.154 | 0.084 | 0.409 | 0.894 | 0.084
— | 47919 | 1990-01-28 | 0.675 | 1.111 | 0.236 | -0.021 | 0.020 | 0.003 | 0.163 | -0.051 | 0.040 | 0.553 | 0.997 | 0.040 | -0.027 | 0.060 | 0.079 | 0.362 | 1.108 | 0.079
— | 47920 | 1990-01-29 | 0.762 | 1.027 | 0.323 | -0.105 | 0.020 | 0.006 | 0.306 | -0.098 | 0.078 | 0.696 | 0.951 | 0.078 | 0.050 | 0.101 | 0.133 | 0.440 | 1.149 | 0.133
— | 47942 | 1990-02-20 | 0.202 | 1.097 | -0.237 | -0.035 | 0.020 | 0.004 | -0.224 | -0.040 | 0.061 | 0.166 | 1.009 | 0.061 | -0.081 | -0.110 | 0.100 | 0.308 | 0.937 | 0.100
— | 47950 | 1990-02-28 | 0.353 | 1.387 | -0.086 | 0.255 | 0.020 | 0.003 | 0.009 | 0.275 | 0.044 | 0.399 | 1.323 | 0.044 | 0.010 | -0.055 | 0.091 | 0.399 | 0.992 | 0.091
— | 47987 | 1990-04-06 | 0.506 | 1.197 | 0.067 | 0.065 | 0.020 | 0.004 | 0.038 | 0.153 | 0.059 | 0.428 | 1.201 | 0.059 | 0.049 | -0.082 | 0.114 | 0.438 | 0.966 | 0.114
— | 47988 | 1990-04-07 | 0.558 | 1.155 | 0.119 | 0.023 | 0.020 | 0.004 | 0.115 | 0.168 | 0.047 | 0.504 | 1.216 | 0.047 | -0.093 | -0.096 | 0.079 | 0.297 | 0.952 | 0.079
— | 47994 | 1990-04-13 | 0.355 | 0.955 | -0.084 | -0.178 | 0.020 | 0.003 | -0.086 | -0.119 | 0.031 | 0.303 | 0.929 | 0.031 | 0.007 | 0.028 | 0.064 | 0.397 | 1.075 | 0.064
— | 48003 | 1990-04-22 | 0.341 | 1.367 | -0.098 | 0.234 | 0.020 | 0.006 | -0.027 | 0.287 | 0.079 | 0.362 | 1.335 | 0.079 | -0.003 | -0.168 | 0.194 | 0.387 | 0.880 | 0.194
— | 48018 | 1990-05-07 | 0.282 | 0.965 | -0.158 | -0.167 | 0.020 | 0.007 | -0.173 | -0.207 | 0.099 | 0.217 | 0.841 | 0.099 | 0.169 | -0.131 | 0.180 | 0.559 | 0.917 | 0.180
— | 48047 | 1990-06-05 | 0.556 | 0.977 | 0.116 | -0.155 | 0.020 | 0.015 | -0.028 | -0.278 | 0.248 | 0.362 | 0.770 | 0.248 | 0.250 | 0.168 | 0.446 | 0.640 | 1.215 | 0.446
— | 48057 | 1990-06-15 | 0.300 | 1.223 | -0.139 | 0.091 | 0.020 | 0.002 | -0.190 | 0.074 | 0.025 | 0.199 | 1.121 | 0.025 | 0.104 | -0.004 | 0.058 | 0.494 | 1.044 | 0.058
— | 48088 | 1990-07-16 | 0.331 | 1.203 | -0.108 | 0.071 | 0.020 | 0.003 | -0.094 | 0.057 | 0.037 | 0.296 | 1.105 | 0.037 | -0.161 | -0.077 | 0.069 | 0.229 | 0.970 | 0.069
— | 48133 | 1990-08-30 | 0.759 | 1.077 | 0.320 | -0.056 | 0.020 | 0.003 | 0.193 | -0.060 | 0.039 | 0.583 | 0.988 | 0.039 | 0.201 | -0.121 | 0.067 | 0.590 | 0.926 | 0.067
— | 48150 | 1990-09-16 | 0.302 | 1.153 | -0.137 | 0.020 | 0.020 | 0.002 | -0.166 | 0.007 | 0.025 | 0.224 | 1.055 | 0.025 | 0.012 | -0.076 | 0.058 | 0.402 | 0.972 | 0.058
— | 48159 | 1990-09-25 | 0.476 | 1.364 | 0.036 | 0.232 | 0.020 | 0.002 | -0.046 | 0.213 | 0.027 | 0.344 | 1.261 | 0.027 | 0.058 | -0.117 | 0.063 | 0.447 | 0.931 | 0.063
— | 48224 | 1990-11-29 | 0.472 | 0.962 | 0.033 | -0.170 | 0.020 | 0.002 | 0.068 | -0.123 | 0.023 | 0.458 | 0.925 | 0.023 | -0.081 | -0.114 | 0.063 | 0.309 | 0.935 | 0.063
— | 48231 | 1990-12-06 | 0.275 | 0.495 | -0.164 | -0.637 | 0.020 | 0.002 | -0.118 | -0.636 | 0.023 | 0.271 | 0.412 | 0.023 | -0.066 | 0.089 | 0.047 | 0.324 | 1.137 | 0.047
— | 48235 | 1990-12-10 | 0.099 | 1.225 | -0.340 | 0.092 | 0.020 | 0.002 | -0.329 | 0.052 | 0.020 | 0.061 | 1.100 | 0.020 | 0.085 | -0.225 | 0.046 | 0.474 | 0.822 | 0.046
Reticon | 48320 | 1991-03-05 | 0.279 | 1.290 | -0.161 | 0.158 | 0.020 | 0.003 | -0.171 | 0.117 | 0.046 | 0.219 | 1.165 | 0.046 | 0.002 | -0.090 | 0.088 | 0.392 | 0.958 | 0.088
— | 48444 | 1991-07-07 | 0.523 | 1.286 | 0.084 | 0.154 | 0.020 | 0.005 | 0.133 | 0.046 | 0.068 | 0.523 | 1.094 | 0.068 | -0.086 | 0.007 | 0.128 | 0.303 | 1.055 | 0.128
— | 48473 | 1991-08-05 | 0.537 | 0.973 | 0.098 | -0.159 | 0.020 | 0.014 | 0.348 | -0.028 | 0.264 | 0.738 | 1.020 | 0.264 | 0.036 | 0.094 | 0.360 | 0.426 | 1.141 | 0.360
— | 48516 | 1991-09-17 | 0.596 | 1.150 | 0.156 | 0.018 | 0.020 | 0.002 | 0.128 | 0.005 | 0.028 | 0.518 | 1.053 | 0.028 | -0.034 | 0.032 | 0.045 | 0.355 | 1.080 | 0.045
— | 48520 | 1991-09-21 | 0.135 | 0.817 | -0.305 | -0.315 | 0.020 | 0.002 | -0.275 | -0.369 | 0.027 | 0.115 | 0.679 | 0.027 | 0.022 | -0.032 | 0.046 | 0.411 | 1.016 | 0.046
— | 48522 | 1991-09-23 | 0.024 | 1.021 | -0.415 | -0.111 | 0.020 | 0.002 | -0.389 | -0.187 | 0.025 | 0.001 | 0.861 | 0.025 | 0.016 | -0.142 | 0.040 | 0.406 | 0.906 | 0.040
— | 48525 | 1991-09-26 | 0.799 | 1.336 | 0.360 | 0.203 | 0.020 | 0.004 | 0.277 | 0.197 | 0.060 | 0.666 | 1.245 | 0.060 | -0.122 | -0.036 | 0.098 | 0.268 | 1.012 | 0.098
— | 48526 | 1991-09-27 | 0.843 | 1.348 | 0.404 | 0.216 | 0.020 | 0.002 | 0.429 | 0.112 | 0.036 | 0.819 | 1.160 | 0.036 | -0.137 | -0.010 | 0.053 | 0.253 | 1.038 | 0.053
— | 48527 | 1991-09-28 | 0.796 | 1.226 | 0.357 | 0.094 | 0.020 | 0.002 | 0.367 | 0.062 | 0.025 | 0.757 | 1.110 | 0.025 | -0.111 | -0.085 | 0.040 | 0.278 | 0.963 | 0.040
— | 48528 | 1991-09-29 | 0.740 | 1.089 | 0.301 | -0.044 | 0.020 | 0.002 | 0.295 | -0.046 | 0.027 | 0.685 | 1.002 | 0.027 | -0.039 | 0.058 | 0.041 | 0.351 | 1.105 | 0.041
— | 48529 | 1991-09-30 | 0.659 | 1.038 | 0.220 | -0.094 | 0.020 | 0.003 | 0.171 | -0.148 | 0.030 | 0.560 | 0.900 | 0.030 | -0.043 | 0.006 | 0.055 | 0.346 | 1.053 | 0.055
— | 48530 | 1991-10-01 | 0.669 | 1.047 | 0.230 | -0.086 | 0.020 | 0.002 | 0.195 | -0.108 | 0.030 | 0.585 | 0.940 | 0.030 | -0.017 | -0.109 | 0.051 | 0.372 | 0.938 | 0.051
— | 49521 | 1994-06-18 | 0.432 | 0.732 | -0.007 | -0.400 | 0.020 | 0.008 | -0.093 | -0.368 | 0.103 | 0.296 | 0.680 | 0.103 | 0.092 | -0.205 | 0.201 | 0.481 | 0.843 | 0.201
— | 49530 | 1994-06-27 | 0.551 | 0.790 | 0.112 | -0.342 | 0.020 | 0.005 | 0.077 | -0.236 | 0.064 | 0.467 | 0.812 | 0.064 | 0.169 | -0.273 | 0.299 | 0.558 | 0.775 | 0.299
— | 49539 | 1994-07-06 | 0.770 | 0.694 | 0.331 | -0.438 | 0.020 | 0.011 | 0.232 | -0.491 | 0.127 | 0.621 | 0.557 | 0.127 | 0.191 | 0.325 | 0.684 | 0.580 | 1.373 | 0.684
— | 49573 | 1994-08-09 | 0.222 | 0.974 | -0.217 | -0.158 | 0.020 | 0.006 | -0.134 | -0.082 | 0.088 | 0.256 | 0.965 | 0.088 | 0.019 | -0.248 | 0.154 | 0.409 | 0.799 | 0.154
— | 49606 | 1994-09-11 | 0.156 | 1.221 | -0.283 | 0.089 | 0.020 | 0.004 | -0.281 | 0.106 | 0.067 | 0.109 | 1.154 | 0.067 | 0.137 | -0.188 | 0.156 | 0.527 | 0.860 | 0.156
— | 49636 | 1994-10-11 | 0.208 | 1.059 | -0.232 | -0.074 | 0.020 | 0.003 | -0.307 | -0.056 | 0.041 | 0.083 | 0.992 | 0.041 | -0.077 | -0.158 | 0.101 | 0.312 | 0.890 | 0.101
— | 49666 | 1994-11-10 | 0.159 | 1.047 | -0.280 | -0.085 | 0.020 | 0.003 | -0.282 | -0.066 | 0.041 | 0.108 | 0.982 | 0.041 | -0.003 | -0.263 | 0.133 | 0.387 | 0.784 | 0.133
— | 49690 | 1994-12-04 | 0.142 | 0.901 | -0.297 | -0.231 | 0.020 | 0.004 | -0.352 | -0.195 | 0.054 | 0.038 | 0.853 | 0.054 | 0.068 | -0.227 | 0.134 | 0.457 | 0.820 | 0.134
CCD | 49777 | 1995-03-01 | 0.291 | 1.241 | -0.125 | 0.121 | 0.006 | 0.005 | -0.084 | 0.170 | 0.032 | 0.304 | 1.213 | 0.032 | 0.037 | 0.002 | 0.058 | 0.427 | 1.053 | 0.058
— | 49779 | 1995-03-03 | 0.308 | 1.140 | -0.109 | 0.018 | 0.006 | 0.003 | -0.086 | 0.085 | 0.018 | 0.301 | 1.127 | 0.018 | -0.044 | -0.035 | 0.028 | 0.346 | 1.015 | 0.028
— | 49785 | 1995-03-09 | 0.051 | 1.169 | -0.366 | 0.048 | 0.006 | 0.002 | -0.358 | 0.041 | 0.019 | 0.029 | 1.084 | 0.019 | -0.093 | 0.051 | 0.041 | 0.298 | 1.102 | 0.041
— | 49788 | 1995-03-12 | 0.285 | 1.577 | -0.132 | 0.455 | 0.006 | 0.003 | -0.114 | 0.363 | 0.035 | 0.273 | 1.406 | 0.035 | -0.085 | -0.041 | 0.053 | 0.305 | 1.010 | 0.053
— | 49891 | 1995-06-23 | 0.600 | 1.110 | 0.184 | -0.010 | 0.006 | 0.003 | 0.226 | 0.049 | 0.021 | 0.614 | 1.092 | 0.021 | -0.026 | -0.025 | 0.067 | 0.365 | 1.027 | 0.067
— | 49913 | 1995-07-15 | 0.339 | 1.247 | -0.077 | 0.127 | 0.010 | 0.003 | -0.004 | 0.145 | 0.017 | 0.384 | 1.188 | 0.017 | 0.012 | -0.059 | 0.042 | 0.403 | 0.993 | 0.042
— | 49940 | 1995-08-11 | 0.283 | 1.055 | -0.133 | -0.066 | 0.010 | 0.003 | -0.072 | 0.023 | 0.017 | 0.316 | 1.066 | 0.017 | -0.017 | 0.021 | 0.050 | 0.375 | 1.074 | 0.050
— | 49965 | 1995-09-05 | 0.482 | 1.381 | 0.066 | 0.261 | 0.010 | 0.003 | 0.030 | 0.249 | 0.019 | 0.418 | 1.293 | 0.019 | -0.063 | -0.032 | 0.076 | 0.328 | 1.019 | 0.076
— | 49977 | 1995-09-17 | 0.666 | 0.961 | 0.249 | -0.160 | 0.010 | 0.002 | 0.262 | -0.226 | 0.017 | 0.650 | 0.817 | 0.017 | -0.138 | -0.007 | 0.050 | 0.254 | 1.047 | 0.050
CCD | 49992 | 1995-10-02 | 0.210 | 1.049 | -0.207 | -0.071 | 0.010 | 0.002 | -0.192 | -0.098 | 0.013 | 0.196 | 0.944 | 0.013 | -0.024 | 0.053 | 0.041 | 0.368 | 1.108 | 0.041
— | 50009 | 1995-10-19 | 0.242 | 1.008 | -0.174 | -0.112 | 0.010 | 0.004 | -0.158 | -0.052 | 0.018 | 0.230 | 0.991 | 0.018 | 0.010 | 0.014 | 0.052 | 0.402 | 1.068 | 0.052
— | 50036 | 1995-11-15 | 0.382 | 1.181 | -0.034 | 0.060 | 0.010 | 0.002 | -0.018 | 0.051 | 0.014 | 0.370 | 1.094 | 0.014 | -0.083 | -0.003 | 0.051 | 0.309 | 1.051 | 0.051
— | 50055 | 1995-12-04 | 0.289 | 1.222 | -0.127 | 0.101 | 0.010 | 0.003 | -0.132 | 0.148 | 0.016 | 0.256 | 1.190 | 0.016 | -0.057 | -0.089 | 0.061 | 0.334 | 0.964 | 0.061
— | 50733 | 1997-10-12 | 0.312 | 1.077 | -0.105 | -0.043 | 0.009 | 0.003 | -0.094 | -0.098 | 0.014 | 0.293 | 0.944 | 0.014 | 0.165 | 0.193 | 0.084 | 0.562 | 1.259 | 0.084
— | 50803 | 1997-12-21 | 0.476 | 1.173 | 0.059 | 0.052 | 0.009 | 0.002 | 0.074 | 0.070 | 0.014 | 0.461 | 1.111 | 0.014 | -0.037 | -0.304 | 0.073 | 0.359 | 0.760 | 0.073
— | 51010 | 1998-07-16 | 0.411 | 1.045 | -0.006 | -0.075 | 0.010 | 0.002 | -0.015 | -0.106 | 0.013 | 0.372 | 0.936 | 0.013 | -0.093 | -0.022 | 0.068 | 0.302 | 1.040 | 0.068
— | 51034 | 1998-08-09 | 0.711 | 0.966 | 0.294 | -0.155 | 0.010 | 0.003 | 0.332 | -0.146 | 0.015 | 0.719 | 0.896 | 0.015 | -0.323 | 0.256 | 0.075 | 0.072 | 1.319 | 0.075
— | 51066 | 1998-09-10 | 0.470 | 1.170 | 0.053 | 0.049 | 0.010 | 0.002 | 0.072 | 0.018 | 0.013 | 0.460 | 1.060 | 0.013 | 0.043 | 0.151 | 0.067 | 0.438 | 1.211 | 0.067
— | 51072 | 1998-09-16 | 0.390 | 1.203 | -0.027 | 0.082 | 0.010 | 0.002 | 0.018 | 0.073 | 0.015 | 0.405 | 1.115 | 0.015 | 0.023 | -0.127 | 0.067 | 0.418 | 0.934 | 0.067
— | 51160 | 1998-12-13 | 0.518 | 1.277 | 0.101 | 0.156 | 0.008 | 0.003 | 0.079 | 0.150 | 0.014 | 0.467 | 1.193 | 0.014 | -0.024 | -0.052 | 0.048 | 0.369 | 1.002 | 0.048
— | 51177 | 1998-12-30 | 0.265 | 1.069 | -0.152 | -0.052 | 0.008 | 0.004 | -0.171 | -0.034 | 0.028 | 0.217 | 1.009 | 0.028 | -0.012 | 0.085 | 0.059 | 0.378 | 1.136 | 0.059
— | 52533 | 2002-09-16 | 0.210 | 0.970 | -0.206 | -0.151 | 0.012 | 0.004 | -0.219 | -0.067 | 0.024 | 0.169 | 0.976 | 0.024 | -0.151 | -0.012 | 0.050 | 0.239 | 1.039 | 0.050
— | 52534 | 2002-09-17 | 0.089 | 0.949 | -0.328 | -0.171 | 0.012 | 0.004 | -0.367 | -0.111 | 0.030 | 0.021 | 0.932 | 0.030 | 0.126 | 0.056 | 0.056 | 0.517 | 1.107 | 0.056
— | 52541 | 2002-09-24 | 0.236 | 1.029 | -0.181 | -0.091 | 0.012 | 0.004 | -0.219 | -0.047 | 0.027 | 0.169 | 0.996 | 0.027 | -0.153 | -0.095 | 0.058 | 0.238 | 0.956 | 0.058
— | 52542 | 2002-09-25 | 0.237 | 1.031 | -0.179 | -0.090 | 0.012 | 0.006 | -0.248 | -0.033 | 0.030 | 0.140 | 1.010 | 0.030 | -0.053 | -0.018 | 0.075 | 0.338 | 1.034 | 0.075
— | 52561 | 2002-10-14 | 0.589 | 1.035 | 0.172 | -0.085 | 0.012 | 0.004 | 0.150 | -0.078 | 0.022 | 0.538 | 0.965 | 0.022 | -0.181 | -0.005 | 0.056 | 0.210 | 1.047 | 0.056
— | 52562 | 2002-10-15 | 0.487 | 1.117 | 0.070 | -0.003 | 0.012 | 0.004 | 0.096 | 0.014 | 0.022 | 0.483 | 1.056 | 0.022 | -0.078 | -0.072 | 0.048 | 0.312 | 0.979 | 0.048
— | 52579 | 2002-11-01 | 0.293 | 0.919 | -0.123 | -0.201 | 0.012 | 0.004 | -0.121 | -0.137 | 0.025 | 0.267 | 0.906 | 0.025 | -0.070 | -0.112 | 0.055 | 0.321 | 0.939 | 0.055
— | 52583 | 2002-11-05 | 0.176 | 0.953 | -0.241 | -0.168 | 0.012 | 0.005 | -0.281 | -0.170 | 0.079 | 0.107 | 0.873 | 0.079 | -0.082 | 0.031 | 0.147 | 0.309 | 1.083 | 0.147
— | 52586 | 2002-11-08 | 0.449 | 0.839 | 0.032 | -0.282 | 0.012 | 0.004 | -0.081 | -0.204 | 0.025 | 0.307 | 0.839 | 0.025 | 0.023 | -0.134 | 0.043 | 0.414 | 0.916 | 0.043

####  2.1.2 H$\alpha$: the pfew Method

In an attempt to quantify the overall behavior of the H$\alpha$ emission line,
we calculated its polarized flux equivalent width (pfew) (Lomax et al., 2015;
Hoffman et al., 1998), which integrates the Stokes $q$ and $u$ values in the
emission line and removes underlying continuum polarization (see Figure 1 for
a graphical representation of the mathematics used the pfew method). These
values are reported in Table 1. In order to do this, we identified two
continuum regions on either side of the line that the pfew method uses to
determine the line’s underlying continuum: 6400 to 6475 Å on the blue side and
6750 to 6800 Å on the red side. Additionally, we integrated over the line
between 6540 and 6600 Å.

In the remainder of this work we use the term pfew to refer to both this
method and the values it produced—listed under H$\alpha$ Continuum and
H$\alpha$ Line columns in Table 1.

###  2.2 WUPPE

We also used three ultraviolet (UV) spectropolarimetric observations of P Cyg
taken with the Wisconsin UV Photo-Polarimeter Experiment (WUPPE). WUPPE had a
resolution of approximately 12 Å and recorded spectropolarimetric data between
1400 and 3200 Å. P Cyg was observed on four separate occasions with WUPPE,
which flew aboard both the Space Shuttle Columbia as part of the STS-35
ASTRO-1 mission and the Space Shuttle Endeavour as part of the STS-67 ASTRO-2
mission. More information about this instrument can be found in Bjorkman et
al. (1993) and Nordsieck et al. (1994).

The ASTRO-1 data were taken on 5 December 1990 and originally published in
Taylor et al. (1991b). The three ASTRO-2 observations, taken on 3, 8, and 12
March 1995, show little to no variability suggesting P Cyg has a constant UV
polarimetric behavior. Therefore, we did not include the ASTRO-1 observation
in our analysis.

##  3 Interstellar Polarization

Figure 2: The stars within $\sim 3^{\circ}$ of P Cygni found in Heiles (1999).
The length of the line going through the position of each star represents the
relative strength of the star’s polarization and the orientation of the line
represents position angle. The yellow star shows the position of P Cygni.

In general, we expect the light coming from spatially unresolved stars to be
unpolarized. It only becomes polarized when those photons scatter off of
something, which is typically either circumstellar dust around a star, free
electrons in the circumstellar medium (CSM) around a star, or dust in the
interstellar medium. Because P Cyg is a hot star, its circumstellar material
is highly ionized. Therefore, the polarization we observe is from a
combination of electron scattering within P Cyg’s CSM and scattering with dust
in the interstellar medium between P Cyg and Earth. Before analyzing our
spectropolarimetric observations of P Cyg, we separate polarization which is
intrinsic to the P Cyg system due to electron scattering in the CSM (hereafter
called intrinsic polarization) from that which is produced by the interstellar
medium (hereafter called interstellar polarization or ISP).

Electron scattering is a gray process that produces a constant polarization
signal in both q and u (or $P_{\%}$ and $\Theta$) with wavelength. However, a
polarization signal from electron scattering can change with time if the
geometry of the scattering region (i.e. region where the free electrons are
located) is changing or the number of free electrons for photons to scatter
off of changes. Therefore, the polarization signal due to electron scattering
from P Cyg may change location in the $q$ vs. $u$ plane (i.e. Figure 3) with
time, but the polarization signal at different wavelengths will lie on the
same point if they were taken simultaneously.

Conversely, the ISP has a wavelength dependence, but is generally thought to
be constant over long periods of time. In the $q$ vs. $u$ plane (i.e. Figure
3) the interstellar polarization signal manifests itself as a vector that
changes with wavelength which adds to the constant electron scattering signal
from P Cyg. In particular, it is the $P_{\%}$ of the ISP that changes with
wavelength, while the $\Theta$ of the ISP is constant with wavelength
(Serkowski et al., 1975; Wilking et al., 1982). Thus one must construct a
model of the wavelength dependence of the $P_{\%}$ of the ISP along the line
of sight to P Cygni—which is caused by light scattering off of dust grains in
the interstellar medium—to subtract from the observed polarization in order to
determine P Cygni’s intrinsic polarization. The wavelength dependence of the
ISP can be modeled with the Serkowski law (Serkowski et al., 1975), as
modified by Wilking et al. (1982):

$P_{\lambda}=P_{max}e^{-K\ln^{2}{\lambda_{max}/\lambda}}$ (3)

where $P_{max}$ is the largest value of $P_{\%}$ due to the ISP,
$\lambda_{max}$ is the wavelength at which that peak occurs, and,

$K=(0.01\pm 0.05)+(1.66\pm 0.09)\lambda_{max},$

where $\lambda_{max}$ is measured in microns (Whittet et al., 1992).

Therefore, to determine the Serkowski law in the direction of P Cygni, we must
find:

  1. 1.

at least one wavelength at which we know the magnitude of the ISP—This allows
us to solve for $P_{max}$ in equation 3.

  2. 2.

the wavelength of peak polarization, $\lambda_{max}$

A summary of techniques and challenges relating to determining the ISP can be
found in Wisniewski et al. (2010).

Traditionally, the most common method to find a known ISP value is to observe
stars near the science target in the sky and find the average direction and
strength of their polarization and take this to be the value of the ISP.
Unfortunately, the Cygnus region is notoriously complex. This is illustrated
in Figure 2, which shows the relative strength and orientation of linear
polarization for the stars in the Heiles (1999) polarimetric catalog within
$\sim 3^{\circ}$ of P Cygni. Taking only the subset of stars at a similar
distance as P Cygni does not help to establish a pattern. It is worth noting
that the Heiles (1999) catalogue does not provide the wavelength at which
polarization was measured. However, since we do not expect a significant
wavelength dependence in the observed position angle of the ISP, it is
generally still valuable to compare nearby stars in this manner. In the case
of P Cygni, the large variations in the position angles of its surrounding
field stars (Figure 2) make this method of determining the ISP less than
ideal.

Figure 3: 76 synthetic $V$-band HPOL observations (see §2.1.1), observed
H$\alpha$ line polarization (see §2.1.2), and our estimate of the ISP at
H$\alpha$ (see §3; Table 1). Figure 4: Comparison of ISP estimates from
Nordsieck et al. (2001), Taylor et al. (1991a), and this paper. These models
are plotted over the mean HPOL CCD data in black (see Section 3), and the mean
of 3 UV observations from the WUPPE polarimeter in grey. The zoom-in in the
upper-righthand corner shows the region around H$\alpha$ in greater detail.

Alternatively, the line center method makes the assumption that strong
emission lines, such as H$\alpha$, are intrinsically unpolarized. In P Cygni,
these emission lines are formed between 3.2 and 9.3 $R_{\star}$ (Avcioğlu,
1984), beyond the region where free-electron scattering is efficient; i.e. the
H$\alpha$ line forms farther away from the star than the region where electron
scattering occurs (Nordsieck et al., 2001; Taylor et al., 1991a; Davies et
al., 2005, 2006). Therefore, the observed polarization of these lines are
likely a direct measurement of the ISP, which we do not expect to vary on the
timescales over which our data were taken. We find evidence for this in Figure
3 which shows the distribution of the 76 synthetic $V$-band observations and
measurements of the polarization in the H$\alpha$ emission line (§2.1.2) in
$q$-$u$ space. It can be seen that H$\alpha$ varies less than the $V$-band
observations ($\sigma_{V}\simeq 0.2$ while $\sigma_{H}\simeq 0.1$), suggesting
it is less affected by a time-variable electron scattering component (i.e.
polarization intrinsic to P Cyg). This makes H$\alpha$ a better estimator of
the ISP, which likely lies near the center of both distributions.

Both T91 and N01 presented models for the ISP, constructed using different
methods and data. To estimate $\lambda_{max}$, T91 made use of the
relationship between the ratio of total-to-selective extinction, $R_{V}$, and
$\lambda_{max}$ in Serkowski et al. (1975):

$R_{V}=(5.6\pm 0.3)\lambda_{max}$ (4)

and assumed the galactic average, $R_{v}=3.1$ (Cardelli et al., 1989), to find
a peak polarization wavelength of 5500 Å. They then employed the line center
method and assumed that H$\alpha$ is intrinsically unpolarized to estimate
$P_{max}$ (1.06%) using their spectropolarimetric observations of P Cygni.
With both $P_{max}$ and $\lambda_{max}$ they were able to construct a model of
the ISP in the direction of P Cygni.

N01 followed a method to similar to T91’s to estimate the ISP. They also
assumed the H$\alpha$ was intrinsically unpolarized and estimated the ISP at
that wavelength. They used UV data from the Wisconsin Ultraviolet Photo-
Polarimeter Experiment (WUPPE) and assumed that the region between 1700-1900 Å
is also unpolarized to get an estimate of the ISP at those wavelengths. With
the combination of those two data points (the H$\alpha$ and UV ISP) they were
able to use the Serkowski law to back out both a $P_{max}$ ($1.17\pm 0.03\%$)
and a $\lambda_{max}$ (5100 Å).

When we subtracted the ISP estimates from T91 and N01 from our HPOL data, we
found that the resulting polarization was strongly wavelength dependent for
most observations. Because P Cygni is a hot star, we expect that its intrinsic
polarization is due to electron scattering, which is a ‘gray’ process that
should have no wavelength dependence. Therefore, we determined our own ISP
estimate to try to resolve this issue.

Using the line center method (despite the problems discussed in §5) we
initially performed a similar analysis, using $pfew$ measurements of
H$\alpha$, eq. 4 and $R_{v}=3.1$. However, we found that the peak wavelength
which best describes our overall dataset is bluer than both T91’s and N01’s
previously used $\lambda_{max}$ estimates (5500Å and 5100Å). However, both T91
and N01 also noted that they observed a persistent decrease of P Cygni’s
intrinsic polarization into the infrared whose strength varies, suggesting
there might be a high amount of free-free absorptive opacity due to clumps at
the base of the system’s wind.

To more rigorously determine $\lambda_{max}$ instead of assuming that
$R_{V}=3.1$, we first found the mean $P_{\%}$ curve of all our CCD
observations. We calculated this by finding the error-weighted mean Stokes q
and u parameters at every wavelength for which we have CCD data and then
combined those mean Stokes parameters into $P_{\%}$ (Figure 4). We then
performed a least squares fit, using equation 3 and SciPy’s curve_fit function
(Jones et al., 2001), to two sections of our mean $P_{\%}$ curve (between 4000
and 5900 Å and between 6200 and 9000 Å) simultaneously, which were chosen to
avoid including data that fell near the edges of the CCD detector that
typically have higher uncertainties. As part of this process, we allowed
$P_{max}$ and $\lambda_{max}$ in equation 3 to vary and found
$\lambda_{max}=4712\pm 51$ Å.

We then estimated the $P_{\%}$ at H$\alpha$ because that should only be due to
the ISP. We used the pfew H$\alpha$ line values tabulated above, found the
mean H$\alpha$ Stokes q and u value weighted inversely by variance, and
converted that into a mean $P_{\%}$ at H$\alpha$. We used this value, equation
3, and our previously determined $\lambda_{max}$ (4712Å) to find
$P_{\%,ISP}=1.06\%$. Because the ISP is not expected to have a wavelength
dependence in its position angle, we found the average position angle of all
of our observations and assumed that was representative of the ISP
($\Theta_{ISP}=34.4^{\circ}\pm 0.7^{\circ}$).

We then use equation 4 to solve for $R_{V}$ in the direction of P Cygni. Given
our value of $\lambda_{max}$ and equation 4, we find that $R_{V}=2.638\pm
0.028$. This lower, when compared to T91, value of $R_{V}$—corresponding to a
smaller size of dust grains in the interstellar medium—is in excellent
agreement with the value of $R_{V}=2.73\pm 0.23$ determined by Turner et al.
(2001). This result can be explained by the dust being processed by ionizing
radiation from P Cygni and other massive stars in the OB association (Draine,
2003), and suggests we have a robust ISP estimate to subtract from our
observations. Our resulting ISP fit is plotted in yellow in Figure 4. We also
over plot the T91 and N01 ISP estimates for comparison.

To be clear, our ISP model itself is not a fit to the data. The fitting done
above is only for determining $\lambda_{max}$. Using $\lambda_{max}$ and the
polarization at H$\alpha$ (which we assume to be the same as the value of the
ISP at H$\alpha$) we solve for $P_{max}$. Plugging $\lambda_{max}$ and
$P_{max}$ into equation 3 yeilds the wavelength dependence of the ISP signal.

It is worth noting that even though we do not use the WUPPE data, our ISP
estimate agrees well with N01’s assumption that the 1700-1900 Å region of P
Cygni is intrinsically unpolarized (i.e. our estimate overlaps with the WUPPE
data in this wavelength regime). The small vertical offset from our ISP fit
can be explained by P Cygni being a “super-Serkowski” object as discussed in
Anderson et al. (1996).

Figure 5: $q$ vs. $u$, for 76 intrinsic synthetic $V$-band observations of P
Cygni after subtracting the ISP estimate described in §3. The yellow square
marks the center of the distribution. The error of that mean location is
smaller than the marker. While the separation between the mean of the $V$-band
observations and the origin is statistically significant it is perhaps most
notable how small this separation is. This indicates symmetry at the base of
the wind in the P Cyg system—at least over the timescale of these
observations. Figure 6: Mean intrinsic $P_{\%}$, $q$, and $u$ taken at each
$\lambda$ bin over time. Means of percent polarization and position angle were
calculated from the means of $q$ and $u$. The grey background indicates
uncertainty of the mean measured at each $\lambda$ across all observations.
$\Theta$ values are not included; since most $P_{\%}$ values are near the
origin, position angle becomes highly uncertain.

##  4 Results

Figure 7: Slopes of lines fit to the $P_{\%}$ of each observation between 4000
and 7500 Å plotted against time with 1$\sigma$ uncertainties. The dashed
horizontal line indicates a slope of 0. The dashed vertical line shows the
cutoff between Reticon observations (on the left) and CCD observations (on the
right). The top panel shows the model (green) fit to the data (black) for the
observation marked with a green dot. Figure 8: The mean, calculated over all
CCD observations, of intrinsic $P_{\%}$, $\Theta$, Stokes $q$, and Stokes $u$
across H$\alpha$. The grey bounded region, between 6540 and 6600 Å shows the
line core from which the ISP at H$\alpha$ was calculated. Figure 9: H$\alpha$
line-profile in percent polarization for 30 CCD observations. Each observation
has a constant value added. The observations are stacked in chronological
order, with the oldest being on the bottom. The grey bounded region, between
6540 and 6600 Å shows the line core from which the ISP at H$\alpha$ was
calculated. Note—the constant added does not scale linearly with time.

We used our ISP estimate determined above ($\lambda_{max}=$ 4712 Å) and
$P_{\%,ISP}=1.06\%$) and equation 3 to calculate the wavelength dependence of
the ISP. We converted that wavelength dependent $P_{\%}$ and
$\Theta_{ISP}=34.4^{\circ}$ into wavelength dependent $q$ and $u$ values for
the ISP. Subtracting those $q$ and $u$ components of the ISP from each
observation results in an estimate of the intrinsic polarization of P Cygni.
After we subtracted the ISP, we applied synthetic V-band filters to our data
(see §2.1.1) to quantify P Cyg’s intrinsic continuum polarization behavior
(Table 1). After this subtraction, the center of the distribution of the
V-band data in $q$-$u$ space is shifted very near to the origin, as can be
seen in Figure 5. The distance between the center of that distribution and the
origin is $0.047\pm 0.002$.

Additionally, we found the intrinsic polarization remaining in P Cyg’s
H$\alpha$ line and continuum (§2.1.2; Table 1). As discussed in §3,
polarization in the H$\alpha$ line appears to be at least partially suppressed
when compared to the V-band, making it a good estimator of the ISP. However,
we note that the H$\alpha$ line does show some variability with time, which
suggests that not 100% of its observed polarization is due to the ISP. Because
we used the mean H$\alpha$ polarization to generate our ISP estimate, the
possibility that we have over or under subtracted the ISP exists. Without a
more reliable way to determine an ISP estimate we cannot evaluate how much of
an effect this has on our intrinsic polarization of P Cyg, but the agreement
of our $R_{v}$ calculated from our ISP estimate with $R_{v}$ values for P Cyg
in the literature suggests our ISP estimate is robust and the temporal
variability in the H$\alpha$ that remains after we performed ISP subtraction
is due to changes in P Cyg’s intrinsic polarization; i.e. on average the
polarization in P Cyg’s H$\alpha$ line is a good estimator of the ISP, but
night-to-night there can be a variable, intrinsic polarization component to
the H$\alpha$ line. Therefore, we still report intrinsic $pfew$ H$\alpha$ line
values in Table 1 for completeness and because the remaining signal in the
lines likely holds key information about P Cyg’s winds.

Figure 6 shows the mean $P_{\%}$, $q$, $u$, and $\Theta$ spectra, taken at
each wavelength bin over time. Overall, these mean spectra are quite flat—as
is expected for polarization caused by free electron scattering. We also see
sharp features, most noteably at H$\alpha$. The magnitude of polarization is
indistinguishable from zero unless averaged over a wide wavelength bin, as is
the case for the $V$-band mean shown in Figure 5. This shows that—at least
when considered over a 13 year period—the ambient wind of P Cygni is nearly
symmetrical on the plane of the sky. Note that uncertainties become
significantly smaller at H$\alpha$, this is likely due to the increased number
of photons received at the emission line. The region around $H\alpha$ is shown
in more detail in Figure 8.

###  4.1 Intrinsic Polarization: the Search for Slope

Both N01 and T91 describe a decrease of polarization levels moving into the
infrared, which may be connected to the free-free absorptive
opacity—$\kappa_{ff}$—in the wind of P Cygni. This observed slope is used in
N01 to estimate the density of the inhomogeneities which are thought to give
rise to intrinsic polarization in P Cygni. In order to investigate this
previously discovered trend we fit a straight line to our data, this serves as
a proxy for the true relation for the attenuation of polarization caused by
$\kappa_{ff}$—$P_{\%}\propto e^{a\lambda^{2}}$. We perform this fit to our
data between 4000 and 7500 Å, after subtracting the ISP estimate found in §3,
using SciPy’s curve_fit function (Jones et al., 2001). This function utilizes
a least-squares fit of the form $y=mx+b$. One example of such a fit is shown
in the upper panel of Figure 7.

We report the results of this analysis in Figure 7 as a slope ($P_{\%}$Å-1)
for each observation with 1 $\sigma$ uncertainties. Because the exact value of
the slope for each individual observation is highly dependent on our ISP
calculation, it is most useful to look for trends across our observations. The
slope of $P_{\%}$ does appear to change between individual observations at a
statistically significant level. However, using Astropy’s (Robitaille et al.,
2013; Price-Whelan et al., 2018) Lomb-Scargle analysis software333All Lomb-
Scargle analysis was done using the default normalization as of Astropy v.
4.0.1 (Lomb, 1976; Scargle, 1982; Press & Rybicki, 1989) on these slopes
showed no statistically significant periodicity, which is consistent with the
stochastic variability displayed in Figure 7.

There do appear to be differences in the slopes of the $P_{\%}$ in our
observations taken with HPOL’s CCD (to the right of the vertical dashed line
in Figure 7) when compared to the Reticon detector (to the left of the
vertical dashed line). The CCD observations have a more consistently negative
slope, which corresponds to a persistent trend of a decrease in the $P_{\%}$
into the infrared consistent with N01’s findings. However, the observations
taken with HPOL’s Reticon detector have a larger overall scatter in their
slope, many of which are consistent with no wavelength dependence at the
3$\sigma$ uncertainty level. Some of our Reticon observations also show an
increasing wavelength dependence into the infrared (positive slope). It is
possible the overall polarization behavior of P Cygni changed when HPOL’s
detector was upgraded, no major changes in the system have been reported to be
observed in other data sets corresponding to this time frame. Still, it cannot
be ruled out that the difference in wavelength dependence between these two
detectors is due to instrumental differences, perhaps due to differences in
sensitivities between them.

###  4.2 H$\alpha$ Polarization

While it is true that intrinsic polarization at H$\alpha$ appears to be
suppressed, as we show in Figure 3, there appear to be more complicated
processes at work. Closer inspection of the H$\alpha$ line reveals a feature
in the polarization spectrum of P Cygni. As shown in Figure 8—though taking
the mean washes out some of this feature—and Figure 9, there is clear
structure at H$\alpha$ in both $P_{\%}$ and $\Theta$. While less easily seen
in the Reticon data, this feature has an amplitude larger than the mean level
of uncertainty of the CCD observations, suggesting that it is astrophysical in
nature. The typical uncertainty of the Reticon observations is large enough
that this feature is not clearly visible in many of those observations.

Figure 8 paints a puzzling picture. The strongest polarization appears near to
the peak of the H$\alpha$ emission line. Moving red-ward there is a strong and
rapid rotation, passing near to the origin in $q$ vs. $u$ space.

This feature shows temporal variations with no clear periodicity with regard
to amplitude or shape (Figure 9). There are no significant correlations
between the amplitude of the feature and the average percent polarization or
position angle across the line.

However, our Figure 7 shows possible night-to-night changes in the wavelength
dependence of P Cygni’s polarization, which was also noted by T91 and N01,
suggesting there are changes in P Cygni’s free-free absorptive wind opacity.
Because it is plausible that this absorption is also affecting the H$\alpha$
line and varying with time, we do not attempt to correct for it. Instead, we
attempt to quantify variability in the line—which may be due to changes in its
underlying absorption—using the methods described below.

###  4.3 Periodicity

Over the past century observers have noted irregular variability in P Cygni
(de Groot, 1969). As Kharadse (1936) put it, “P Cygni — Nova of 1600 — is one
of the most remarkable stars of the Northern sky. It has attracted the
attention of astronomers by an unusual character of light variation.” More
recently, however, evidence for periodicity in these unusual light variations
has been presented by various authors (de Groot et al., 2001a, b; van Genderen
A.~M., 2002; Michaelis et al., 2018; Kochiashvili et al., 2018). The periods
found range from $\sim 17$ day microvariations (de Groot et al., 2001b, a) to
$\sim 4.7$ years (Michaelis et al., 2018). As will be discussed in §5.1, there
are many fascinating phenomenae that can cause periodic signals in
polarimetric data. Here we detail our search for such signals.

Figure 10: Lomb-Scargle periodograms of $P_{\%}$ data. Frequencies picked by
the prewhitening procedure described in §4.3 are marked by vertical lines.
Fainter lines indicate locations of frequencies which are harmonics or linear
combinations of the darker, base frequencies (marked with $\dagger$ in Table
4.3). Errorbars at the top of each vertical line indicate error in the
frequency of the peak selected. Figure 11: Lomb-Scargle periodogram analyzing
simulated data which has the same temporal spacing of observations in Table 1,
based on a sinusoidal model with a period of 475 days. The horizontal dashed
line shows the $\frac{1}{1,000}$ false alarm level.

We use Astropy (Robitaille et al., 2013; Price-Whelan et al., 2018) to
calculate the Lomb-Scargle periodogram (Lomb, 1976; Scargle, 1982; Press &
Rybicki, 1989) for all available sets of $P_{\%}$ data from Table 1. Some sets
of $\Theta$ data were found to have periodicities, these periods were
significantly less robust under the prewhitening procedure described below.

We first search for peaks in the periodograms with low false alarm probability
(FAP) under the null hypothesis of white noise (Horne & Baliunas, 1986). This
results in multiple detected periods with FAP $\leq\frac{1}{1000}$. However,
not all of these periodicities may be the result of astrophysical phenomena.
As discussed extensively in VanderPlas (2018), there are many subtleties in
interpreting the Lomb-Scargle periodogram. Therefore, we also conduct a
prewhitening procedure, similar to that described in Dorn-Wallenstein et al.
(2019). In short, we iteratively select the highest peak in the periodogram,
fit a sinusoid to the data (allowing the frequency to vary within the
resolution of the periodogram), subtract the sinusoid, and recalculate the
periodogram until we reach a minimum in the Bayesian Information Content of
the fit (Schwarz, 1978). Formal errors on the frequency, amplitude, and phase
of the sinusoids extracted in each stage of prewhitening are calculated
following Lucy & Sweeney (1971) and Montgomery & Odonoghue (1999). This
results in a list of frequencies that uniquely describe the data to within the
noise. Frequencies extracted via prewhitening are listed—in the order that
they were detected by the prewhitening algorithm—in Table 4.3, and the full
periodograms are shown in Figure 10. We additionally search for harmonic
frequencies—$f_{i}=nf_{\circ}$ where n is an integer—and combination
frequencies of the form $f_{i}=f_{1}+f_{2}$.

These periodograms have relatively poor signal-to-noise ratios (SNRs), and
many of the extracted frequencies have large false alarm probabilities,
meaning that it is difficult to say whether any particular period is “real” or
simply a statistical artefact. However, several periods are more convincing
than others. For example, the 97 day period found in the $pfew$ continuum
$P_{\%}$ data has a small FAP (of order $10^{-4}$), is the base of multiple
harmonic/combination frequencies, and is also detected in the intrinsic
$V$-band $P_{\%}$ data.

Data |  Frequency [days-1] | Frequency Error |  Period [days] | SNR | False Alarm Probability
---|---|---|---|---|---
Observed $V$-Band | 0.001629† | 0.000031 | 614.0 | 5.0 | 0.048
0.016990† | 0.000038 | 58.9 | 5.6 | 0.129
0.030053† | 0.000055 | 33.3 | 3.7 |  $\sim$1
0.037419 | 0.000041 | 26.7 | 5.1 |  $\sim$1
0.036741 | 0.000061 | 27.2 | 4.9 |  $\sim$1
0.023797 | 0.000045 | 42.0 | 5.7 |  $\sim$1
0.004683† | 0.000045 | 213.5 | 5.3 |  $\sim$1
0.034663 | 0.000040 | 28.8 | 4.6 |  $\sim$1
0.006739† | 0.000067 | 148.4 | 4.9 | 0.874
0.031080† | 0.000055 | 32.2 | 9.4 |  $\sim$1
0.029281† | 0.000064 | 34.2 | 5.4 | 0.647
0.016665† | 0.000074 | 60.0 | 5.5 | 0.083
0.034099† | 0.000071 | 29.3 | 5.0 |  $\sim$1
Intrinsic $V$-Band | 0.034717 | 0.000035 | 28.8 | 4.6 | 0.001
0.026951 | 0.000032 | 37.1 | 7.0 | 0.023
0.018596 | 0.000036 | 53.8 | 4.1 | 0.049
0.010299† | 0.000038 | 97.1 | 6.9 | 0.206
0.020022 | 0.000049 | 49.9 | 4.7 | 0.927
0.031412 | 0.000042 | 31.8 | 5.1 | 0.117
0.035521 | 0.000047 | 28.2 | 6.0 | 0.025
0.032523 | 0.000069 | 30.7 | 5.5 |  $\sim$1
0.011359 | 0.000055 | 88.0 | 4.7 | 0.744
0.038363 | 0.000061 | 26.1 | 4.1 |  $\sim$1
0.030532 | 0.000058 | 32.8 | 7.4 |  $\sim$1
0.012354† | 0.000045 | 80.9 | 6.4 |  $\sim$1
0.004669† | 0.000039 | 214.2 | 4.9 |  $\sim$1
0.001084† | 0.000044 | 922.9 | 5.6 |  $\sim$1
0.016677† | 0.000059 | 60.0 | 5.0 | 0.811
0.018210 | 0.000061 | 54.9 | 6.1 | 0.455
0.014242 | 0.000066 | 70.2 | 5.5 | 0.152
0.023261 | 0.000058 | 43.0 | 4.1 | 0.948
0.003896† | 0.000200 | 256.6 | 5.2 |  $\sim$1
Observed $pfew$ Continuum | 0.021624† | 0.000037 | 46.2 | 5.1 | 0.018
0.034627† | 0.000038 | 28.9 | 5.4 | 0.280
0.027373† | 0.000038 | 36.5 | 6.4 |  $\sim$1
0.034168† | 0.000057 | 29.3 | 5.2 |  $\sim$1
0.006258† | 0.000044 | 159.8 | 6.3 |  $\sim$1
0.033576 | 0.000052 | 29.8 | 8.4 |  $\sim$1
0.029280† | 0.000083 | 34.2 | 4.6 |  $\sim$1
0.038049† | 0.000085 | 26.3 | 4.4 |  $\sim$1
0.026358† | 0.000053 | 37.9 | 6.0 |  $\sim$1
0.018158† | 0.000064 | 55.1 | 7.1 |  $\sim$1
Intrinsic $pfew$ Continuum | 0.010285† | 0.000035 | 97.2 | 5.0 | 0.0001
0.037901† | 0.000045 | 26.4 | 5.5 | 0.0007
0.020591 | 0.000033 | 48.6 | 5.0 | 0.012
0.034052† | 0.000053 | 29.4 | 4.3 | 0.591
Observed $pfew$ Line | 0.001637† | 0.000041 | 611.0 | 6.2 | 0.965
Intrinsic $pfew$ Line | 0.000720† | 0.000025 | 1388.2 | 3.3 | 0.118
0.038670† | 0.000039 | 25.9 | 7.7 | 0.162

Note. — † mark base frequencies—which are $not$ a harmonic or linear
combination of other frequencies.

We also note that several of the base periods (those marked with a $\dagger$
in Table 4.3) have multiple harmonics of the same order. This could be a
result of errors in either the base or harmonic frequencies. In fact, the
presence of any harmonics may be evidence that any true periodic oscillations
are not sinusoidal. If the polarization in P Cygni arises from the ejection of
asymmetries into the ambient stellar wind, as is described in T91 and N01, we
should not expect the time-series $P_{\%}$ curve to be a sinusoid. Under these
conditions that curve would peak sharply when such an asymmetry is ejected
from the photosphere.

####  4.3.1 Cadence Test

To determine if this observed periodicity is a result of the temporal spacing
of the observations, we perform a Lomb-Scargle periodicity search on simulated
data which lacks any of the periods that we have found. We claim that if we
detect a similar period in this test, then the periods found in Table 4.3 can
be attributed to how the data was sampled, not an astrophysical process.

To create a set of simulated data we used a series of periodic functions with
various periods (for example $\sin(\frac{2\pi t}{475})$), as well as flat
distributions. We calculated the values of these functions at the time values
of our observations (MJD) then added Gaussian noise defined by errors at each
observation in Table 1.

However, as is exemplified in Figure 11, all incorrect periods found in such
data sets are either clustered around the true period or a harmonic of the
true period. Based on this result, we conclude that the periodicities must be
either astrophysical in nature or a result of more complicated
instrumental/analytical effects. The likelihood of our findings being due to
HPOL is minimal as HPOL is a well characterized instrument and spurious
periodicities have not been found in other HPOL data sets. But, as can be seen
in Table 4.3, there are different periods found in intrinsic polarization as
compared to observed polarization. This could be a sign that the process of
subtracting out the interstellar polarization signal (see §3) may either
obscure or insert periodicities. For example, in $q$-$u$ space, the $P_{\%}$
of a circle centered at the origin would be constant, whereas one shifted away
from the origin would result in a sine wave pattern in $P_{\%}$.

Additionally, converting the MJD values in Table 1 to barycentric Julian date
had no effect on any of the results discussed here.

####  4.3.2 pfew Periodicity

By searching for periodicity in pfew H$\alpha$ line data, we can test whether
there is any non-stochastic variability present in P Cyg’s H$\alpha$ emission
line itself, as opposed to the V band continuum or the pfew continuum. Figure
10 shows our resulting Lomb-Scargle periodograms. The evidence for any
significant periods is much weaker than in the continuum data. This suggests
that any periodicities are not due to changes in the H$\alpha$ emission region
of P Cygni.

###  4.4 Ellipticity in $q$-$u$ Space: The Mauchly Test

The position angle of P Cygni’s intrinsic polarization appears to move
randomly, with no strong preference for one orientation. In other words,
whatever mechanism is responsible for P Cygni’s mass loss appears to eject
mass asymmetries stochastically in all directions. This can be seen from a
$q$-$u$ plot. A random distribution of position angles will produce a circular
distribution of points in q vs. u space, as appears to be the case in Figures
3 and 5.

Here, we attempt to quantify the statistical significance of ellipticity of
the V band q-u distribution using the Mauchly test (Mauchly, 1940a) to
determine if their position angle is truly random (i.e. circular in q-u space)
or if there is a slightly preferred position angle (i.e. elliptical in q-u
space). The Mauchly test has been used extensively in the field of
biostatistics (Davis, 2003; Myers et al., 2010; Crowder & Hand, 2017), but has
thus far not been used in the context of polarimetry.

In the original paper by Mauchly, he seeks to identify significant ellipticity
in a collection of two-dimensional data points. A collection of data points in
two dimensions may appear elliptically distributed even when their parent
population is actually circular, simply due to random sampling. Thus, Mauchly
asks “What is the probability of obtaining, from a circular normal population,
a random sample of N points for which the ellipticity is as great or greater
than that actually obtained in the given sample of N points?” (Mauchly,
1940b).

To help answer this question, Mauchly defines the ellipticity statistic:

$\mathcal{L}_{e}=\frac{2\sigma_{q}\sigma_{u}\sqrt{1-r^{2}}}{\sigma_{q}^{2}+\sigma_{u}^{2}},$
(5)

where $r$ is Pearson’s correlation coefficient between the distribution in $q$
and $u$, and where $\sigma^{2}_{q}$ and $\sigma^{2}_{u}$ are the sample
variances. By using $r$, $\sigma^{2}_{q}$, and $\sigma^{2}_{u}$, this test
assumes that the distribution can be described by a two dimensional Gaussian,
defined by the covariance matrix,

$K_{qu}=\begin{pmatrix}\sigma^{2}_{q}&r\\\ r&\sigma^{2}_{u}\end{pmatrix}.$

Note that if the distribution of points is circular, then
$\sigma^{2}_{q}=\sigma^{2}_{u}$, $r=0$, and $\mathcal{L}_{e}=1$. Conversely,
if the distribution of points is elliptical, then
$\sigma^{2}_{q}\neq\sigma^{2}_{u}$, $r\rightarrow 1$, and
$\mathcal{L}_{e}\rightarrow 0$.

Assuming that $N$ points are drawn from a circular distribution, the statistic
$\mathcal{L}_{e}$ has the distribution

$f(\mathcal{L}_{e})=(N-2)\mathcal{L}_{e}^{N-3}d\mathcal{L}_{e}.$ (6)

If the actual value of $\mathcal{L}_{e}$ calculated from the data is
$\mathcal{L}_{e}$, then the probability that a value as small or smaller than
$\mathcal{L}_{e}$ is found from a random sample of $N$ points is

$P(\mathcal{L}_{e})=\int_{0}^{\mathcal{L}_{e}}(N-2)\mathcal{L}_{e}^{N-3}d\mathcal{L}_{e}=\mathcal{L}_{e}^{N-2}.$
(7)

When $P(\mathcal{L}_{e})$ is close to one, then the null hypothesis cannot be
rejected (i.e. the parent distribution is likely circular) (Mauchly, 1940b).
We take the position that a value of $P(\mathcal{L}_{e})$ which is less than
$\alpha=0.05$ provides sufficient evidence to reject the null hypothesis.

To illustrate how these numbers are interpreted, Table 2 includes three
general examples in the last three rows. The first is a distribution which is
very close to being circular, however, the small amount of ellipticity present
is statistically significant. The next is an extremely elliptical distribution
which would appear to be a line. But, there is a 50% chance that points drawn
from a circular distribution would appear more elliptical. The last example is
an elliptical distribution which is very statistically significant; this
example would be similar to what we would see for a significant preferred
position angle.

Figure 12: White crosses show the position of observed $V$-band observations
of P Cygni in $q$ vs. $u$ space, with the size of the cross indicating
uncertainty. Underlying colormap shows the 2D Gaussian distribution derived
from $r$, $\sigma_{q}$, and $\sigma_{u}$ (discussed in §4.4). Contours show
equally spaced slices in likelihood, the lightest colors corresponding to a
peak in the probability distribution function.

Using all $q$-$u$ data sets from Table 1 we apply this Mauchly test. For the
observed $V$-band values, plotted in Figure 12, the sample standard deviations
and Pearson correlation coefficient are
$\sigma_{q}=0.202;~{}\sigma_{u}=0.177;~{}r=0.115$. Given these values we find
that $\mathcal{L}_{e}=0.985$ and $P(\mathcal{L}_{e})=0.320$. This is not
significant evidence for a preferred position angle and is consistent with the
findings of T91 and N01. In fact, analyzing any other data set in Table 1 in
this manner yields similar results. The fact that the observed $pfew$-line
observations do not show evidence for a preferred position angle actually
strengthens the argument for using it as an ISP probe, since deviations will
then be evenly distributed about the desired value.

####  4.4.1 A Test of the Mauchly Test

To see if the Mauchly Test is capable of finding a preferred position angle we
use the method described above on a system with a known preferred position
angle. Here we use HPOL data from Hoffman et al. (1998); Lomax et al. (2012)
of $\beta$ Lyrae—a self-eclipsing binary with mass transfer via an accretion
disk. The presence of this disk creates net linear polarization with a
preferred position angle perpendicular to the disk.

Using 69 non-ISP removed, $V$-band polarimetric observations from Lomax et al.
(2012), we find the distribution shown in the left panel of Figure 13, with
$\sigma_{q}=0.089;~{}\sigma_{U}=0.099;~{}r=-0.203$. This corresponds to
$\mathcal{L}_{e}=0.975$ and $P(\mathcal{L}_{e})=0.178$. This would seem to
suggest that the Mauchly test would not be capable of detecting the present
preferred position angle.

However, there is obviously a set of outliers towards the bottom of the
figure. These points are identified as outliers by Lomax et al. (2012), and
all three observations were taken on the same night. When these outliers are
discounted, the model of this distribution changes to that of the right panel
of Figure 13. Using this subset of data, we derive the values
$\sigma_{q}=0.089;~{}\sigma_{u}=0.074;~{}r=-0.443$. This corresponds to
$\mathcal{L}_{e}=0.880$ and $P(\mathcal{L}_{e})=0.0003$. which is
significantly below our threshold of $\alpha=0.05$. Thus, we can reject the
null hypothesis that the distribution is circular.

This example outlines an important limitation of this test. In small data-
sets, the presence of even a small number of outliers can have an impact on
the estimated value of $\mathcal{L}_{e}$ and the results of the associated
significance test.

Figure 13: Stokes $q$-$u$ data for a test system with a preferred position angle, $\beta$ Lyr from Lomax et al. (2012), see §4.4.1. White crosses represent individual $V$-band observations, with the size of the cross indicating uncertainty, while the underlying colormap shows the distribution derived from $r$, $\sigma_{q}$, and $\sigma_{u}$ (discussed in §4.4). Contours of the color map are equally spaced slices in likelihood, the lightest colors corresponding to a peak in the probability distribution function. The left panel shows the distribution with the three outliers included, the right panel shows the distribution without those observations. Table 2: Ellipticity  Data | $\mathcal{L}_{e}$ | $P(\mathcal{L}_{e})$
---|---|---
P Cygni | 0.986 | 0.354
$\beta$ Lyrae | 0.975 | 0.178
$\beta$ Lyrae (adjusted) | 0.880 | 0.0003
Significant Circle | 0.95 | 1e-05
Insignificant Line | 0.01 | 0.5
Significant Ellipse | 0.5 | 1e-05

##  5 Discussion

These results allow us to infer properties of the material in the near
circumstellar environment around P Cygni. Reporting on 13 years of
observations of P Cygni we can revisit the findings of previous
spectropolarimetric studies.

This paper supports previous claims from T91, N01, and Hayes (1985) that
asymmetries in the system lack a preferred polarimetric position angle. Not
only is the distribution of observations in $q$-$u$ space nearly circular, but
the small amount of observed ellipticity is statistically insignificant as
shown in Table 2. Previous studies investigating the resolved nebula around P
Cygni (Leitherer & Zickgraf, 1987; Barlow et al., 1994; Nota et al., 1995;
Meaburn et al., 1996) have reported that this nebula is roughly spherical in
shape, with some clumping. This is not entirely expected as more than 50% of
resolved galactic LBVs have bipolar nebulae (Weis, 2012). Our findings suggest
that this lack of observed directionality in the circumstellar environment
starts at the base of the wind. This means that this mode of mass loss is
either truly lacking any directional preference, or that the preferred
direction is hidden by projection effects. It is possible, for example, that
the mass inhomogeneities could be ejected solely from the equator of the star,
but if the equatorial plane is close to the plane of the sky then it would be
nearly impossible to tell.

It is more difficult to interpret the wavelength dependent feature observed at
H$\alpha$. This is troubling, particularly in light of the importance of using
the line center method for ISP estimates of early-type stars. The higher
levels of polarization in the line (shown in Figure 8) could potentially skew
measurements of the ISP in studies of P Cygni and other strong emission-line
objects. While this feature has a relatively small amplitude, and it likely
does not have a significant impact on the analysis in this work, we cannot be
certain as to the extent of its impact until it is explained.

###  5.1 Problematic Polarimetric Periodicities

The series of periodicities in the polarimetric data poses a serious question.
What could cause these periodicities? While it is unlikely that all
instrumental effects can be ruled out, HPOL was a well-characterized
instrument. Therefore, for the remainder of this section we shall attempt to
explain this result as a truly astrophysical phenomenon.

We did not find strong evidence for periodicity in the H$\alpha$ emission line
(using $pfew$ as described in §2.1.2) or in the slopes of the individual
observations’ $P_{\%}$ versus wavelength data. This suggests that the cause of
these periods is likely not due to changes in the free-free absorptive opacity
of P Cygni’s wind, which is thought to cause a general decrease in its
$P_{\%}$ into the infrared.

These results could be explored with great depth, and likely deserves their
own paper, but here we will discuss the major points, and possible
explanations for the periodic effects observed in these data. In particular,
we focus the majority of this discussion on the 97 day period in the continuum
polarization found in §4.3, as that is the period with the strongest false
alarm probability and has proven quite robust over different iterations of
analysis.

####  5.1.1 P Cygni Absorption

It is interesting to note that there are periodicities which appear in the
pfew Continuum $P_{\%}$ that are not present in the $V$-Band $P_{\%}$ data.
Given how close these two regions of continuum are, this is somewhat
unexpected. If these differences are astrophysical in nature it points to a
difference not in the H$\alpha$ emission line, but to a difference in what is
happening at H$\alpha$ in continuum emission.

A possible explanation could be the eponymous P Cygni absorption component of
the H$\alpha$ line. Since our implementation of the pfew method does not
account for this absorption component, it should have some impact on the
continuum values tabulated in Table 1. Therefore, it is possible that we are
detecting periodicity in the polarization of the absorption component itself.

####  5.1.2 Pulsations

The 97 day period listed in Table 4.3 lines up with $\sim 100$-d type periods
found by de Groot et al. (2001a) in a study of the photometric behavior of P
Cygni over multiple decades. The authors attempt to explain this mode of
variability as being a result of radial pulsations. They find an anti-
correlation between $m_{V}$ and $B-V$; the star becomes redder as it gets
brighter, as would be expected for radial pulsations where the increased
brightness is the result of an increase in radius and subsequent drop in
effective temperature. While it is not likely that radial pulsations would
directly affect the net polarization that we observe, inhomogeneities in
pulsation-driven mass loss from the star (expected when the star is in its
cool state, with a more extended atmosphere and lower surface gravity) could
lead to periodic variations in polarization.

Hydrodynamic models of LBV pulsations by Jiang et al. (2018) have also
predicted pulsations with periods of order $\sim$day timescales. While it is
less likely that the data presented here would be sensitive to such short
periods, these fast pulsations may combine with other effects to produce the
periodicity which we see.

####  5.1.3 Rotation

It is possible that inhomogeneities are being sourced from a stationary
feature on the photosphere, such as a hotspot, corotating interaction region,
or some other localized site of asymmetric mass loss. Using a radius of
76$R_{\sun}$ and $v\sin{i}=35$ kms-1 (Najarro et al., 1997), the rotational
period of P Cygni would be roughly 110 days. Since free-electron scattering
produces maximum polarization for a 90∘ scattering angle, we would observe the
effect of asymmetries most strongly when this feature is at the plane of the
sky.

However, we would expect to see peaks in polarization twice per rotation
period for this scenario to make sense in the context of our data set—which is
a higher frequency than most of the base frequencies we have found.
Additionally, it is not certain that a feature such as a hotspot or mass loss
asymmetry can be stable over a 13 year period.

####  5.1.4 A Companion

The potential binarity of LBVs is a contentious issue, but one explanation for
the period found in this data is that P Cygni has a binary companion. Periodic
variability in polarimetric data arising from companions is not new. For
example, Wolf-Rayet-O-star and Roche-lobe overflow binary systems show
variability in this manner (see Brown et al. 1978; Moffat & Piirola 1993;
Hoffman et al. 1998; Lomax et al. 2015).

Assuming the masses of P Cygni and its possible companion are each between 20
and 60 $M_{\sun}$, we can place some rough limits on the geometry of the
system. For these masses, and a 97 day orbital period, we find a semi-major
axis of roughly 1-2 AU, or about 3-6 times the radius of P Cygni (Najarro,
2001).

However, the observations in $q$-$u$ space do not behave like typical binary
systems. Figure 14 shows the somewhat erratic path that H$\alpha$ continuum
observations take as a function of phase when phase folded by a 97 day period,
which shows the most orderly path of all periods in Table 4.3. Typically,
polarization in binary systems creates elliptical patterns in q-u space (Brown
et al., 1978). The path in Fig. 14 is not obviously elliptical; this may be
due to the irregular cadence and very long-term data set compared with this
short period.

While this does not prove that P Cygni is not a binary, it does cast doubt on
a simple binary model of the P Cygni system.

Figure 14: pfew Continuum $q$-$u$ observations, phase folded to a 97 day
period, and re-binned to 10 phase bins. Solid line shows path between adjacent
bins in phase-space.

###  5.2 The Mauchly Test in Polarimetry

The shape of distributions in $q$-$u$ space encodes important information
about the geometry of objects with time varying polarization. For example, N01
discusses how polarimetric observations of Be stars are generally constrained
to a line in $q$-$u$ space. The Mauchly Test would find these systems quickly
and computationally efficiently.

It is important to note two strengths of the Mauchly test:

  * •

This test works even in the case of unresolved features.

  * •

Since this technique only measures the shape of a distribution, as long as the
data being analyzed is at one wavelength, ISP-correction has no effect.

Additionally, for objects that are already known to have a preferred position
angle, this statistic can quantify how “preferred” it really is. In the case
of a star with a circumstellar disk, ellipticity of polarimetric observations
could yield the inclination angle of the disk with respect to the line of
sight. Here, $\mathcal{L}_{e}=0$ corresponds to viewing the disk edge on and
$\mathcal{L}_{e}=1$ corresponds to a face on alignment.

The eagle-eyed readers may ask, “But what about $\beta$ Lyr? That is an
eclipsing binary system (Hoffman et al., 1998; Lomax et al., 2012), but your
analysis found $\mathcal{L}_{e}=0.88$. Are you claiming that we are somehow
seeing the $\beta$ Lyr system face on?”. Thankfully we are not making that
claim. There is some nuance in how we must interpret this statistical tool.
$\beta$ Lyr, for which we found $\mathcal{L}_{e}=0.88$ and
$P(\mathcal{L}_{e})=0.0003$, is an excellent example. These values are
measurements of our data, not of the parent distribution. What these values
tell us is that it is very unlikely the the ellipticity of the parent
distribution is greater than $0.88$. We have placed an upper bound on the true
value of $\mathcal{L}_{e}$.

##  6 Conclusion

Using 76 observations taken over 13 years we have been able to update and
confirm the conclusions of past work on the circumstellar environment around P
Cygni.

  * •

P Cygni produces intrinsic polarization, indicating time varying asymmetries
in the geometry of the circumstellar material.

  * •

Using the line center method at H$\alpha$ and by fitting the Serkowski law we
find a $\lambda_{max}$ of 4595 Å and a $P_{max}$ of 1.23% for the ISP. This
implies that in the direction of P Cygni $R_{V}=2.638\pm 0.028$.

  * •

We utilize the Mauchly ellipticity statistic to find that the ellipticity of
the distribution of observations in $q$-$u$ space is small and statistically
insignificant. Thus we do not see evidence for a preferred position angle.

  * •

H$\alpha$, while being less polarized than the continuum on average, contains
a wavelength dependent polarization feature. Given its appearance across
multiple HPOL observations, and in Davies et al. (2005), this feature is most
likely astrophysical in nature.

  * •

There are a set of statistically significant periods in various aspects of our
data. These periodicities are in the $V$-Band, and H$\alpha$ Continuum data,
but not in the H$\alpha$ emission line polarization.

The feature found in H$\alpha$ is particularly intriguing, for many reasons.
It provides a cautionary tale for those wishing to use the line-center method.
While we do not believe this structure has had a significant impact on the ISP
estimate derived in §3, it exemplifies how the base assumptions of the line-
center method—that strong emission lines in hot stars are intrinsically
unpolarized—can break down. This strange feature is also a reminder of how
much is left to discover, even about an object like P Cygni which has been
studied for more than 400 years, and even about its H$\alpha$ line which has
become the prototype for an entire class of similar spectral features. In this
sense, this study joins a long and distinguished astronomical tradition, as
Wilson (1936) wrote after finding variability in multiple absorption features
in the spectra of P Cygni taken over the course of 20 years at Lick
Observatory, “…although these variations have hitherto escaped observation,
they are by no means minor, but constitute one of the major features of the
problem.”

More complete and detailed observations of P Cygni will lead to a better
understanding of its spectropolarimetric behaviour. The current available data
on P Cygni shows its long term behavior. However, it would be incredibly
valuable to obtain higher cadence data, particularly, over the course of the
longer periods discussed in §4.3. Such a campaign could reveal potential
correlations between the different aspects discussed here and provide a more
unified view of the mass-loss of P Cygni. Additionally, observing the ejection
of an inhomogeneity, and the $\sim$14 day period afterwards—which Hayes (1985)
called the characteristic timescale of polarimetric changes in P Cygni—could
reveal much about the nature of the mass-loss in this system.

We gratefully acknowledge the contributions of the large team of HPOL and
WUPPE observers from the University of Wisconsin who helped to acquire and
reduce these data. We also thank Inger Olovsson, curator at the Collections of
Skokloster Castle, for her timely and enthusiastic assistance in locating and
confirming the existence of Blaeu’s globe chronicling the discovery of P
Cygni. This research has made use of the SIMBAD database, operated at CDS,
Strasbourg, France. Some of the data presented in this paper were obtained
from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by
the Association of Universities for Research in Astronomy, Inc., under NASA
contract no. NAS5-26555. Support for MAST for non-HST data is provided by the
NASA Office of Space Science via grant no. NNX09AF08G and by other grants and
contracts. This work was funded by Washington NASA Space Grant Consortium,
NASA Grant #NNX15AJ98H and a Cottrell Scholar Award to EML by the Research
Corporation for Scientific Advancement. The authors acknowledge that the work
presented was largely conducted on the traditional land of the first people of
Seattle, the Duwamish People past and present, and honor with gratitude the
land itself and the Duwamish Tribe.

## References

  * Anderson et al. (1996) Anderson, C. M., Weitenbeck, A. J., Code, A. D., et al. 1996, The Astronomical Journal, 112, 2726, doi: 10.1086/118217
  * Avcioğlu (1984) Avcioğlu, K. 1984, Astrophysics and Space Science, 106, 341, doi: 10.1007/BF00650359
  * Barlow et al. (1994) Barlow, M. J., Drew, J. E., Meaburn, J., & Massey, R. M. 1994, Monthly Notices of the Royal Astronomical Society, 268, L29
  * Bessell (1990) Bessell, M. S. 1990, Publications of the Astronomical Society of the Pacific, 102, 1181, doi: 10.1086/132749
  * Bjorkman et al. (1993) Bjorkman, K. S., Meade, M. R., Nordsieck, K. H., et al. 1993, The Astrophysical Journal, 412, 810, doi: 10.1086/172963
  * Blaeu (1602) Blaeu, W. J. 1602, Skokloster Castle
  * Brown et al. (1978) Brown, J. C., McLean, I. S., & Emslie, A. G. 1978, Astronomy and Astrophysics, 68, 415
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, The Astrophysical Journal, 345, 245, doi: 10.1086/167900
  * Coyne & Gehrels (1967) Coyne, G. V., & Gehrels, T. 1967, Astronomical Journal, 72, 887
  * Coyne et al. (1974) Coyne, G. V., Gehrels, T., & Serkowski, K. 1974, Astronomical Journal, 79, 581
  * Crowder & Hand (2017) Crowder, M. J., & Hand, D. J. 2017, Analysis of Repeated Measures (Routledge), doi: 10.1201/9781315137421
  * Davidson et al. (2014) Davidson, J. W., Bjorkman, K. S., Hoffman, J. L., et al. 2014, Journal of Astronomical Instrumentation, 03, 1450009, doi: 10.1142/S2251171714500093
  * Davies et al. (2005) Davies, B., Oudmaijer, R. D., & Vink, J. S. 2005, Astronomy & Astrophysics, 439, 1107, doi: 10.1051/0004-6361:20052781
  * Davies et al. (2006) Davies, B., Oudmaijer, R. D., & Vink, J. S. 2006, in Stars with the B[e] Phenomenon. ASP Conference Series., ed. M. Kraus & A. S. Miroshnichenko, 173
  * Davis (2003) Davis, C. S. 2003, Statistical Methods for the Analysis of Repeated Measurements, Springer Texts in Statistics (New York, NY: Springer New York), doi: 10.1007/b97287.  http://link.springer.com/10.1007/b97287
  * de Groot & Mart (1988) de Groot, & Mart. 1988, Irish Astronomical Journal, 18, 163
  * de Groot (1969) de Groot, M. 1969, Bulletin of the Astronomical Institutes of the Netherlands, 20, 225
  * de Groot et al. (2001a) de Groot, M., Sterken, C., & van Genderen, A. M. 2001a, Astronomy & Astrophysics, 376, 224, doi: 10.1051/0004-6361:20010960
  * de Groot et al. (2001b) de Groot, M., Sterken, C., & Van Genderen, A. M. 2001b, Astronomy and Astrophysics, 376, 224, doi: 10.1051/0004-6361:20010960
  * Dorn-Wallenstein et al. (2019) Dorn-Wallenstein, T. Z., Levesque, E. M., & Davenport, J. R. A. 2019, The Astrophysical Journal, 878, 155, doi: 10.3847/1538-4357/ab223f
  * Draine (2003) Draine, B. T. 2003, The Astrophysical Journal, 598, 1017, doi: 10.1086/379118
  * Groh et al. (2014) Groh, J. H., Meynet, G., Ekström, S., & Georgy, C. 2014, Astronomy & Astrophysics, 564, A30, doi: 10.1051/0004-6361/201322573
  * Hayes (1985) Hayes, D. P. 1985, Astrophysical Journal, 84, 487
  * Heiles (1999) Heiles, C. 1999, The Astronomical Journal, 119, 923, doi: 10.1086/301236
  * Hoffman et al. (1998) Hoffman, J. L., Nordsieck, K. H., & Fox, G. K. 1998, The Astronomical Journal, 115, 1576, doi: 10.1086/300274
  * Horne & Baliunas (1986) Horne, J. H., & Baliunas, S. L. 1986, The Astrophysical Journal, 302, 757, doi: 10.1086/164037
  * Humphreys & Davidson (1979) Humphreys, R. M., & Davidson, K. 1979, The Astrophysical Journal, 232, 409
  * Humphreys & Davidson (1994) —. 1994, Publications of the Astronomical Society of the Pacific, 106
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Israelian & de Groot (1999) Israelian, G., & de Groot, M. 1999, Space Science Reviews, 90, 493, doi: 10.1023/A:1005223314464
  * Jiang et al. (2018) Jiang, Y. F., Cantiello, M., Bildsten, L., et al. 2018, Outbursts of luminous blue variable stars from variations in the helium opacity, Nature Publishing Group, doi: 10.1038/s41586-018-0525-0
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., & and others. 2001, SciPy: Open source scientific tools for Python.  http://www.scipy.org/
  * Kharadse (1936) Kharadse, E. K. 1936, \zap, 11, 304
  * Kochiashvili et al. (2018) Kochiashvili, N., Beradze, S., Natsvlishvili, R., et al. 2018, Astrophysics, 61, 22, doi: 10.1007/s10511-018-9512-2
  * Leitherer & Zickgraf (1987) Leitherer, C., & Zickgraf, F. 1987, Astronomy and Astrophysics, 174, 103
  * Lomax et al. (2012) Lomax, J. R., Hoffman, J. L., Elias II, N. M., Bastien, F. A., & Holenstein, B. D. 2012, The Astrophysical Journal, 750, 59, doi: 10.1088/0004-637X/750/1/59
  * Lomax et al. (2015) Lomax, J. R., Nazé, Y., Hoffman, J. L., et al. 2015, Astronomy & Astrophysics, 573, A43, doi: 10.1051/0004-6361/201424468
  * Lomb (1976) Lomb, N. R. 1976, Astrophysics and Space Science, 39, 447, doi: 10.1007/BF00648343
  * Lucy & Sweeney (1971) Lucy, L. B., & Sweeney, M. A. 1971, \aj, 76, 544, doi: 10.1086/111159
  * Mauchly (1940a) Mauchly, J. W. 1940a, The Annals of Mathematical Statistics, 11, 204, doi: 10.1214/aoms/1177731915
  * Mauchly (1940b) —. 1940b, Journal of Geophysical Research, 45, 145, doi: 10.1029/TE045i002p00145
  * Meaburn et al. (1996) Meaburn, J., Lopez, J. A., Barlow, M. J., & Drew, J. E. 1996, Monthly Notices of the Royal Astronomical Society, 283, L69, doi: 10.1093/mnras/283.3.L69
  * Michaelis et al. (2018) Michaelis, A. M., Kashi, A., & Kochiashvili, N. 2018, \na, 65, 29, doi: 10.1016/j.newast.2018.06.001
  * Moffat & Piirola (1993) Moffat, A. F. J., & Piirola, V. 1993, The Astrophysical Journal, 413, 724, doi: 10.1086/173040
  * Montgomery & Odonoghue (1999) Montgomery, M. H., & Odonoghue, D. 1999, Delta Scuti Star Newsletter, 13, 28
  * Myers et al. (2010) Myers, J. L., Well, A. D., & Lorch, R. F. J. 2010, Research Design and Statistical Analysis, 3rd edn. (Routledge)
  * Najarro (2001) Najarro, F. 2001, in P Cygni 2000: 400 Years of Progress, ASP Conference Proceeding (Astronomical Society of the Pacific), 294.  https://ui.adsabs.harvard.edu/abs/2001ASPC..233..133N/abstract
  * Najarro et al. (1997) Najarro, F., Hillier, D. J., & Stahl, O. 1997, Astronomy & Astrophysics, 326, 1117
  * Nook et al. (1990) Nook, M. A., Cardelli, J. A., & Nordsieck, K. H. 1990, The Astronomical Journal, 100, 2004, doi: 10.1086/115655
  * Nordsieck et al. (1994) Nordsieck, K. H., Code, A. D., Anderson, C. M., et al. 1994, in X-Ray and Ultraviolet Polarimetry, ed. S. Fineschi, Vol. 2010, 2–11.  http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=935224
  * Nordsieck & Harris (1996) Nordsieck, K. H., & Harris, W. 1996, in Polarimetry of the Interstellar Medium, ed. W. G. Roberge & D. C. B. Whittet.  http://articles.adsabs.harvard.edu/pdf/1996ASPC...97..100N
  * Nordsieck et al. (2001) Nordsieck, K. H., Wisniewski, J., Babler, B. L., et al. 2001, P Cygni 2000: 400 Years of Progress, ASP Conference Proceeding, 233, 261
  * Nota et al. (1995) Nota, A., Livio, M., Clampin, M., & Schulte-Ladbeck, R. 1995, The Astrophysical Journal, 448, 788, doi: 10.1086/176006
  * Press & Rybicki (1989) Press, W. H., & Rybicki, G. B. 1989, The Astrophysical Journal, 338, 277, doi: 10.1086/167197
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, The Astronomical Journal, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Richardson & Mehner (2018) Richardson, N. D., & Mehner, A. 2018, Research Notes of the AAS, 2, 121, doi: 10.3847/2515-5172/aad1f3
  * Robitaille et al. (2013) Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, Astronomy & Astrophysics, 558, A33, doi: 10.1051/0004-6361/201322068
  * Scargle (1982) Scargle, J. D. 1982, The Astrophysical Journal, 263, 835, doi: 10.1086/160554
  * Schwarz (1978) Schwarz, G. 1978, Annals of Statistics, 6, 461, doi: 10.1214/AOS/1176344136
  * Serkowski et al. (1975) Serkowski, K., Mathewson, D. L., & Ford, V. L. 1975, The Astrophysical Journal, 196, 261, doi: 10.1086/153410
  * Smith et al. (2003) Smith, N., Gehrz, R. D., Hinz, P. M., et al. 2003, The Astronomical Journal, 125, 1458
  * Smith & Hartigan (2006) Smith, N., & Hartigan, P. 2006, Astrophysical Journal, 638, 1045
  * Taylor et al. (1991a) Taylor, M., Nordsieck, K. H., Schulte-Ladbeck, R. E., & Bjorkman, K. S. 1991a, The Astronomical Journal, 102, 1197, doi: 10.1086/115947
  * Taylor et al. (1991b) Taylor, M., Code, A. D., Nordsieck, K. H., et al. 1991b, The Astrophysical Journal, 382, L85, doi: 10.1086/186218
  * Turner et al. (2001) Turner, D. G., Welch, G., Graham, M., et al. 2001, Journal of the American Association of Variable Star Observers (JAAVSO), 29, 73
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * van Genderen A.~M. (2002) van Genderen A.~M., S. C. D. G. M. 2002, Journal of Astronomical Data, 8, 8
  * VanderPlas (2018) VanderPlas, J. T. 2018, The Astrophysical Journal Supplement Series, 236, 16, doi: 10.3847/1538-4365/aab766
  * Weis (2012) Weis, K. 2012, in Proceedings of a Scientific Meeting in Honor of Anthony F. J. Moffat, ed. L. Drissen, C. Robert, N. St-Louis, & A. F. J. Moffat, Vol. 465, Auberge du Lac Taureau, St-Michel-Des-Saints, Qu{\’e}bec, Canada, 213.  https://ui.adsabs.harvard.edu/abs/2012ASPC..465..213W/abstract
  * Whittet et al. (1992) Whittet, D. C. B., Martin, P. G., Hough, J. H., et al. 1992, The Astrophysical Journal, 386, 562, doi: 10.1086/171039
  * Wilking et al. (1982) Wilking, B. A., Lebofsky, M. J., & Rieke, G. H. 1982, The Astronomical Journal, 87, 695, doi: 10.1086/113147
  * Wilson (1936) Wilson, O. C. 1936, \apj, 84, 296, doi: 10.1086/143765
  * Wisniewski et al. (2010) Wisniewski, J. P., Draper, Z. H., Bjorkman, K. S., et al. 2010, Astrophysical Journal, 709, 1306, doi: 10.1088/0004-637X/709/2/1306
  * Wolff et al. (1996) Wolff, M. J., Nordsieck, K. H., & Nook, M. A. 1996, The Astronomical Journal, 111, 856, doi: 10.1086/117833
  * Wood et al. (1996a) Wood, K., Bjorkman, J. E., Whitney, B. A., & Code, A. D. 1996a, The Astrophysical Journal, 461, 828, doi: 10.1086/177105
  * Wood et al. (1996b) —. 1996b, The Astrophysical Journal, 461, 828, doi: 10.1086/177105