# Giant Outbursts of Luminous Blue Variables and the Formation of the Homunculus Nebula Around 𝜂 Carinae

## 1 Introduction

Current stellar evolution models predict that Galactic stars initially more massive than 25–30$`M_{}`$ lose more than 50% of their initial mass, and stars above 30–35$`M_{}`$ more than 80% (Maeder 1992, Woosley, Langer & Weaver 1993). Much of this mass loss is thought to occur during a short-lived and highly unstable stage preceding the Wolf-Rayet stage, that may be observed in the form of the luminous blue variables (LBVs; cf. Maeder 1989, Pasquali et al. 1997), which are located close to the Humphreys-Davidson (HD) limit in the HR diagram, beyond which no normal stars are observed (Humphreys & Davidson 1979). It is strongly debated what produces the giant LBV outbursts observed in these stars (Langer et al. 1994, Nota & Lamers 1997). Among the potential mechanisms (see Humphreys & Davidson 1994, and references therein), the idea that massive stars reach their Eddington luminosity close to the HD limit is particularly appealing since the Eddington limit appears to be located very close to the HD limit in the HR diagram (cf. Davidson 1971, Lamers & Noordhoek 1992). 

Time-varying massive star winds are able to produce circumstellar nebulae with kinematic properties similar to those observed around LBVs (Frank, Balick & Davidson 1995; Nota et al. 1995; García-Segura, Mac Low & Langer 1996, Frank, Ryu, & Davidson 1998). The morphology of the nebula around the extraordinary star $`\eta `$Carinae may give essential clues for the understanding of these giant outbursts. During its eruption from 1840 to 1860 A.D., $`\eta `$Car — today a telescopic object — was the second brightest star in the sky (van Genderen & The 1984). Recent observations by the Hubble Space Telescope have revealed in spectacular detail the resulting circumstellar nebula (Humphreys & Davidson 1994, Morse et al. 1998), the hourglass-shaped inner part of which is known as the Homunculus (Meaburn, Wolstencroft & Walsh 1987, Allen 1989). This is now the best studied example of the bipolar structures often observed around LBVs such as $`\eta `$Car (Nota et al. 1995).

## 2 Winds of rotating stars near the Eddington limit

It has been proposed (e.g. Maeder 1989) that sufficiently luminous stars, after core hydrogen exhaustion, may arrive at or exceed their Eddington limit 

$$\mathrm{\Gamma }L/L_{\mathrm{Edd}}=1,$$
(1)   
as their Eddington luminosity $`L_{\mathrm{Edd}}=4\pi cGM/\kappa `$ drops below their actual luminosity $`L`$ as the opacity coefficient $`\kappa `$ increases. The theory of radiation-driven stellar winds predicts that the mass loss rate will increase as $`\dot{M}(1\mathrm{\Gamma })^\mu `$ with $`\mu >0`$ for stars approaching the Eddington limit (cf. Castor, Abbott & Klein 1975). This mass loss rate formally diverges for $`\mathrm{\Gamma }1`$ (though see Owocki & Gayley 1997), suggesting that the strong mass loss associated with $`\mathrm{\Gamma }1`$ may be related to giant outbursts of LBVs.

However, rotation reduces the luminosity required for all external forces to balance each other at the stellar surface (Langer 1997). The Eddington limit $`\mathrm{\Gamma }<1`$ should then be replaced by a criterion that we will call the $`\mathrm{\Omega }`$-limit,

$$\mathrm{\Omega }v_{\mathrm{rot}}/v_{\mathrm{crit}}<1,$$
(2)   
with $`v_{\mathrm{crit}}^2v_{\mathrm{esc}}^2/2=GM(1\mathrm{\Gamma })/R`$, $`M`$ and $`R`$ being the stellar mass and radius, and $`v_{\mathrm{esc}}`$ being the polar escape velocity. $`\mathrm{\Omega }=1`$ implies that centrifugal and radiation force balance gravity at the equator, while at higher latitudes gravity still dominates. The feedback of rotation on the local surface luminosity is neglected here, since according to the generalized von Zeipel theorem (Kippenhahn 1977) the radiation flux at the equator may be either reduced or enhanced, depending on the internal rotation law; this may have been overlooked in the recent criticism of the $`\mathrm{\Omega }`$-limit by Glatzel (1998). Here, we apply the result of Friend & Abbott (1986) that the mass-loss rate of rotating hot stars depends on $`\mathrm{\Omega }`$ as $`\dot{M}(1\mathrm{\Omega })^\nu `$ with $`\nu 0.43`$.

We divide the evolution of a star approaching the $`\mathrm{\Omega }`$-limit — that is, going through an outburst — into three phases. In the first phase, before the star reaches the $`\mathrm{\Omega }`$-limit ($`\mathrm{\Omega }<1`$), it has the fast, energetic wind expected of a luminous blue star. In the second phase, reaching the $`\mathrm{\Omega }`$-limit ($`\mathrm{\Omega }1`$) has three consequences for the wind: the mass loss rate is much higher than before; the mass flux increases strongly at latitudes close to the equator; and the bulk of the wind is slow since the equatorial escape velocity is almost zero. In the third phase the outburst is over, the stellar radius has decreased, and the configuration is similar to the first phase. The smaller value of $`\mathrm{\Gamma }`$ before and after the outburst has two consequences: larger wind velocities due to larger escape velocities and smaller values of $`\mathrm{\Omega }`$ leading to more spherical winds.

To compute the latitudinal dependence of the wind properties of a star close to critical rotation ideally requires multi-dimensional models of the star and its outflowing atmosphere, which are not available. However, Langer (1997, 1998) argued that the stellar flux and the radius might still vary only weakly from pole to the equator in very luminous stars. Therefore, we applied equations similar to those found by Bjorkman & Cassinelli (1993, BC) for winds of rotating stars in the limit of large distance from the star: 

$$v_{\mathrm{}}(\theta )=\zeta v_{\mathrm{esc}}\left(1\mathrm{\Omega }\mathrm{sin}\theta \right)^\gamma ,$$
(3)   
$$(4\pi r^2\rho )_{\mathrm{}}(\theta )=\frac{\alpha }{2}\delta \dot{M}_0\left(1\mathrm{\Omega }\mathrm{sin}\theta \right)^\xi /v_{\mathrm{}}(\theta ),$$
(4)   
where we set the parameters defined in BC to $`\zeta =1`$ ,$`\gamma =0.35`$, and $`\xi =0.43`$. The correction factor $`\delta `$ is introduced to ensure that the total stellar mass loss rate $`\dot{M}_0`$ obeys $`\dot{M}_0=v_{\mathrm{}}(\theta )\rho (\theta )\mathrm{sin}\theta d\theta d\varphi `$ at the inner boundary of our grid (cf. Table 1). $`v_{\mathrm{}}`$ is the terminal wind velocity, and $`(4\pi r^2\rho )_{\mathrm{}}`$ the terminal wind density times $`4\pi r^2`$, as function of the polar angle $`\theta `$. The quantity $`\alpha `$ is defined by

$$\alpha =\left(\mathrm{cos}\varphi ^{}+\mathrm{cot}^2\theta \left(1+\gamma \frac{\mathrm{\Omega }\mathrm{sin}\theta }{1\mathrm{\Omega }\mathrm{sin}\theta }\right)\varphi ^{}\mathrm{sin}\varphi ^{}\right)^1,$$
(5)   
with $`\varphi ^{}=\mathrm{\Omega }\mathrm{sin}\theta v_{\mathrm{crit}}/(2\sqrt{2}v_{\mathrm{}}(\theta ))`$. Eq. (5) differs from the corresponding quantity defined by BC in their implicit formula (26) (cf. also Owocki et al. 1994, Ignace et al. 1996). This difference came along originally through a misinterpretation of BC’s equations, i.e., in equations (3) to (5) and in the equation defining $`\varphi ^{}`$, the $`\theta `$s were taken as $`\theta _0`$s, the initial co-latitude of the streamline. We note that since for wind compressed zone models near the equator it is $`\theta _0\theta `$, our models are similar to wind compressed zone models when $`\mathrm{\Omega }\mathrm{\Omega }_{\mathrm{th}}1`$, where $`\mathrm{\Omega }_{\mathrm{th}}`$ is the threshold value for the formation of a wind-compressed disk. Since $`\varphi ^{}<\pi /2`$ for $`\mathrm{\Omega }0.995`$, our formulation has the advantage of avoiding the formation of wind compressed disks for large $`\mathrm{\Omega }`$, a structure which can not be numerically resolved in our calculations, whose properties cannot be well predicted, and whose very existence has even been questioned (Owocki et al.1994, 1996). At the same time, wind density and velocity distributions obtained from our approach are similar to those derived from the formalism of BC, provided that $`\mathrm{\Omega }_{\mathrm{th}}1`$ and we choose $`\mathrm{\Omega }_{\mathrm{BC}}=\mathrm{\Omega }\mathrm{\Omega }_{\mathrm{th}}`$, where $`\mathrm{\Omega }_{\mathrm{BC}}`$ is the value to be inserted in BC’s equations. We shall see that the exact nature of the latitude dependence of the wind properties is not essential for our main results, as long as a dense wind with enhanced mass loss rate close to the equator occurs between two phases of an energetic, more or less spherical wind. 

We simulate the LBV outburst phenomenon by assuming the wind properties to be constant during each phase. For the pre-outburst wind, which is only used to initialize the numerical grid and to which our results are insensitive, we took a mass-loss rate of $`\dot{M}=10^3M_{}\mathrm{y}r^1`$, a wind final velocity $`v_{\mathrm{}}=450\text{ km\hspace{0.17em}s}\text{-1}`$ and $`\mathrm{\Omega }=0.53`$. For the post-outburst wind, we used two different sets of parameters. The one that we prefer appears to reasonably reproduce the observed morphology of the Homunculus, with $`\dot{M}=1.7\times 10^4M_{}\mathrm{y}r^1`$, $`v_{\mathrm{}}=1800\text{ km\hspace{0.17em}s}\text{-1}`$, and $`\mathrm{\Omega }=0.13`$. For comparison we also computed a model that represents $`\eta `$ Car’s presently observed wind parameters of $`\dot{M}=3\times 10^3M_{}\mathrm{y}r^1`$, $`v_{\mathrm{}}=800`$ km s<sup>-1</sup>, and $`\mathrm{\Omega }=0.3`$ (Davidson et al. 1995). The high wind velocity of our preferred model, $`v_{\mathrm{}}=1800\text{ km\hspace{0.17em}s}\text{-1}`$, corresponds to a stellar radius $`R21R_{}`$ or $`\mathrm{log}T_{\mathrm{eff}}4.7`$ for an O star wind with $`\zeta =3`$ at $`\mathrm{log}L/L_{}=6.4`$, implying that the star strongly contracted after the episode of mass-loss.

The parameters of the outburst wind largely determine the morphology of the resulting nebula. To compute these parameters, we assumed the following stellar properties: $`M=80M_{}`$, $`\mathrm{log}L/L_{}=6.4`$, and $`\mathrm{log}T_{\mathrm{eff}}=4.2`$, implying $`R=210R_{}`$, in agreement with observational estimates (cf. Humphreys & Davidson 1994). The final wind velocity follows from the escape velocity for a specified value of $`\mathrm{\Gamma }`$. For our preferred model we took a mass loss rate of $`\dot{M}=7\times 10^3M_{}\mathrm{y}r^1`$ and obtained a Homunculus mass of roughly $`0.15M_{}`$ (van Genderen & The 1984), while for our model assuming present-day wind parameters, we used an increased outburst mass loss rate of $`5\times 10^2M_{}\mathrm{y}r^1`$ such that we obtained a nebula mass of $`1M_{}`$ (Humphreys & Davidson 1994).

## 3 Hydrodynamic models

We perform two-dimensional hydrodynamic simulations of the wind interaction; first results have already been reported by García-Segura, Langer & Mac Low (1997). We use the hydrocode ZEUS-3D developed by M. L. Norman and the Laboratory for Computational Astrophysics. ZEUS-3D is a finite-difference, fully explicit, Eulerian code descended from the code described by Stone & Norman (1992). We used spherical coordinates for our simulations, with a symmetry axis at the pole, and reflecting boundary conditions at the equator and the polar axis. See García-Segura et al. (1996) for further details about our numerical method. Our models have grids of $`200\times 360`$ zones, with a radial extent of 0.125 pc, and an angular extent of $`90^{}`$. The innermost radial zone lies at $`r=9.7\times 10^{15}`$cm.

We compute the hydrodynamic evolution of the circumstellar gas, starting our computations at a time $`t=1840`$yr, and run them until $`t=1995`$yr. Our outburst scenario leads to a characteristic distribution of the circumstellar gas. The initial fast wind blows a stellar wind bubble which forms the background for the subsequent development of the nebula. During outburst, the wind becomes slow and dense, and the stellar rotation concentrates it toward the equatorial plane (BC; Ignace et al. 1996). When the final fast wind starts in the center of this nebula, it sweeps up the dense wind from the outburst into a thin, radiatively cooled shell that fragments due to dynamical instabilities (García-Segura et al. 1996). The shell expands more easily into the lower density wind at the poles, producing a double-lobed structure, as shown in Figures 1 and 2.

In Figure 1, we show six models computed with various values of $`\mathrm{\Omega }`$ and $`\mathrm{\Gamma }`$ for the outburst wind given in Table 1, and otherwise using the parameters of our preferred model. We find three major results: First, the nebular shape appears nearly independent of the Eddington factor, $`\mathrm{\Gamma }`$. Second, the nebula is strongly confined in the equatorial plane only for the case with nearly critical rotation ($`\mathrm{\Omega }=0.98`$). Finally, in this case, we obtain a structure very similar to that of the Homunculus, with two almost spherical lobes and an equatorial density enhancement that is an expanded relic of the outburst wind. In Figure 2 we show the results of a computation identical to that of Figure 1a — in particular with $`\mathrm{\Omega }=0.98`$ — but using twice the radial resolution. We find it striking how well this model, without much fine tuning, not only reproduces the large-scale, bipolar morphology, but also the small-scale turbulent structure seen in high-resolution observations of the Homunculus (Humphreys & Davidson 1994, Morse et al. 1998). 

We also compute a model at the same resolution as Figure 2 using $`\eta `$ Car’s presently observed wind parameters. We find that the large scale shape of the resulting nebula, shown in Figure 3, is almost identical to that shown in Figure 2, but the higher wind densities cause the wind termination shock to be strongly radiative. This changes the Vishniac instabilities (Vishniac 1983) seen in Figure 2 into ram-ram-pressure instabilities (Vishniac 1994, García-Segura et al. 1996), which have a much spikier morphology and a shorter wavelength, inconsistent with the observed structures. We emphasize that only the properties of the post-eruption wind are responsible for this feature, not the larger shell mass obtained in this case. This result might imply that the current wind is not representative of the wind over the last 140 years. Instead, $`\eta `$ Car might have been smaller and hotter in the recent past, with a post-outburst wind that has become slower and more dense with time. This is consistent with its gradual visual brightening over the last 140 years (Humphreys & Davidson 1994), and suggests that it is evolving towards another giant eruption.

## 4 Discussion

Our work extends previous hydrodynamic models for bipolar LBV nebulae (Nota et al. 1995, Frank et al. 1995, 1998, Dwarkadas & Balick 1998) by relating the outburst to the properties of evolving, massive, post-main sequence stars. In contrast to Nota et al. (1995) and Frank et al. (1995), who concluded that a strong equatorial density enhancement must have existed before the outburst occurred, we obtain the two lobes, including their small-scale structure, and the equatorial density enhancement self-consistently as a consequence of the evolutionary state of the star. Frank et al. (1998) used an arbitrary non-spherical wind during the post-outburst phase to produce a bipolar nebula. However, such a wind will only produce a bipolar shape if the wind termination shock is strongly radiative and therefore momentum conserving, a condition we have shown to be inconsistent with the small-scale morphology of the Homunculus. Dwarkadas & Balick (1998) introduce instead a ring-like density distribution, again without relating it to the underlying star.

Our result appears to be quite general, because all stars, even slow rotators, must by definition arrive at critical rotation if they approach their Eddington limit. The strong dependence of the nebula shape on $`\mathrm{\Omega }`$ shown in Figure 1, as well as the clear bipolar nature of virtually all LBV nebulae (Nota et al. 1995), lends strong support to the idea that LBV’s are stars approaching their Eddington limits that reach critical rotation and lose large amounts of mass quickly. In fact, a general mechanism for giant LBV outbursts is needed in order to understand the absence of stars beyond the Humphreys-Davidson limit in the HR diagram, and the bipolar nature of most LBV nebulae. Therefore, even though bipolar nebulae may also form from interacting binary stars (e.g. Han, Podsiadlowski & Eggleton 1995), and binarity has been repeatedly proposed also for $`\eta `$ Car (van Genderen, de Groot & The 1994, Damineli, Conti & Lopes 1997), it appears useful to continue to pursue single star models. The conjecture that the Homunculus nebula around $`\eta `$ Car is a paradigm rather than a freak has recently been supported by its strong similarity to the nebula around the LBV HR Carinae, as found by Weis et al. (1997) and Nota et al. (1997).

NL is grateful to Gloria Koenigsberger for enabling an extended visit to the UNAM Astronomy Institute, and to many colleagues there for their extremely warm hospitality. We thank J. Fliegner, S. Owocki, M. Peimbert, S. White and in particular J. Bjorkman for useful comments and discussions, and M. L. Norman and the Laboratory for Computational Astrophysics for the use of ZEUS-3D. The computations were performed at the Pittsburgh Supercomputing Center, the Supercomputer Center of the Universidad Nacional Autónoma de México and the Rechenzentrum Garching of the Max-Planck-Gesellschaft. This work was partially supported by the Deutsche Forschungsgemeinschaft through grants La 587/15-1 and 16-1 and by the US National Aeronautics and Space Administration.     

Figure Captions