H\alpha Variability of P Cygni

# The Hα Variations of the Luminous Blue Variable P Cygni: Discrete Absorption Components and the Short S Doradus Phase

N. D. Richardson11affiliation: Center for High Angular Resolution Astronomy, Department of Physics and Astronomy, Georgia State University, P. O. Box 4106, Atlanta, GA 30302-4106; richardson@chara.gsu.edu, gies@chara.gsu.edu , N. D. Morrison22affiliation: Ritter Astrophysical Research Center, Department of Physics and Astronomy, University of Toledo, 2801 W. Bancroft, Toledo, OH 43606; nmorris@utnet.utoledo.edu , D. R. Gies11affiliation: Center for High Angular Resolution Astronomy, Department of Physics and Astronomy, Georgia State University, P. O. Box 4106, Atlanta, GA 30302-4106; richardson@chara.gsu.edu, gies@chara.gsu.edu , N. Markova33affiliation: Institute of Astronomy with National Astronomical Observatory, Bulgarian Academy of Sciences, P.O. Box 136, 4700 Smoljan, Bulgaria; nmarkova@astro.bas.bg , E. N. Hesselbach44affiliation: University of Notre Dame, Department of Physics, 225 Nieuwland Science Hall, Notre Dame, IN 46556; ehesselb@nd.edu , and J. R. Percy55affiliation: Department of Astronomy and Astrophysics, University of Toronto, Toronto, ON, Canada, M5S 3H4; jpercy@utm.utoronto.ca
###### Abstract

P Cygni is a prototype of the Luminous Blue Variables (or S Doradus variables), and the star displays photometric and emission line variability on a timescale of years (known as the “short S Doradus phase” variations). Here we present new high resolution H spectroscopy of P Cyg that we combine with earlier spectra and concurrent -band photometry to document the emission and continuum flux variations over a 24 y time span. We show that the emission and continuum fluxes vary in concert on timescales of 1.6 y and longer, but differ on shorter timescales. The H profile shape also varies on the photometric timescales, and we describe the observed co-variations of the emission peak and absorption trough properties. We argue that the episodes of photometric and emission brightening are caused by increases in the size of the emission region that are related to variations in wind mass loss rate and outflow speed. We find evidence of blueward accelerating, Discrete Absorption Components (DACs) in the absorption trough of the H profile, and these features have slower accelerations and longer durations than those observed in other lines. The DAC strengths also appear to vary on the photometric timescales, and we suggest that the propagation of the DAC-related wind structures is closely related to changes in the overall wind mass loss rate and velocity.

stars: variables — stars: early-type — stars: individual (P Cyg, HD 193237) — stars: winds, outflows — stars: circumstellar matter — stars: mass loss
slugcomment: Accepted by AJ

## 1 Introduction

Luminous Blue Variables (LBVs or S Doradus variables) are evolved, massive stars. LBVs are characterized by large mass loss rates and variability on multiple timescales. The two “prototypical” Galactic LBVs are Carinae and P Cygni, and they probably represent different extremes of mass loss rate within the scheme of LBV evolution (Israelian & de Groot, 1999). One of the defining criteria of the LBVs is the observation of a large scale eruption, when the star brightens by several magnitudes. The quiescent times between these eruptions may last centuries. In addition to such rare, giant eruptions, these stars also display lesser photometric and spectroscopic variations on other timescales (e.g.  Humphreys & Davidson 1994). van Genderen (2001) defines the S Doradus (SD-) phase to be the moderate, long-term, brightening and fading phases. There are two types of these phases, short and long, with similar characteristics. The short SD-phase is typically on the order of years ( 10 years), while the longer SD-phase is on a timescale of decades. These phases are thought to originate from changes in the star’s photosphere, and both may have the same physical driving mechanism. These long-term variations are observed to differ from cycle to cycle, both in duration and amplitude (Sterken et al., 1997).

P Cygni (HD 193237, HR 7763, Nova Cyg 1600) remains one of the most fascinating objects in the sky. It was discovered during its first recorded great eruption in 1600 by Willem Janszoon Blaeu, a Dutch chart-maker and mathematician. During this eruption, the star brightened to about 3rd magnitude for about six years, and then it faded from visibility by 1626. It rose again in 1654 to about the same maximal brightness, where it remained for five years. The star faded after this, and although its long-term variability is poorly documented, the star has been slowly brightening to its current magnitude of about 4.8 (Israelian & de Groot, 1999). The slow brightening may reflect evolutionary changes (de Groot & Lamers 1992; Lamers & de Groot 1992; Langer et al. 1994).

Long-term photometric monitoring of P Cygni began in the 1980s when Percy & Welch (1983), Percy et al. (1988), and de Groot (1990) embarked on extended observing campaigns. Their observations showed that the variations often occur on three characteristic timescales: a short 17 day variation similar to the Cygni type variations observed in hot supergiants, a 100 day “quasi-period” similar to that observed in other LBVs, and a long-term cycle (years) attributed to a short-SD phase (de Groot et al. 2001; Percy et al. 2001).

According to Israelian & de Groot (1999), comprehensive spectral monitoring of P Cygni was started by Luud (1967) and Markova (1993), among others. The first long-term spectroscopic monitoring campaign of P Cygni was presented in seminal papers by Markova et al. (2001a,b). They found evidence of co-variability of the H emission line strength and Johnson photometry, indicating a short-SD phase with a quasi-period of 7 years, although their observations did not fully cover two cycles. The variations were attributed to inversely correlated changes in effective temperature and radius, maintaining a nearly constant luminosity. A similar cycle time was found by de Groot et al. (2001), who reported on photometric variations which were consistent with a timescale of 5.5 to 8.5 years.

In addition to the large scale variations in emission strength, Markova (2000) found that there are at least four other kinds of line profile variability in the spectrum of P Cygni. The most striking of these is the long documented appearance of blueward-migrating, absorption sub-features that are called Discrete Absorption Components (DACs: Israelian & de Groot 1999; Markova 2000). These are generally observed in low and intermediate excitation state lines in the optical (Markova 2000) and UV spectrum (Israelian et al. 1996). They are frequently detected in the upper sequence of the H Balmer lines (principal quantum number ; Markova 2000), but to our knowledge, DACs have not been reported before now for the absorption component of H. DACs are often (but not always) narrow (FWHM km s) and may be unresolved in low dispersion spectra. The DACs tend to appear over a radial velocity range of to km s with an acceleration of to km s d. A recurrence timescale of  d is sometimes observed (Kolka 1983; Markova 1986a; Israelian et al. 1996; Kolka 1998). These accelerations are much slower and the timescales are much longer than those associated with DACs in the winds of O-stars (Kaper et al. 1999). The DACs in the spectrum of P Cygni may form in outward moving and dense shells (Kolka 1983; Lamers et al. 1985; Markova 1986a; Israelian et al. 1996), in spiral-shaped co-rotating interaction regions (CIRs; Cranmer & Owocki 1996; Markova 2000), or in dense clumps in the wind (Lépine & Moffat 2008).

In this paper, we present new high resolution H spectroscopy, which we combined with previous measurements by Markova et al. (2001a) to explore the characteristics of P Cygni’s short SD-phase. We also compare this with archival Johnson photometry and new observations obtained by AAVSO observers. Section 2 describes our observations. In Section 3, we present our analysis of long-term variations of the continuum and the H equivalent width. We describe the H profile morphology changes and DAC propagations in Section 4. Our discussion and conclusions are presented in Section 5.

## 2 Observations

We obtained 126 new spectroscopic observations of P Cygni using the Ritter Observatory 1 m telescope and échelle spectrograph (Morrison et al., 1997) between 1999 June 7 and 2007 October 30. These high resolving power () spectra were reduced by standard techniques with IRAF666IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.. Observations collected prior to 2007 were taken using the setup described in Morrison et al. (1997). These observations record a 70 Å range in the order centered on H, and they typically have a signal-to-noise ratio between 50 and 100 per resolution element in the continuum. Observations made during the calendar year 2007 were taken with the same spectrograph, except the camera was a Spectral Instruments 600 Series camera, with a front-illuminated Imager Labs IL-C2004 41004096 pixel sensor (1515 micron pixels). To maintain consistency with the older observations, the camera was operated with 22 pixel binning. The newer observations recorded a larger portion of the order centered on H and typically reached a signal to noise ratio between 50 and 100 per resolution element in the continuum. We trimmed these spectra so that the wavelength range was the same as in the older data. The spectra taken after 2002 September have poor wavelength calibration due to problems with the Th-Ar lamp. In order to use these spectra for kinematical measurements, the telluric HO lines in the vicinity of H were fitted to improve the solution. This worked for most cases, but the errors associated with the telluric re-calibration are roughly km s, compared with the errors for earlier data of km s.

We collected -band photometry from three sources. The first was from Markova et al. (2001a). This provided concurrent photometry for the H data previously published. The second set came from Percy et al. (2001)777Available for download at http://schwab.tsuniv.edu/papers/paspc/pcyg/pcyg.html. These observations also ended at nearly the same time as the first data set. Finally, we downloaded the photoelectric photometry in the -band from the American Association of Variable Star Observers (AAVSO). The AAVSO data are helpful in understanding the long-term trends, and we only used data where measurements of the check and comparison stars differed from expected values by less than 0.05 mag. The errors in the AAVSO measurements are typically around 0.01 mag, comparable to those of Markova et al. (2001) and Percy et al. (2001). The combined set contains 3142 measurements from 1985 to 2009.

## 3 The Long-Term Photometric and Hα Equivalent Width Variability

Markova et al. (2001a) found that the H line flux, obtained by correcting the observed equivalent widths for the changing continuum, varied in concert with the -band flux over the period from 1989 to 1999 (see their Fig. 3). Here we extend their work by considering the long-term photometric and H variations through 2007. Figure 1 shows the large time span of available photometry. The light curve over this interval shows rather modest, mag variations, consistent with the star’s classification by van Genderen (2001) as a “weak-active” LBV. It exhibits the kind of variability associated with a short SD-phase, similar to that reported by Markova et al. (2001a). The short SD-phase is most evident in the data prior to 2000, when the star experienced two fadings of mag (Markova et al., 2001a; Percy et al., 2001). The amplitude of this long-term variability decreased in subsequent years, which indicates that the properties of the short SD-phase change with time. We made a fit of the very-long-term trend and found a brightening rate of mag century (overplotted in Fig. 1). This rate is consistent with the very-long-term trend of mag century documented by de Groot & Lamers (1992).

Figure 2 presents the “dirty” discrete Fourier Transform (Roberts et al. 1987) of the 24.5 y -band photometry with the long-term brightening (Fig. 1) removed. There are no individual significant peaks in the periodogram, but there is a general tendency for more power to appear at the longer timescales (lower frequencies). Thus, the longer timescales of the short SD-phase variability tend to dominate the light curve.

We measured H emission strength for both the new and originally reported spectra (Markova et al. 2001a) for a total of 158 measurements covering the interval from 1994 to 2007. Equivalent widths of the full H profile (including both the blue absorption and large emission component) were measured in the same manner as done by Markova et al. (2001a) in order to keep the data sets mutually consistent. The only improvement is that telluric HO lines in the vicinity of H were removed by means of a template fitting procedure (telluric) in IRAF. This correction resulted in equivalent width increases of less than 2%. This is much smaller than the typical measurement error of 6% (as found by comparing equivalent widths from closely spaced observations, d, where the variability of this star is minimal). Since the available wavelength range around H alpha does not extend beyond the electron scattering line wings to the actual continuum levels, a multiplicative constant was used to retrieve the full equivalent width of the line. This correction, (net) = 1.096 (Ritter), accounts for unseen line wing flux and unmeasured flux lying below our continuum placement (over an integration range of 6531.5 to 6593.5Å) and is identical to that adopted by Markova et al. (2001a) for the Ritter data. The heliocentric Julian dates and net adjusted equivalent widths are tabulated in columns 1 and 2 of Table 1.

The actual line flux of H can be estimated by correcting the measured equivalent widths for the changing continuum. In order to make this transformation, we averaged the -band measurements made within 20 days of each spectroscopic measurement. This time span was chosen to cancel any small but fast variations and to include enough measurements for a reliable average. We compared all the photometry measurements to a benchmark =4.8 to remain consistent with the flux correction adopted by Markova et al. (2001a). The equivalent widths were corrected using the relationship

 Wλ(corr)=Wλ(net)10−0.4(V(t)−4.8)).

The averaged magnitudes and flux corrected equivalent widths are given in columns 3 and 4 of Table 1. If no magnitude was available within days, then no correction was applied, which affects 18 of our measurements. These correction factors are usually small () and comparable to the photometric scatter within each time window.

We show the temporal variations in the flux corrected equivalent widths in Figure 3. The plot includes earlier measurements from Markova et al. (2001a), the new measurements from Table 1, and some additional measurements from 2005 to 2007 from Balan et al. (2010). There are two maxima (occurring around 1992 and 2002) that are separated by years, which is longer than the reported lengths of the short SD-phase found by Markova et al. (2001a), de Groot et al. (2001), or Percy et al. (2001). Furthermore, the rise and fall around the peak in 1992 are steeper than that for the 2002 peak. There is also ample evidence of faster variability within each observing season that appears to be unrelated to the longer term trends.

A visual comparison of the -band photometry in Figure 1 with the flux-corrected H equivalent widths in Figure 3 immediately shows some variations in common. We found that the relative flux (from the time interpolated magnitude) is positively correlated with the corrected equivalent width. A linear fit yields a slope of , confirming the visual impression of co-variability. In order to compare directly the photometry and H equivalent widths, we removed the long-term linear trend from the photometry (Fig. 1), performed a running average of the photometry differences using a Gaussian weighting scheme parametrized by a Gaussian FWHM, transformed the resulting flux differences into a variation in Å units according to the correlation slope given above, and then added the mean equivalent width to the final result. We made a number of trial comparisons by varying the adopted Gaussian FWHM to smooth the photometry, and the best fit with FWHM = 598 d is shown as a solid line in Figure 3. For completeness, we also plot a similar running average of the H equivalent widths as a dashed line. The agreement between the temporally smoothed photometric and flux corrected H variations is striking and it appears to confirm the positive correlation first noted by Markova et al. (2001a). The fact that smoothing parameter values smaller than 598 d yield worse fits suggests that the co-variations are less correlated on shorter timescales. Taken at face value, this result indicates that the continuum and H emission fluxes sample structures in the wind in different ways (probably because of different sites of formation in the outflow). Finally, we note that the correlation also exists between the running averages of the continuum flux and the uncorrected equivalent widths (net), so the covariations are unrelated to the flux correction procedure.

## 4 Hα Profile Morphology Variability

### 4.1 Emission Component Changes

The large H equivalent width variations described in §3 are accompanied by changes in the morphology of the profile. We present two individual spectra in Figure 4 that represent the extrema of the equivalent widths observed (a minimum at HJD 2,450,004, plotted with a dashed-dot line, and a maximum at HJD 2,452,070 plotted with a solid line). It is clear that the profile experiences a change in the peak emission intensity, the line width, the net profile velocity, and the shape of the blue absorption trough. For each of the spectra collected at Ritter Observatory we measured the peak intensity above continuum level, , which we corrected for the changing continuum level in the same manner as the equivalent width (§3), the FWHM (profile width) of the emission portion of the profile above continuum, and a relative radial velocity , derived by cross-correlating each profile against an unweighted average of all the spectra obtained at Ritter Observatory. We chose to use a cross-correlation technique because this method is model-free and is most sensitive to the steep emission line wings, resulting in a measure similar to a FWHM bisector velocity. The resulting measurements of FWHM and are shown for these two profiles in Figure 4 with horizontal and vertical lines, respectively.

All these measurements are given in columns 5, 6, and 7 of Table 1, and are plotted as a function of time in the three panels of Figure 5. We see that times of strong emission (for example, HJD 2,452,500; see Fig. 3, MJD 5.25) correspond to profiles with the largest peak intensity and smallest FWHM. We also see a small radial velocity shift that is correlated with the long-term variations. The profile had the largest (most positive) velocity when the line flux at the position of the emission peak was strongest, which was also when the profile showed the smallest FWHM (Fig. 5). This is likely due to changes in the P Cygni absorption component. When the profile has the most emission, the blue absorption portion appears to shift to a more positive velocity and removes more of the blue side of the emission peak (see Fig. 4), and thus, the net radial velocity tends toward a larger (more positive) value at those times. We find that the measurement errors for and the FWHM are about km s, which adds in quadrature to the calibration errors discussed in §2, to yield net errors of approximately km s for most of the data, and km s for data taken after 2002 September. The errors for are on the order of .

Kashi (2010) has suggested P Cygni is a binary system with a fainter B-type companion and that small long-term radial velocity variations due to reflex motion might be observed in extended high resolution spectroscopic observations. This cannot be the explanation for the changes we observe, since the H emission is formed over a volume that is much larger in radius than the predicted semimajor axis of the putative orbit of the P Cygni primary star. Wind gas leaving the star at any instant would have a Keplerian orbital component that decreases with distance from the center of mass. As the gas packet moves out to the radius where H becomes optically thin and emits the photons we observe (at and larger; see below), the radial outflow component will increase by radiative driving while the orbital motion component will drop with distance to conserve angular momentum. Thus, at the large distance of line formation, the gas motion will be almost completely radial. If the putative companion is to be found from radial velocity variations of this star, then detailed analyses of photospheric or wind lines formed very close to the star will need to be analyzed. Further, these radial velocity variations are not strictly periodic, and cannot be considered orbital motion. Lastly, given the method of measuring these velocities, the measured radial velocity is at least partially due to morphological changes in the H line profile.

### 4.2 Blue Absorption Changes

We used all of the H spectra from Ritter Observatory, including those that were measured by Markova et al. (2001a), to investigate the variations in the blue absorption trough of the P Cygni profile. This portion of the profile is especially interesting as it is formed in the outflowing gas along our line of sight to the star. In order to emphasize the relative changes in line absorption, we first formed a reference, average high-intensity, minimum-absorption spectrum, as follows. At each wavelength step, we ordered the time-series by intensity and then constructed the average of all the intensities falling between the 90th and 95th percentiles at that wavelength step. This removed any spurious peaks caused by cosmic rays from contaminating the minimum absorption average. We then divided each of the spectra by this reference spectrum to form a matrix of quotient spectra. Since we are interested in the variability of the central absorption, and not that of the far wings, and because the line wings never reach the continuum in the region we recorded, the quotient spectra had a depressed continuum. We then re-normalized these to a unit continuum (outside of the velocity region km s). These spectra are illustrated in a gray-scale dynamical spectrum in Figure 6. In this figure, we present each quotient spectrum as a function of radial velocity and time with a gray-scale intensity between the minimum value (black; in the quotient) and maximum value (white; 1.75 in the quotient) based upon a linear time interpolation between the nearest observations (indicated by arrows). The low absorption reference spectrum is displayed for comparison in a panel below the gray-scale image. For simplicity, these quotient spectra were not corrected for the variable continuum flux since we are interested in both emission and absorption changes.

We need to bear in mind that the low absorption spectrum was formed by different subsamples at each wavelength point, and this has important consequences for the appearance of the dynamical spectrum. For example, inspection of Figure 4 shows that blue absorption core can extend to high negative velocities (dash-dot line) while at other times the blue absorption is limited to moderate velocities (solid line). Thus, the construction of the low absorption spectrum will be dominated by the latter examples in the vicinity of the blue absorption edge, and in our collection of H spectra, the more extended blue absorption occurred much more frequently. Consequently, the quotient spectra in Figure 6 appear to be dominated by a blue absorption feature, near km s, except near HJD 2,452,500 (MJD 5.25) when the blue edge moved to a more positive velocity. This feature is due entirely to our selection of spectra from around HJD 2,452,500 in making the low absorption spectrum.

We also see evidence in Figure 6 of several blueward moving, absorption features (primarily between and km s) that appear morphologically similar to the Discrete Absorption Components (DACs) observed previously in other spectral lines (Israelian & de Groot 1999; Markova 1986a, 2000). Figure 7 is a montage of a subset of the quotient spectra. It shows how the DAC (center) moved progressively blueward over this time span ( d). There are times where the regions between successive DACs appear bright in the dynamical spectrum, which correspond to those cases (usually sparsely sampled in time) where the flux was higher than the mean in the 90% to 95% part of the flux distribution that defined the minimum absorption spectrum. Finally, we see in the gray-scale image of Figure 6 the long-term variations in peak intensity () and wing extension (FWHM) that are associated with the short SD-phase, equivalent width variations (Fig. 5).

We measured the variability of the quotient spectra by calculating the standard deviation at each pixel of velocity space. This standard deviation spectrum is shown in Figure 8. There is a broad feature centered at rest velocity which is associated with the varying peak height of the profile. Another feature is present at km s, which could be caused by either the variations in the profile width (FWHM) or red emission wing variability from traveling bumps (Markova 2000). The largest two features are from the DACs (seen as a broad peak around km s) and from the variations present near the blue edge of the absorption core (visible as a peak centered at km s).

Figures 4, 6, 7, and 8 all demonstrate that there are significant changes observed in the profile near the blue edge of the absorption core. The blue-edge velocity of a P Cygni profile is generally set by the position where the absorption core rises to intersect with the local continuum, and this velocity corresponds to the fastest moving gas projected against the disk of the star. However, in the case of H in the spectrum of P Cygni, this location is often poorly constrained because absorption may extend blueward with a shallow slope (see Fig. 4). Thus, we decided instead to document the kinematical changes near the blue edge by measuring the position of the absorption core flux minimum, , which is normally found near km s where the slope of the profile changes sign abruptly. We determined this position by finding the zero crossing in the numerical derivative of a smoothed version of the spectrum. The S/N ratio was sufficient in all our spectra that the zero of the derivation was always well-defined. This estimate of the minimum flux velocity (min) is given in column 8 of Table 1, and the errors in are comparable to the errors associated with the emission line kinematic measurements (§4.1). This velocity is probably related to the wind speed at a location in the wind where H ceases to be optically thick.

It is difficult to measure the radial velocities of the DACs because their morphologies vary and because the absorption may consist of multiple components. We decided to measure a centroid for the DACs wherever possible by means of the relationship

 Vr(DAC)=∫v2v1vr(1−Q(vr)) dvr∫v2v1(1−Q(vr)) dvr

where represents the quotient spectrum and is the radial velocity. We adopted a velocity range of km s and km s based upon the strongest regions of the standard deviation spectrum plotted in Figure 8. Typical errors in these measurements (from the scatter in densely sampled regions of the time series) are km s. This approach worked for most of the spectra, but some low contrast features were not measured correctly, and are omitted from Table 1 and our analysis. During some epochs, there were multiple DACs present, so represents a weighted average of multiple components. The DAC radial velocities are given in column 9 of Table 1. Column 10 lists a relative equivalent width for the DAC measured by a direct numerical integration of between and . The typical errors associated with these equivalent width measurements are on the order of 5%, which is similar to the errors associated with the equivalent widths of the profile (§3).

We show the temporal evolution of the DAC and blue-minimum flux velocities in Figure 9. At some epochs, a long progression of DAC velocities is present. The most well defined sequence began near HJD 2,450,700, and lasted about 800 days (Fig. 7). This timescale is much longer than that of a typical DAC observed in an O star, where the progression is measured in hours or days. As this particular DAC was relatively narrow and recorded in many spectra over its duration, it is an optimal feature to measure the H DAC acceleration. From a simple linear fit, shown overplotted in Figures 7 and 9, we measured the acceleration to be km s d. For comparison, the observed DACs in the spectrum of a normal hot supergiant of similar spectral type, Ori (HD 37128; B0 Ia), have an acceleration of km s d (Prinja et al. 2002).

We used the Lomb-Scargle periodogram method (Scargle 1982) to search for a recurrence time in the appearance of the DACs, and we derived a cycle time of 1700 d. This recurrence timescale is shown in Figure 9 where we overplot the linear acceleration of the DAC shown in Figure 7 for three additional epochs over the time span of the Ritter data. There is some evidence that a DAC progression is seen at each of these four epochs. The major deviations from the expected velocity trends occur when multiple components are present. For example, there were two DACs present between HJD 2,451,700 and 2,452,400, and the velocity centroid we measured represents a blend of these components. Neither of these two DACs occurred at the expected recurrence time in the 1700 d cycle.

Our work represents the first detection of DACs in the blue absorption trough of H, and their properties differ from those observed in other spectral lines (Israelian et al. 1996; Markova 1986a, 2000). For example, the recurrence timescale of 1700 d is much larger than the 200 d interval found in earlier work, and the acceleration we measure is about a factor of 10 smaller than measured by others. We suspect that these differences are due to the large optical depth of the H line compared to that of other lines where DACs have been investigated. This will mean that the radius of optical depth unity is larger for H (Najarro et al. 1997), and consequently any wind structures formed at smaller radii will have no affect on the H line formation. We suspect that observations of other, less optically thick lines are more sensitive to the detection of DACs formed at smaller radii in the wind of P Cygni where more and faster accelerating structures may exist.

## 5 Discussion

The -band and H variations we observe need to be interpreted in the context of current models for the star and its wind. Langer et al. (1994) discuss the atmospheric properties of P Cyg, and they argue that the star is in the LBV phase where the temperature and helium abundance are increasing, and the mass and luminosity are decreasing, as the star evolves towards the Wolf-Rayet phase. Langer et al. emphasize the earlier conclusion from Pauldrach & Puls (1990) that the atmospheric parameters are close to a bi-stability point where the mass loss rate can change by an order of magnitude with small changes in radius and/or luminosity, which may explain the great eruptions observed in prior centuries. The atmospheric parameters are well established through a detailed quantitative spectroscopic analysis by Najarro et al. (1997) and Najarro (2001), who find that He is overabundant and that the mass loss rate is high (  y including wind clumping effects) and wind terminal velocity is low ( km s). Najarro et al. (1997) derive a systemic velocity of km s, and thus our minimum measurement of km s is consistent with their estimate of km s as this velocity measurement is related to . They estimate that the continuum forming radius is , which for a distance of 1.8 kpc implies an angular diameter of mas. On the other hand, Najarro et al. (1997) predict that the emitting size of H will be much larger because of its greater optical depth. For example, their models show that there is a local maximum in the wind temperature distribution (presumably where the recombination processes that form H peak) near (see their Fig. 5b). The corresponding angular size for H of  mas agrees well with the range of  mas from H interferometry by Balan et al. (2010). Thus, we need to keep in mind that the H variations reflect changes over a much larger spatial scale in the wind than those observed in the -band flux.

Variations in the H emission equivalent width are related to changes in both the mass loss rate and the wind velocity. In a very simplified approach, we can assume that most of the H flux originates in the optically thick region projected on the sky,

 f=2πr2τF(T)

were is the boundary separating the optically thick and thin regimes, is the monochromatic surface flux, and is the wind temperature at (Najarro et al. 1997). If we assume that the wind is approximately isothermal at this physical location (a reasonable choice: see Fig. 5b in Najarro et al. 1997), then the emission flux variations are due to changes in the projected size of the optically thick region,

 △f/f=2△rτ/rτ.

Thus, we expect that the relative variations in angular size will be only half as large as the emission equivalent width variations, which is probably consistent with the lack of measurable size changes in the H interferometric measurements (Balan et al. 2010).

The H optical depth is dependent on the electron density squared since the emission is a recombination process. Thus, we expect that the optical depth unity boundary will always be defined by the location in the wind with a specific characteristic density, . We assume that has an approximately constant value so that the effective boundary will vary as fluctuations in the wind mass loss rate and velocity define the radius where the density reaches . According to the mass continuity equation, is related to this density by

 r2τ=˙M4πρτv

where is the mass loss rate and is the wind velocity at the radial distance . We can differentiate the mass continuity equation to express the radius variation in terms of the changes in and ,

 2rτ△rτ=14πρτ△[˙M/v],

which we divide by to obtain

 2△rτrτ=△[˙M/v][˙M/v].

Since we argued above that the flux also varies as , we can then use the relation above to re-write the fractional flux variation in terms of logarithmic changes in and ,

 △lnf=△ln˙M−△lnv.

Since we have observational data on the variations in emission strength and wind velocity, we can rearrange this equation to solve for the mass loss variations as a function of flux variations,

 △ln˙M/△lnf=1+△lnv/△lnf.

Puls et al. (1996) present a much more detailed analysis of the dependence of the emission equivalent width on the wind parameters of hot, massive stars. However, in the limit of high optical depth, their expression for the emission flux (their eq. 41) leads to a similar relation,

 △ln˙M/△lnf=34(1+△lnv/△lnf).

We found in the previous section that the H equivalent width appears to vary inversely with two quantities related to wind dynamics, the H emission peak FWHM and the blue minimum flux velocity (min). Figure 10 quantifies this relationship. The upper panel shows the inverse correlation between the emission peak FWHM and flux corrected H equivalent width, and if we take FWHM as a proxy for the wind speed, then a linear fit of natural logarithms of these measures gives . We caution that the FWHM is also influenced by the absorption component of H, and we showed above (Fig. 4) that the absorption component moves inward towards the line core when the emission is strong. Consequently, the apparent decrease in FWHM as the emission increases probably results both from a wind speed decrease and a blue wing decline due to blending with the absorption component. The minimum flux velocity is perhaps a more direct measurement of wind speed (at least in our line of sight), and we show in the lower panel of Figure 10 the co-variations of the difference between (min) and the systemic velocity of P Cyg, km s (Najarro et al. 1997), as a function of the corrected equivalent width. A fit of the logarithmic slope here gives a smaller estimate of .

If we adopt the minimum flux co-variation result as representative of the wind velocity component of variability, then the mass loss rate variation we derive from the relation above has an emission flux dependence of . Omitting the bottom and top 10% of the distribution of (corr), the derived range in emission strength in our observations of probably implies mass loss rate changes of . Markova et al. (2001a) used the optically thick relation from Puls et al. (1996) to arrive at an estimate of for the mass-loss variation amplitude. We showed above that the relation from Puls et al. carries a factor of that is missing from our simple analysis, and if we use the Puls et al. relation instead, we arrive at a mass loss rate variation of , confirming the estimate from Markova et al. (2001a). In this scenario, the H emission variations result from changes in the effective emission radius caused by variations in the mass loss rate and wind velocity. During episodes when the mass loss rate is higher and the wind is slower, the projected size of the emission region increases leading to larger H emission flux. The same process probably causes the -band variations, but the fractional radius variations must be smaller at the continuum forming radius because we found in last section that so that , i.e., the continuum size variations are only as large as those in H. It is possible that the changes in the mass-loss rate are caused by a change in the luminosity of the star which would propagate through the wind, and be observed in both the -band brightness and the emission line flux of H.

Finally we return to the relationship of the DACs to the short SD-phase variations. We found a trend with the DAC strength and the H emission equivalent width. We show in Figure 11 the temporal behavior of DAC quotient equivalent width (DAC) (Table 1, column 10) along with scaled, running averages of the H equivalent width and -band flux (Fig. 3). We rescaled the amplitude of the H flux by and the photometric light curve was rescaled by . These curves were then shifted vertically to match the points. We see that the DAC strength variations track both the H and -band flux variations, suggesting that the DACs are related in some way to the short SD-phase changes. For example, we see that some of the strongest DACs were observed when the H emission was strong (around HJD 2,452,500; MJD 5.25 in our plots) and the DACs were very faint or absent when the emission was weak (around HJD 2,450,500).

The DAC phenomenon is primarily observed in the UV wind lines of hot stars (Kaper et al. 1999; Puls et al. 2008), and, in fact, DACs are only rarely seen in the H profile where wind-related variations are usually due to changes in the dense and slower moving wind close to the star (Kaper et al. 1997; Markova et al. 2005). The only comparable DAC observed in H was in HD 92207 (Kaufer et al. 1996) with a lifetime of d and only one DAC observed, so the recurrence time is unknown. This star is cooler than P Cygni and is a “normal” supergiant, compared to the luminous blue variable nature of P Cyg. The DACs observed in the UV wind lines of O-stars are first seen at velocities of 0.2 to 0.4 and then they migrate blueward to on time scales of a day or so, exhibiting accelerations that are much slower than expected for the wind velocity law (Kaper et al. 1999). Many of these same features are seen in the DACs in H for P Cygni, although they occur on vastly longer time scales. For example, the wind flow time scales as , and while the wind gas will accelerate from to in 0.5 d for an O-star like  Per (Kaper et al. 1999), it will take some 43 d for the wind of P Cygni. However, this flow time is very short compared to the longevity of the DACs ( d), so we are led to the same conclusion found for the O-stars, namely that the DACs represent some kind of perturbation in the wind through which the gas flows. Changes in wind velocity and/or mass loss rate can cause shocks and create structure in the wind, and these structures produce density enhancements and/or velocity plateaus that imprint DACs in the wind lines (Fullerton & Owocki 1992; Runacres & Owocki 2002; Puls et al. 2008).

Current theory suggests that the DACs are related to large scale spiral features in the wind known as co-rotating interaction regions that are formed at the intersection of fast and slower outflows, which develop from some inhomogeneity in the stellar photosphere (for example, pulsation or spots; Cranmer & Owocki 1996; Lobel & Blomme 2008). In these models, it is the slow transit of these equatorially centered regions across the photosphere that creates the DACs in the absorption cores but has little influence on the emission parts of the wind line (Dessart 2004; Lobel & Blomme 2008). However, in the case of P Cygni, we find that emission parts do appear to strengthen when DACs are prominent (Fig. 11), and this indicates that the wind perturbation profoundly affects both wind gas surrounding the star and the wind gas projected against the photosphere. Thus, we suggest that the structures causing the DACs in P Cygni may be more spherically symmetric than assumed in the geometry of the co-rotating interaction regions. In some models the seed perturbation occurs at a fixed longitude on the star, so that a new wind structure appears each time the star rotates (although for the best studied case of HD 64760, Lobel & Blomme 2008 argue that the originating spots must rotate some five times slower than the star in order to fit models to the observations). If this is the case for P Cygni, then the 1700 d DAC recurrence time may be related to the star’s rotational period. Markova (2000) found an upper limit of 100 d for the rotational period based upon estimates of the stellar radius and the projected rotational velocity. However, the line broadening in early supergiants may be dominated by turbulence rather than rotation (Howarth et al. 1997, 2007; Markova & Puls 2008), so a longer rotational period remains a possibility. However, regardless of the origin of the DACs, their close relation with the emission line and continuum flux variations (Fig. 11) suggests that much of the short-SD variability is caused by propagating structural perturbations in the outer atmosphere of the star.

The discovery of the short SD-phase variations in P Cygni and its relationship to the wind velocity and DAC occurrence in the wind is a new observational result that warrants future investigation both for this star and other LBVs. The changing characteristics over these long timescales may eventually lead to a better understanding of the LBV stage of evolution and the underlying physics of their winds and circumstellar environments. We are currently pursuing a three year, spectroscopic monitoring program of Galactic and Magellanic Cloud LBVs. Such long-term observations will reveal if DACs are present in all LBVs and will show whether the variations found in P Cygni are a general phenomena among LBVs.

We are grateful to an anonymous referee for important and helpful comments on an earlier version of this paper that greatly improved the quality of the paper. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. Support for Ritter Astrophysical Research Center during the time of the observations was provided by the Fund for Astrophysical Research, the American Astronomical Society Small Grants Program as well as the National Science Foundation Program for Research and Education with Small Telescopes (NSF-PREST) under grant AST-0440784 (NM). This work was also supported by the National Science Foundation under grants AST-0606861 and AST-1009080 (DG). Institutional support has been provided from the GSU College of Arts and Sciences and from the Research Program Enhancement fund of the Board of Regents of the University System of Georgia, administered through the GSU Office of the Vice President for Research. N. Markova acknowledges the Bulgarian NFSR grant DO 02-85. We are grateful for all this support. Facilities: Ritter Observatory AAVSO

## References

• Balan et al. (2010) Balan, A., et al. 2010, AJ, 139, 2269
• Cranmer & Owocki (1996) Cranmer, S., & Owocki, S. 1996, ApJ, 462, 469
• de Groot (1990) de Groot, M. 1990, in Properties of Hot Luminous Stars (ASP Conf. Ser. 7), ed. C. D. Garmany (San Francisco: ASP), 165
• de Groot & Lamers (1992) de Groot, M., & Lamers, H. J. G. L. M. 1992, Nature, 355, 422
• de Groot et al. (2001) de Groot, M., Sterken, C., & van Genderen, A. M. 2001, A&A, 376, 224
• Dessart (2004) Dessart, L. 2004, A&A, 423, 693
• Fullerton & Owocki (1992) Fullerton, A. W., & Owocki, S. P. 1992, in Nonisotropic and Variable Outflows from Stars (ASP Conf. Ser. 22), ed. L. Drissen, C. Leitherer, & A. Nota (San Francisco: ASP), 177
• Howarth et al. (1997) Howarth, I. D., Siebert, K. W., Hussain, G. A. J., & Prinja, R. K. 1997, MNRAS, 284, 265
• Howarth et al. (2007) Howarth, I. D., et al. 2007, MNRAS, 381, 433
• Humphreys & Davidson (1994) Humphreys, R. M., & Davidson, K. 1994, PASP, 106, 1025
• Israelian & de Groot (1999) Israelian, G., & de Groot, M. 1999, Space Sci. Rev., 90, 493
• Israelian et al. (1996) Israelian, G., de Groot, M., Parker, J. W., & Sterken, C. 1996, MNRAS, 283, 119
• Kaper et al. (1999) Kaper, L., Henrichs, H. F., Nichols, J. S., & Telting, J. H. 1999, A&A, 344, 231
• Kaper et al. (1997) Kaper, L., et al. 1997, A&A, 327, 281
• Kashi (2010) Kashi, A. 2010, MNRAS, 405, 1924
• Kaufer et al. (1996) Kaufer, A., et al. 1996, A&A, 305, 887
• Kolka (1983) Kolka, I. 1983, Publ. Acad. Sci. Estonian SSR, 32, 51
• Kolka (1998) Kolka, I. 1998, in Cyclical Variability in Stellar Winds, ed. L. Kaper & A. W. Fullerton (Berlin: Springer-Verlag), 111
• Lamers & de Groot (1992) Lamers, H. J. G. L. M., & de Groot, M. 1992, A&A, 257, 153
• Lamers et al. (1985) Lamers, H. J. G. L. M., Korevaar, P., & Cassatella, A. 1985, A&A, 149, 29
• Langer et al. (1994) Langer, N., Hamann, W.-R., Lennon, M., Najarro, F., Pauldrach, A. W. A., & Puls, J. 1994, A&A, 290, 819
• Lépine & Moffat (2008) Lépine, S., & Moffat, A. F. J. 2008, AJ, 136, 548
• Lobel & Blomme (2008) Lobel, A., & Blomme, R. 2008, ApJ, 678, 408
• Luud et al. (1967) Luud, L. S. 1967, Sov. Astron., 11, 211
• Markova (1986a) Markova, N. 1986a, A&A, 162, L3
• Markova (1986b) Markova, N. 1986b, Ap&SS, 123, 5
• Markova (1993) Markova, N. 1993, Ap. Space Sci., 201, 61
• Markova (2000) Markova, N. 2000, A&AS, 144, 391
• Markova et al. (2001a) Markova, N., Morrison, N., Kolka, I., & Markov, H. 2001a, A&A, 376, 898
• Markova & Puls (2008) Markova, N., & Puls, J. 2008, A&A, 478, 823
• Markova et al. (2005) Markova, N., Puls, J., Scuderi, S., & Markov, H. 2005, A&A, 440, 1133
• Markova et al. (2001b) Markova, N., Scuderi, S., de Groot, M., Markov, H., & Panagia, N. 2001b, A&A, 366, 935
• Morrison et al. (1997) Morrison, N. D., Knauth, D. C., Mulliss, C. L., & Lee, W. 1997, PASP, 109, 676
• Najarro et al. (1997) Najarro, F., Hillier, D. J., & Stahl, O. 1997, A&A, 326, 1117
• Najarro (2001) Najarro, F. 2001, in P Cygni 2000: 400 Years of Progress (ASP Conf. Ser. 233), ed. M. de Groot & C. Sterken (San Francisco: ASP), 133
• Pauldrach & Puls (1990) Pauldrach, A. W. A., & Puls, J. 1990, A&A, 237, 409
• Percy et al. (2001) Percy, J. R., Evans, T. D. K., Henry, G. W., & Mattei, J. A. 2001, in P Cygni 2000: 400 Years of Progress (ASP Conf. Ser. 233), ed. M. de Groot & C. Sterken (San Francisco: ASP), 31
• Percy & Welch (1983) Percy, J. R., & Welch, D. L. 1983, PASP, 95, 491
• Percy et al. (1988) Percy, J. R., et al. 1988, A&A, 191, 248
• Prinja et al. (2002) Prinja, R. K., Massa, D., & Fullerton, A. W. 2002, A&A, 388, 587
• Puls et al. (2008) Puls, J., Vink, J. S., & Najarro, F. 2008, A&A Rev., 16, 209
• Puls et al. (1996) Puls, J., et al. 1996, A&A, 305, 171
• Roberts et al. (1987) Roberts, D. H., Lehar, J., & Dreher, J. W. 1987, AJ, 93, 968
• Runacres & Owocki (2002) Runacres, M. C., & Owocki, S. P. 2002, A&A, 381, 1015
• Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835
• Sterken et al. (1997) Sterken, C., van Genderen, A. M., & de Groot, M. 1997, in Luminous Blue Variables: Massive Stars in Transition (ASP Conf. Ser. 120), ed. A. Nota & H. Lamers (San Francisco: ASP), 35
• van Genderen (2001) van Genderen, A. M. 2001, A&A, 366, 508
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters