© ESO, 2011

1. Introduction

Luminous blue variables (LBVs) are unstable massive stars whose evolutionary state is under heavy debate. While they have been thought to be in a short, transitory phase between an O-type and Wolf-Rayet (WR) star (Humphreys & Davidson 1994), recent observations suggest that LBVs might actually be the final stage before the star explodes as a type-II SN (e. g. Kotak & Vink 2006; Smith et al. 2007; Gal-Yam & Leonard 2009).

LBVs are characterized by severe mass loss (~10^-6−10^-3   M_⊙   yr^-1), slow winds (terminal speeds  ~50−500   km   s^-1), photometric and spectroscopic S Doradus-type variations on timescales of years, and are often surrounded by circumstellar nebulae (Humphreys & Davidson 1994). During the S Dor-type variations, the stellar radius changes by up to a factor of 10, while the effective temperature changes in the range 8000 − 25   000   K (van Genderen 2001). There is evidence for variability in the bolometric luminosity of up to a factor of two in the few LBVs that have been analyzed in detail with modern radiative transfer codes (Groh et al. 2009; Clark et al. 2009).

During the S-Dor variations, LBVs may eventually cross the so-called bi-stability limit (Vink & de Koter 2002; Groh et al. 2011), which occurs in line-driven winds at around T_eff ~ 21   000 K (Lamers et al. 1995; Vink & de Koter 2002). Model predictions by Pauldrach & Puls (1990) suggested that wind bi-stability plays an important role in a shell-ejection scenario of P Cygni, in the sense that small variations (~5%) in L _⋆  or R _⋆  cause significant changes in the wind properties (see also Najarro et al. 1997; Hillier et al. 1998).

Today, bi-stability refers to mass loss and wind-velocity variations resulting from changes in the Fe ionization (Vink et al. 1999). Given that the defining characteristic of LBVs are their S-Dor variations (Humphreys & Davidson 1994), there is an expectation that LBVs could change their mass-loss properties due to the bi-stability whilst they travel back and forth across the upper Hertzsprung-Russell diagram (Vink & de Koter 2002; Lamers 1997). Indeed, Groh et al. (2011) confirmed that the prototype LBV AG Carinae crosses the bi-stability limit, but were unable to identify whether a gradual change or a jump in the wind properties occurs.

It is also the property of having variable mass loss and wind terminal speed that first revealed LBVs as candidate supernova (SN) progenitors (Kotak & Vink 2006). This assertion was strengthened by the detection of double-troughed P Cygni absorption components in the narrow Hα portion of the interacting supernova SN2005gj (Trundle et al. 2008). Such multiple-absorption components have been known, for several decades, to exist in Hα spectra of confirmed LBVs (Fig. 1; e.g., Stahl et al. 1983, 2001, 2003; Leitherer et al. 1994). The phenomenon might be the result of shells ejected by the LBV, or caused by radiative transfer effects (Hillier 1992), or due to variable winds, as would be expected from the physics of the bi-stability jump.

The main aim of this Letter is to explain the origin of the multiple P-Cygni absorption components seen in LBVs. Explaining this phenomenon is particularly relevant as it is as yet far from proven that LBVs really are SN progenitors. A number of hints in this direction have been suggested, but they all have their specific issues (see e.g. Dwarkadas 2011, and references therein). Most significantly, stellar evolution theory does not predict the cores of LBV stars to be advanced enough to be anywhere close to core-collapse (Langer et al. 1994; Meynet & Maeder 2003). It is thus important to scrutinize all suggestions involving a direct link between LBVs and SNe. Before being able to attack the line profiles seen in SN2005gj, which involve the complexities inherent to both progenitor mass loss and the subsequent SNe, it seems logical to first attempt to explain the more widespread existence of multiple absorption components in LBVs themselves.

2. Radiative transfer modeling method

We employ a modified version of the spherical symmetric, fully line-blanketed, non-LTE radiative transfer code cmfgen (Hillier & Miller 1998) to analyze the effects of time-dependent wind parameters on the stellar spectrum. Each model is specified by the hydrostatic core radius R_⋆, luminosity L_⋆, mass-loss rate Ṁ, wind volume filling factor f, wind terminal velocity v_∞, and chemical abundances Z_i of the included species. We compute time-dependent models with the following stellar and wind parameters: T_eff in the range of 13   000 to 25   000   K, Ṁ in the range of 0.2 to 2 × 10^-5   M_⊙   yr^-1, M = 60   M_⊙, L _⋆  = 10⁶   L_⊙, f = 0.25, and chemical abundances typical of LBVs, with a He mass fraction of Y = 0.62, N/N_⊙ = 7, C/C_⊙ = 0.08, O/O_⊙ = 0.025, and solar abundance for the iron group elements. These parameters are roughly characteristic for AG Car.

Here we briefly describe our implementation of time-dependent density and velocity structures. The reader is referred to Groh et al. (2008) and Groh & Hillier (2011, in prep.) for further details. At the moment, no time dependence is accounted for in the solution of the radiative transfer and rate equations, since the recombination timescale is much shorter than the flow timescale for the wind regions studied here. As the wind hydrodynamics is currently not yet solved for, in order to account for time-dependent outflows, we allow for arbitrary v(r) and ρ(r) stratifications as inputs into cmfgen. We consider cases when the star crosses the bi-stability limit, evolving between two epochs: the initial epoch 1, when T_eff = 23   000 K and a fast wind is present with v_∞,1, and epoch 2, when T_eff decreases and a slow wind is present, according to v_∞,2 < v_∞,1. For simplicity, we keep Ṁ fixed, so that ρ(r) is computed following the equation of mass conservation.

We explore two scenarios as to whether the double P-Cygni profiles can be caused by the time variability of wind parameters: a gradual change in v_∞, or an abrupt instantaneous variability.

For models with gradual changes in v_∞(t), for simplicity we assume that, starting at t = 0, v_∞ decreases linearly with time. For a given time Δt after t = 0, we compute the distance traveled by material ejected between t = 0 and Δt. The result is a velocity law that increases linearly with the distance up to a point where the initial wind is reached. From that point onwards, v(r) follows the standard β-type, asymptotically reaching v_∞,1 (Fig. 2a).

For models with an abrupt change, we assume that at t = 0, v_∞ switches instantaneously from v_∞,1 to v_∞,2. To compute a model for a given Δt, we evaluate the distance traveled by the last particles of the initial wind (r₁), as well as the distance traveled by the first particles of the final wind (r₂). Following a suggestion by Owocki (2010, priv. comm.), we assume that the initial wind characterizes the region from r₁ up to the outer boundary of the computational domain, whilst a β-type law characterized by v_∞,1 and β₁ is assumed. From the star up to r₂, we assume that v(r) is given by the velocity law of the final wind according to v_∞,2 and β₂. We further assume that the region in between the two winds, i.e. between r₂ and r₁, v(r) and ρ(r) are linear, and join the two winds smoothly. We defer to future hydrodynamical computations for a more accurate treatment of this region, but we note that initial tests using somewhat different assumptions may lead to small changes in the resulting P-Cygni line profiles. These changes are not anticipated to affect the conclusions reached here.

For both gradual and abrupt models, we assume that the stellar wind is at t = 0 in steady-state with v_∞,1 = 150   km   s^-1 and β₁ = 1 and, as t → ∞, evolves to a final state with v_∞,2 = 70   km   s^-1 and β₂ = 1. The choice of v_∞ values is based on LBV spectral data, whilst the factor of two change in v_∞ is adopted from theoretical predictions (Vink et al. 1999).

3. Time-dependent P-Cygni line profiles: gradual vs. abrupt changes in the wind parameters

Figure 2b presents the P-Cygni line profiles calculated for both a gradual and an abrupt variation of the wind parameters as a function of time, here computed at Δt = 120 days. A gradual change in v_∞(t) causes a linear variation of the velocity with distance (Fig. 2a). At Δt = 120 days, the P-Cygni profile deviates only slightly from a steady-state model with the same stellar parameters. To test the robustness of our results, we computed models for different Δt, as well as stellar and wind parameters all with gradual changes, but we still could not obtain double P-Cygni absorption components.

On the other hand, we found that an abrupt change in v_∞ causes significant effects in the wind velocity law (Fig. 2a). While the velocity law of the outer wind is similar to that of the gradually changing as well as steady-state models, the inner part of the wind follows a β-type law asymptotically approaching v_∞,2. The abrupt change in v_∞ modifies not only the velocity and density structures, but also causes the appearance of a double-troughed P-Cygni absorption component (Fig. 2b). The double P-Cygni absorption component arises as a result of the presence of a small amount of material with velocities in between v_∞,1 and v_∞,2, causing a very low Hα line opacity between these two velocities. Models computed with a variable Ṁ across the bi-stability jump show a similar qualitative behavior as the ones shown here, implying the abrupt change in v_∞ is the key phenomenon responsible for the double-troughed P-Cygni absorption components. We are not able to tell whether there is a jump in Ṁ as well.

Figure 3 presents the temporal evolution of the Hα absorption profile for selected Δt, showcasing how the absorption evolves from the initial fast wind in steady-state (panel a) to the final slow wind in steady-state (panel f, red dashed line). We find that the P-Cygni absorption component from the fast wind becomes weaker, as the slow wind takes over, disappearing at around Δt = 240 days. The P-Cygni absorption component from the slow wind first appears rather weakly at around Δt = 90 days, and increases in strength as a function of time. Note that even when no absorption from the fast wind is seen, its weak emission still influences the spectral morphology, even at Δt = 960 days.

So when can the double P-Cygni absorption components be seen after the abrupt change in v_∞? The duration of the double P-Cygni absorption period depends on the characteristic flow time (dependent on R_⋆/v_∞ and β) of both winds, as well as on the location of the line forming region of interest. For the parameter space explored here, Hα forms at a distance from  ~2 to 10 R _⋆ , and the double P-Cygni absorption component is seen from Δt around 90 to 180 days, with its maximum strength around 120 days. Higher values of β increase the flow timescale and allow the double P-Cygni absorption components to be observed during a more extended period of time.

4. The dependence of the double P-Cygni absorption profiles on T_eff and $\dot{\mathit{M}}$

An abrupt change in v_∞ is required to reproduce double P-Cygni absorption profiles (Fig. 2). The next relevant question is: for which range of stellar and wind parameters does a double P-Cygni absorption profile appear in Hα?

Here we focus on changes in T_eff and Ṁ, keeping the other stellar and wind parameters constant as well as Δt = 120 days. We found that the morphology of the double P-Cygni absorption profiles is strongly sensitive to both T_eff and Ṁ, and that a zoo of absorption line profiles are predicted by our models (Fig. 4).

We note that an increase in T_eff from 13   000   K to 18   400   K causes less absorption by the initial, fast wind, resulting in weaker high-velocity absorption components. Models with T_eff > 18   400   K have H ionized throughout the wind, which dramatically reduces the population of the n = 2 level related to the Hα absorption. As a consequence, these hotter models show extremely weak or no P-Cygni absorption, since the time-dependent velocity structure is not probed by Hα if H is ionized throughout the wind.

In addition, our model predict weak blue-shifted emission between the two absorption components. This blue-shifted emission arises in the initial, fast wind, from non-core rays. According to our models, the intensity of the blue-shifted emission increases as T_eff increases.

Perhaps surprisingly, our models predict that, for a given L _⋆  and T_eff, double-troughed P-Cygni absorption appears in Hα only for a relatively narrow range of Ṁ. We found that there is an optimal value of Ṁ (Ṁ_opt) when the contrast between the two absorption components is maximal. The double-troughed P-Cygni absorptions are present when Ṁ is about  ~± 30% around Ṁ_opt. For T_eff < 18   000   K, we found four possible regimes of the Hα absorption morphology (Fig. 5), depending on Ṁ:

• Ṁ > 1.3   Ṁ_opt: absorption from the fast wind only, since the emission from the fast wind from non-core rays fills out the absorption from the slow wind;

• 0.7   Ṁ_opt < Ṁ < 1.3   Ṁ_opt: double P-Cygni absorption, i.e. absorption from both winds;

• 0.3   Ṁ_opt < Ṁ < 0.7   Ṁ_opt: absorption from the slow wind only since, for lower Ṁ, the fast wind is not optically-thick enough in Hα to produce detectable absorption;

• Ṁ < 0.3   Ṁ_opt: no P-Cygni absorption, since both winds are not optically-thick enough in Hα.

We note that the ranges given above are only indicative and valid for Δt = 120 days, since the relative location of the four different regions described above varies with Δt. In reality, the actual observational detection of the double P-Cygni absorption should also depend on the spectral resolution and the signal-to-noise ratio of the observations, amongst other factors.

5. Discussion

We demonstrated via time-dependent, radiative transfer modeling with cmfgen that abrupt wind-velocity changes are needed to reproduce the observed double P-Cygni absorption components seen in a plethora of LBV spectra. The most promising manner causing such abrupt changes in v_∞ is by means of the bi-stability jump, as discussed in Vink & de Koter (2002).

We reinforce that the modeling performed here is not exhaustive, given the vast parameter space to be explored and the model limitations. For instance, the full hydrodynamics and time-dependent solution of the radiative transfer and statistical equilibrium equations are not yet accounted for. Furthermore, we consider an idealized scenario of an instantaneous abrupt change in v_∞, whereas in reality, the changes could occur over extended time periods of days or weeks. Constraints on the timescale could be obtained via time-dependent modeling of observational data of LBVs that actually cross the bi-stability limit, e.g. AG Car. These topics will be the subject of future studies.

We note that the empirical identification of a luminous massive star as a bona-fide member of the LBV class is a challenging exercise. This is mainly due to the long timescales (years to decades) necessary to monitor the irregular spectral variability. As the detection of double P-Cygni absorption components is a strong indication that the candidate object is an LBV undergoing S-Dor type variability cycles – with its abrupt changes in v_∞ on timescales of years – we suggest that the sheer presence of such double P-Cygni absorption components may become an efficient, novel way to detect Galactic and extra-galactic LBVs.

We found that the morphology of the double P-Cygni absorption components depends strongly on Ṁ and T_eff. This opens up the possibility to constrain these quantities without the lengthy process of non-LTE radiative transfer modeling: once L_⋆ is known, for given Δt, v_∞,1, and v_∞,2 could be derived empirically.

If observations at multiple epochs are available, we should also be able to trace the LBV mass-loss history using a similar modeling strategy as employed here, as the Hα line forms over an extended stellar-wind region in LBVs. This will allow us to use the double P-Cygni absorption components as a tool to look back in time, and to be able to estimate LBV stellar parameters from one single epoch. Although in this Letter we focused on Hα, motivated by the observational data available, we note that our time-dependent models predict similar double-troughed P-Cygni absorption in many other spectral lines formed at different wind regions, such as most of the Balmer lines, Fe ii, N ii, He i λ10830, among others. Analyses of these lines should lead to a much improved understanding of the mass-loss history of LBVs.

Interestingly, the morphology of the absorption-line profiles seen in Fig. 4 closely resemble those found in the interacting SN 2005gj (Trundle et al. 2008). This strengthens the link between LBVs and SNe that based their assertions on the bi-stability mechanism during S Dor cycles (Kotak & Vink 2006; Trundle et al. 2008), and the LBV/SN link in general (Smith et al. 2007; Gal-Yam & Leonard 2009). We expect that our results will encourage further observational and modeling efforts to establish under which conditions LBVs explode as SNe.

Acknowledgments

We thank John Hillier for making CMFGEN available. We acknowledge discussions with Rubina Kotak, Stan Owocki, and Tom Madura. J.H.G. is thankful for the warm hospitality during his visit to Armagh Observatory.

References

All Figures