Hydrodynamical simulations and similarity relations for eruptive mass loss from massive stars

TL;DR

This paper uses 1D hydrodynamical simulations to study eruptive mass loss in massive stars, revealing similarity relations that describe ejecta properties and their dependence on energy injection location, with implications for LBV outbursts and other stellar phenomena.

Contribution

It introduces a new similarity framework for modeling eruptive mass loss in massive stars, connecting simulation results with analytic scaling relations.

Findings

Ejecta velocity and density follow exponential similarity forms.

Total ejecta mass and energy scale with energy injection parameters.

Bulk ejecta speed is comparable to the star's surface escape velocity.

Abstract

Motivated by the eruptive mass loss inferred from Luminous Blue Variable (LBV) stars, we present 1D hydrodynamical simulations of the response from sudden energy injection into the interior of a very massive ($100 \, M_\odot$) star. For a fiducial case with total energy addition set to a factor $f=0.5$ of the net stellar binding energy, and applied within the stellar envelope, we detail the dynamical response that leads to ejection of the outermost $7.2 \, M_\odot$. We find that the ejecta's variations in time $t$ and radius $r$ for the velocity $v$, density $\rho$, and temperature $T$ are quite well fit by similarity forms in the variable $r/t \approx v$. Specifically the scaled density follows a simple exponential decline $\rho t^{3} \sim \exp (-r/v_{\rm o} t)$. This `exponential similarity' leads to analytic scaling relations for total ejecta mass $\Delta M$ and kinetic energy $\Delta K$ that agree well with the hydrodynamical simulations, with the specific-energy-averaged speed related to the exponential scale speed $v_{\rm o}$ through ${\bar v} \equiv \sqrt{2 \Delta K/\Delta M} = \sqrt{12} \, v_{\rm o}$, and a value comparable to the star's surface escape speed, $v_{\rm esc}$. Models with energy added in the core develop a surface shock breakout that propels an initial, higher-speed ejecta ($>$5000km s$^{-1}$), but the bulk of the ejected material still follows the same exponential similarity scalings with ${\bar v} \approx v_{\rm esc}$. A broader parameter study examines how the ejected mass and energy depends on the energy-addition factor $f$, for three distinct model series that locate the added energy in either the core, envelope, or near-surface. We conclude by discussing the relevance of these results for understanding LBV outbursts and other eruptive phenomena, such as failed supernovae and pulsational pair instability events.

Figures

Tables

Table 1. Table 1: Input parameters and results for fiducial model of energy addition to stellar envelope.

Input parameters:
initial mass	$M_{\ast }$	100 $M_{\odot }$
initial radius	$R_{\ast }$	35 $R_{\odot }$
surface escape speed	$v_{\mathrm{esc}}$	1044 km s-1
gravitational binding energy	$E_g$	$-3.31\times {10}^{51}$ erg
internal energy	$E_i$	$+2.65\times {10}^{51}$ erg
nets initial binding energy	$E_{\ast }=E_g+E_i$	$-6.65\times {10}^{50}$ erg
energy addition factor	$f$	$0.5$
added energy	$-f{E}_{\ast }$	$+3.32\times {10}^{50}$ erg
radius range for added energy	$r$	7-35 $R_{\odot }$
mass range for added energy	$M_r$	$>76M_{\odot }$
Results:
ejecta mass	$\Delta M$	7.2 $M_{\odot }$
ejecta kinetic energy	$\Delta K$	$0.93\times {10}^{50}$ erg
energy-averaged speed	$\overline{v}$	1125 km s-1
similarity speed	$v_o$	325 km s-1
similarity density	${\rho }_o{t}_o^3$	$1.65\times {10}^{10}$ g cm-3 s3
similarity temperature	$T_o{t}_o^{1.08}$	$2.6\times {10}^{11}$ K s1.08

Table 2. Table 2: For the solid curves fits to the mass data in figure 12 , parameters for the fitting function ( 28 ), for energy addition near the surface, in the envelope, and in the core.

Heating region	$M_{\mathrm{bot}}$	$f_b$	$n$
$M_r>89M_{\odot }$	$11M_{\odot }$	0.55	1.25
$M_r>76M_{\odot }$	$24M_{\odot }$	1.00	1.75
	$100M_{\odot }$	1.70	3.50

Equations

$$\frac{\partial \rho }{\partial t}+\frac{1}{r^2}\frac{\partial (\rho v{r}^2)}{\partial r}=0,$$

$$\frac{Dv}{Dt}=\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial r}=-\frac{1}{\rho }\frac{\partial P}{\partial r}-\frac{G{M}_r}{r^2},$$

$$P=P_{\mathrm{rad}}+P_{\mathrm{gas}}=\frac{a_{\mathrm{rad}}{T}^4}{3}+\frac{\rho kT}{\mu },$$

$$\frac{\partial (\rho v^2/2+E)}{\partial t}+\frac{1}{r^2}\frac{\partial [(\rho v^2/2+E+P)v{r}^2]}{\partial r}=-\rho v\frac{G{M}_r}{r^2},$$

$$\frac{DE}{Dt}=-\left(E+P\right)\nabla \cdot \mathbf{v}=-\frac{E+P}{r^2}\frac{\partial (v{r}^2)}{\partial r}.$$

$$B_{\mathrm{n}}\equiv \frac{v^2}{2}-\Phi +h,$$

$$\Phi \equiv \frac{G{M}_r}{r}+4\pi G{\int }_r^{\infty }{r}^{\prime }\rho (r^{\prime })\text{d}{r}^{\prime },$$

$$h\equiv \frac{E+P}{\rho }=\frac{5}{2}\frac{kT}{\mu }+\frac{4a_{\mathrm{rad}}{T}^4}{3\rho }$$

$$\frac{D{B}_{\mathrm{n}}}{Dt}=\frac{1}{\rho }\frac{\partial P}{\partial t}-\frac{\partial \Phi }{\partial t},$$

$$\rho (r,t)\approx \frac{{\rho }_{\mathrm{o}}{t}_{\mathrm{o}}^3}{t^3}{e}^{-r/v_{\mathrm{o}}t},$$

$$\Delta M(r,t)$$

$$\Delta M\equiv 8\pi {\rho }_{\mathrm{o}}{t}_{\mathrm{o}}^3{v}_{\mathrm{o}}^3=7.16{\mathrm{M}}_{\odot },$$

$$m(x)\equiv \left(1+x+\frac{x^2}{2}\right)e^{-x}=\frac{M_{\ast }-M_r}{\Delta M},$$

$$v(M_r)=v_{\mathrm{o}}\chi \left(\frac{M_{\ast }-M_r}{\Delta M}\right).$$

$$\Delta K$$

$$\bar{v}\equiv \sqrt{\frac{2\Delta K}{\Delta M}}=\sqrt{12}{v}_{\mathrm{o}}=1125\mathrm{km}{\mathrm{s}}^{-1}\approx v_{\mathrm{esc}}.$$

$$E_{\mathrm{fin}}\approx (1-f)E_{\ast }-\Delta K=-4.2\times 10^{50}\mathrm{erg}.$$

$$\rho t^3=d_{\mathrm{p}}{\left(\frac{1000\mathrm{km}{\mathrm{s}}^{-1}}{v}\right)}^9;v>v_{\mathrm{b}},$$

$$\Delta M_{\mathrm{tail}}$$

$$m(x,x_{\mathrm{b}})$$

$$\Delta K\approx (1-f)E_{\ast }-E_{\mathrm{fin}},$$

$$\frac{DE}{Dt}=-3\frac{E+P}{t}.$$

$$\delta \equiv -\frac{D\ln E}{D\ln t}=3\left(1+\frac{P}{E}\right)=3+\frac{1}{1-\beta /2}.$$

$$E(r,t)=\frac{E_{\mathrm{o}}{t}_{\mathrm{o}}^{4.33}}{t^{4.33}}{e}^{-r/v_{\mathrm{o}}t},$$

$$E=3P_{\mathrm{rad}}+\frac{3}{2}{P}_{\mathrm{gas}}\approx \frac{9}{2}{P}_{\mathrm{rad}}=\frac{3}{2}{a}_{\mathrm{rad}}{T}^4.$$

$$T(r,t)\approx \frac{T_{\mathrm{o}}{t}_{\mathrm{o}}^{1.08}}{t^{1.08}}{e}^{-r/4v_{\mathrm{o}}t},$$

$$T_{\mathrm{o}}{t}_{\mathrm{o}}^{1.08}={\left(\frac{E_{\mathrm{o}}{t}_{\mathrm{o}}^{4.33}}{\frac{3}{2}{a}_{\mathrm{rad}}}\right)}^{1/4}\approx 2.7\times 10^{11}\mathrm{K}{\mathrm{s}}^{1.08}.$$

$$\frac{\Delta M}{M_{\mathrm{bot}}}(f)={\left[1+{\left(\frac{f_{\mathrm{b}}}{f}\right)}^{2n}\right]}^{-1/2},$$

Full text

1 Introduction and Background

Hot, luminous, massive stars lose mass through strong stellar winds, driven by line-scattering of the star’s continuum radiation (Puls et al., 2008). But for many such massive stars, the overall mass loss may actually be dominated by transient mass ejection episodes associated with eruptive Luminous Blue Variable (eLBV) stars (Smith & Owocki, 2006). The most prominent and extensively studied eLBV is the very luminous ($L\approx 5\times 10^{6}{\rm L}_{\odot}$), very massive ($M>100{\rm M}_{\odot}$) star $\eta$ Carinae, which during the “Great Eruption” of the 1840’s is inferred (Smith et al., 2003) to have ejected $>10{\rm M}_{\odot}$, resulting in the famous Homunculus nebula observed today. The associated mass loss rate of $>1{\rm M}_{\odot}$ yr*-1* is well above what can be explained by a line-driven wind ($\lesssim 10^{-4.5}{\rm M}_{\odot}$ yr*-1*, Smith & Owocki, 2006); but it could in principle be achieved by a quasi-steady continuum-driven wind associated with the star’s super-Eddington luminosity during this epoch (Shaviv, 2000; Owocki et al., 2004, 2017; van Marle et al., 2008; Quataert et al., 2016).

However, a recent extensive empirical study by Smith et al. (2018) argues that many features of the Great Eruption and the resulting nebula can not be explained by such a quasi-steady, super-Eddington wind alone, but must also include a relatively sudden explosive component more analogous to a supernova (SN) explosion than a steady wind. Indeed, much as inferred for the class of pre-SN eLBV’s (see review by Smith, 2014, and references therein), the explosive mass loss likely interacts with a more prolonged, quasi-steady, enhanced wind outflow prior to the eruption. Such pre-SN eruptions may be triggered by core instabilities in the final stages before core collapse (Chen et al., 2014), which lead to energy deposition into the envelope, e.g. by gravity waves (Piro, 2011; Quataert & Shiode, 2012; Fuller & Ro, 2018), or into the core itself, e.g. by pair-instability pulsation (Barkat et al., 1967; Kasen et al., 2011; Yoshida et al., 2016; Takahashi et al., 2016; Woosley et al., 2007; Woosley & Heger, 2015; Woosley, 2017), or by late-shell-burning instabilities (Smith & Arnett, 2014). Recent studies by Gieles et al. (2018) of the formation of supermassive stars by multiple mergers in the very dense cores of globular clusters also posit a strong mass loss that regulates the maximum mass to be of order 1000 ${\rm M}_{\odot}$.

For $\eta$ Carinae, the bipolar form of the Homunculus, along with its irregular equatorial “skirt”, has led to broad conjecture (Podsiadlowski, 2013; Portegies Zwart & van den Heuvel, 2016; Smith et al., 2018) that the eruption might have stemmed from the merger of two massive stars (Gallagher, 1989; Morris & Podsiadlowski, 2006, 2007, 2009; MacLeod et al., 2018). In such a merger scenario, the orbital angular momentum leads during a common envelope phase (Ivanova et al., 2013; Podsiadlowski, 2013; Justham et al., 2014) to mechanical ejection of the equatorial skirt, while the bipolar Homunculus stems from the energy deposited into the envelope from the merging of the stellar cores. This energy addition leads to a prompt ejection of mass, followed perhaps by a strong quasi-steady wind driven by the super-Eddington luminosity of the post-ejection star. While both explosive SNe and steady stellar winds have been quite extensively studied through detailed dynamical models, there have been relatively limited attempts to study the dynamics of such sudden mass eruptions (see, however, Dessart et al., 2010; Ro & Matzner, 2013, 2017).

A full simulation of such a stellar merger should be 3D, to account for the inspiralling of the two stars, the associated spin-up of the common envelope, and eventual merging of the stellar cores. This could naturally lead to a distinctly non-spherical form of the ejected nebula, including perhaps both a thin equatorial ejecta along with an associated pinching of the mass ejection into a bipolar form, much as seen in the skirt plus Homunculus of $\eta$ Carinae. But to build insight for such complex and computationally expensive 3D models, we elect here to focus first on just the effect of the energy addition associated with the merger of the stellar cores, within the context of idealized 1D, spherically symmetric models.

Specifically, the paper here carries out a series of 1D time-dependent hydrodynamical simulations of the disruption that results when a generic energy addition of unspecified origin is introduced impulsively into a massive star and core. A key feature is that, unlike supernova explosions, the level of energy addition here is taken to be only an order unity factor $f$ of the stellar binding energy ${E_{\ast}}$, thus leading to only partial disruption of the stellar envelope. While some mass fraction is ejected to full escape from the star, the bulk of the star’s mass remains bound by the stellar gravity, and so falls back onto a somewhat inflated surviving star (somewhat analogous to what is thought to occur in failed explosions from core collapse, e.g. Lovegrove & Woosley, 2013; Coughlin et al., 2018). This is intended as just a first step toward future models grounded in specific mechanisms for energy addition, and in the case of merger models, including 2D and 3D dynamics that account explicitly for the role of angular momentum transport in mass ejection and envelope spin-up.

Section 2 describes the basic method and setup for our fiducial model of a $M_{\ast}=100{\rm M}_{\odot}$ star with impulsive addition of a fraction $f=1/2$ of its binding energy, deposited over the stellar envelope above the core. The key results in section 3 detail how this unbinds the outermost $\sim 7.2\%$ of the star’s mass into an expanding ejecta, with the rest remaining gravitationally bound and thus falling back onto an inflated stellar envelope. The analysis in section 4 then shows that the long-term evolution of this ejected mass can be well described by relatively simple, analytic similarity forms for the velocity, density and temperature. Section 5 explores the effect of adding energy within (vs. above) the core, for cases with both half and double the stellar binding energy. Section 6 extends this to a broad parameter study for how the mass and kinetic energy of the ejecta scales with energy addition factor $f$ for core, envelope and near-surface locations of the energy addition. We conclude (section 7) with a discussion of the implications of our results for understanding observed LBV’s, and an outline for future work.

2 1D Simulation of Envelope Ejection from Sudden Energy Addition

2.1 Conservation equations

The simulations here use a one-dimensional (1D), spherically symmetric implementation of the multi-dimensional Eulerian hydrodynamics code described by Hirai et al. (2016). In its full multi-D form a key distinguishing feature of this code lies in its fast and accurate gravity solver; but in the 1D form used here, the gravity is trivially solved in spherical symmetry. The code is otherwise quite similar to other hydro codes (e.g. PLUTO, Mignone et al., 2007) that use Godunov-type schemes to integrate finite-volume forms of the spherical equations for conservation of mass, momentum and energy, and thereby follow the time ($t$) and radius ($r$) evolution of mass density $\rho$, radial velocity $v$, and temperature $T$. Conservation of mass takes the form,

[TABLE]

while momentum conservation yields an equation of motion of the form,

[TABLE]

with $M_{r}$ the mass contained within radius $r$, and $G$ the gravitation constant.

Here $\partial/\partial t$ and $D/Dt\equiv\partial/\partial t+{\bf v}\cdot\nabla$ are respectively the (Eulerian) partial derivative at a fixed radius, and the total (Lagrangian) time variation in a frame co-moving with the fluid. In our 1D simulations, the latter represents the time variation for a fixed mass coordinate $M_{r}$. The total pressure includes components from both radiation and gas,

[TABLE]

with $k$ Boltzmann’s constant, $\mu$ the mean molecular weight (in g), and $a_{\rm rad}$ the radiation constant. The mean molecular weight distribution over mass coordinate is fixed throughout the simulation, with no mixing assumed.

For energy conservation, we evolve the adiabatic form

[TABLE]

where the total internal energy again includes contributions from both radiation and a monatomic ideal gas111For a monatomic gas the internal energy density is $\frac{1}{2}P_{\rm gas}$ for each of the $n=3$ degrees of freedom, which also sets the adiabatic index $\gamma=(n+2)/n=5/3$., $E=3P_{\rm rad}+\frac{3}{2}P_{\rm gas}$.

This adiabatic form ignores radiative diffusion and transport, since in the very optically thick flows that develop here this is much smaller than competing terms from energy advection. Alternatively, we can use the mass and momentum eqns. (1) and (2) to recast this energy conservation (4) into a form for just the internal energy $E$,

[TABLE]

This form proves useful for interpreting the computed evolution of internal energy and temperature (see Appendix A).

2.2 Fiducial model and numerical specifications

To explore the nature of mass ejection in eLBV stars like $\eta$ Carinae, let us first focus on a fiducial massive-star model with the parameters given in Table 1. This starts with a $M_{\ast}=100{\rm M}_{\odot}$ stellar model evolved with the public stellar evolution code MESA (v10398; Paxton et al., 2011, 2013, 2015, 2018) up to the time that the central hydrogen mass fraction has dropped to $X=0.2$. The resulting He-enriched core contains about 76% of the star’s mass, but is only about 20% of its $R_{\ast}\approx 35{\rm R}_{\odot}$ surface radius. Because of the prominence of radiation pressure, the star’s total internal energy, $E_{\rm i}=2.65\times 10^{51}$ erg, is a substantial fraction (more than the 50% for virial equilibrium of a pure ideal gas) of the gravitational binding energy, $-E_{\rm g}=3.31\times 10^{51}$ erg; the associated net binding energy of the whole star is thus ${E_{\ast}}=-6.65\times 10^{50}$ erg. (See Table 1.)

sOur approach idealizes the complex merger processes by assuming that the structure of the merger product is somewhat similar to that of this fiducial model, with a combined mass similar to $\eta$ Carinae. The central hydrogen fraction corresponds to a late case A merger, but we do not wish here to specify the details of the pre-merger binary system, or of the merger process itself. Instead, within this assumed 1D structure, we focus only on one key aspect of a stellar merger, namely the energy addition to the merged envelope. The actual amount and distribution of the additional energy will depend on the details of the merger process. We do not aim here to reproduce this, but rather explore the overall effects of varying the location of energy deposition from the core, to the envelope, to the near-surface.

The energy directly available from a merger is given by the sum of the binding energies of the two stars, plus their orbital energy, minus the binding energy of the assumed 100${\rm M}_{\odot}$ star. Neglecting for simplicity any net change in binding energy as relatively small, the orbital kinetic energy scales with the product of masses divided by the final separation $a$ at the onset of the dynamical stage, $E_{\rm kin}\approx GM_{1}M_{2}/2a$. For the optimal case of two equal-mass stars with $M_{1}=M_{2}=50M_{\odot}$, and taking a final separation in the range $a=(7-35)R_{\odot}$ that corresponds to the core to surface radius of our merged star, this orbital energy deposition lies in the range $E_{\rm kin}\approx(1.3-6.5)\times 10^{50}$ erg, which is an order-unity factor $f$ of the net stellar binding energy of our assumed model. In principle, even more energy could arise from the difference in the binding energies, or even from sudden hydrogen burning if original envelope material is mixed into the hot, merged core (Podsiadlowski et al., 2010).

In our simulations we first map the stellar model on the centre of a spherical grid, with a dilute atmosphere of negligible dynamical influence placed around the star. At an initial time we impulsively add some fraction $f$ of the stellar binding energy to the stellar envelope above this core, distributed in proportion to the mass. The analysis below focuses on detailed results for a fiducial $f=0.5$ model, i.e. with an added energy $\Delta E\approx+3.3\times 10^{50}$ erg; but, as discussed in section 6, we have also computed a full grid of models from $f=0.1$ to $2.0$ in steps of $0.1$, with three model series with energy deposition in either the core, envelope, or near surface. As noted in the introduction, while motivated by the specific notion of a merger for the giant eruption in $\eta$ Carinae, this choice of a broad range for the energy fraction and location mean that results could also be relevant for eruptive mass loss from a variety of specific mechanisms.

To follow the evolution of mass ejected well beyond the initial stellar radius $R_{\ast}\approx 35{\rm R}_{\odot}$, our numerical grid extends from the stellar center $r=0$ to a large maximum radius $r=R_{\rm max}=4000{\rm R}_{\odot}\gg R_{\ast}$, using $nr=25,000$ radial cells with spacing that increases with grid index $i$ as $\Delta r_{i+1}=1.000186\,\Delta r_{i}$, with $\Delta r_{0}=0.0072\,{\rm R}_{\odot}$. We apply outgoing boundary conditions at the outer boundary radius $R_{\rm max}=4000{\rm R}_{\odot}$. In order to follow mass outflow to this large radius, our fiducial model is evolved to a very long time, $t_{\rm max}=4050$ ks, with time steps set to a maximum of $0.9$ times the Courant time (Courant et al., 1967) .

3 Simulation Results for Fiducial Model

3.1 Mass ejection vs. fallback

Figures 1 illustrate the radial evolution of mass shells over the initial $300$ ks. The impulsive introduction of energy at initial time $t=0$ is applied over the radial range ($r=7-35\,{\rm R}_{\odot}$), distributed in proportion to the mass. The mass contours show that the initial response of nearly all the mass – even in the core with $M_{r}<76\,{\rm M}_{\odot}$ – is to expand in radius. For most of the mass, this expansion is followed by a partial fallback to a final radius that remains somewhat above the initial radius.

The outer $\sim 7.2{\rm M}_{\odot}$ of the envelope escapes completely from the star. The right panel plots contours of mass measured from $r$ outward, i.e. $M_{\rm ro}\equiv M_{\ast}-M_{r}$, using now a linear radius scale from $r=0-700{\rm R}_{\odot}$, and time plotted over the full $t=0-4$ Ms of the simulation. This shows that a parcel with $M_{\rm ro}=7.5{\rm M}_{\odot}$ rises to $r>500{\rm R}_{\odot}$ before falling back; indeed, we shall see that even the still-rising parcels with $M_{\rm ro}=7.3$ and 7.4 ${\rm M}_{\odot}$ have a slightly negative net total energy, and so will eventually fall back. But parcels with $M_{\rm ro}\lesssim 7.2{\rm M}_{\odot}$ ($7.2$% of the stellar mass) do completely escape.

3.2 Variation of Bernoulli energy with mass and time

To explore further how this mass ejection is achieved, the left panel of figure 2 shows a contour plot of velocity $v$ vs. mass $M_{r}$ over the early time interval $t=0-50$ ks. As denoted by the colorbar, the colors from red to blue indicate positive (outflowing) velocities in the range $v=0-1900$ km/s; the black in the lower right corner denotes mass and time where $v<0$, i.e. infall.

These velocity contours clearly show two propagating discontinuities, which we identify with direct and rebound shocks that arise from the sudden initial deposition of energy in the mass range $M_{r}=76-100\,{\rm M}_{\odot}$.

The right panel of figure 2 shows a corresponding contour plot of the total specific energy, also known as the Bernoulli energy,

[TABLE]

where the gravitational potential is defined by

[TABLE]

and

[TABLE]

is the total specific enthalpy due to both gas and radiation. This is defined in terms of the associated total internal energy $E$ and total pressure $P$, and given in the second equality in terms of density $\rho$, temperature $T$, and associated constants. The colorbar shows the numerical value of $B_{\rm n}$ (in units of $G{\rm M}_{\odot}/{\rm R}_{\odot}$), with white to blue being positive, and yellow to red being negative.

The red boldface contour more explicitly marks the divide between the upper-right region in time and mass, where this total energy is positive ($B_{\rm n}>0$), vs. the lower regions where it is negative ($B_{\rm n}<0$). This roughly separates mass parcels that will escape the star from those that will fall back to it.

For a set of mass parcels with mass coordinate in the range $M_{r}=89-99{\rm M}_{\odot}$ (as marked in the right legend), figure 3 plots the corresponding Lagrangian variations in time $t$ of this Bernoulli energy $B_{\rm n}(t,M_{r})$. Note that following the boosts from the two shocks, $B_{\rm n}$ for each parcel remains nearly constant, but with just a slow decline for bound parcels.

The reason for this decline can be understood through analysis of the time-dependent equations for conservation of mass and energy, eqns. (1) and (4). With some manipulation, these can be recast to show that

[TABLE]

where again $D/Dt$ represents the total Lagrangian variation co-moving with the flow.

In a steady-state flow with all $\partial/\partial t=0$, this means $B_{\rm n}=$ constant along the flow of a mass parcel. But if the right-hand-side of eqn.  (9) is negative, then $B_{\rm n}$ will decline in time within a mass parcel, as seen in figure 3 for the bound parcels. At later times, these intrinsic time variations vanish, and so the Bernoulli energy over the full range of our model, extending to the maximum radius $R_{\rm max}=4050{\rm R}_{\odot}$ and final time $t_{\rm fin}=4050$ ks, does indeed approach a fixed constant value for each mass parcel of the ejecta. Specifically, solving $B_{\rm n}(\Delta M,t_{\rm fin})\equiv 0$ from our simulation gives $\Delta M=7.18{\rm M}_{\odot}$ as an estimate for the total ejecta mass in this fiducial model.

4 Similarity relations

4.1 Free-expansion velocity of ejecta

Over this full range in radius $r=0-4000{\rm R}_{\odot}$ and time $t=0-4.05$ Ms, the top panel of figure 4 plots contours of the velocity $v$. The resulting overall fan-like form of these contours suggests that this velocity variation can be approximated by a similarity form $v(r,t)\approx r/t$. The bottom panel of figure 4 shows line plots of $v$ vs. $r/t$ for the selected time snapshots denoted in the righthand legend. Note that indeed, for all but the very early time $t=50$ ks, the velocity approaches the similarity form $v(r,t)=r/t$ marked by the dashed line.

This similarity form for velocity is just what is found in free-expansion, such as occurs in supernova explosions, or even in the Hubble law for the expanding universe. However, it is important to emphasize that in this case it applies only approximately for unbound material that escapes from the star, and then only at large $r/t$, and for times $t\gg 100$ ks.

4.2 Exponential similarity relation for density

Since in such a free expansion, the specific volume increases as $r^{3}\sim t^{3}$, one expects the density to decline as $\rho\sim 1/t^{3}$. Figure 5 thus shows plots of the time-scaled density $\rho t^{3}$ vs. the similarity variable $r/t$. Note that this scaled density is not constant, but declines with $r/t$, approaching the dashed line at large $r/t$ and larger times $t$. On this semi-log scale, such a linear decrease indicates an exponential decline, represented here by the dashed line for the quoted fit function,

[TABLE]

where we find the visually best fit parameters $\rho_{\rm o}t_{\rm o}^{3}\approx 1.65\times 10^{10}$ g cm*-3* s3 and $v_{\rm o}\approx 325$ km/s. Such an exponential similarity relation for density has also been noted previously for SNIa explosions (Dwarkadas & Chevalier, 1998), but in the present context of eruptions we argue below (sections 5 and 6) that this may be a signature of the role of gravity in slowing the ejecta and keeping some of the stellar mass bound.s Appendix A shows that analogous exponential similarity relations hold here for the internal energy and temperature.

4.3 Mass distribution of ejecta velocity

We can use the density fitting formula (10) to derive an associated estimate for the time and space variation of ejected mass,

[TABLE]

where the total ejected mass is

[TABLE]

with the latter evaluation applying the above fitting values for $\rho_{\rm o}t_{\rm o}^{3}$ and $v_{\rm o}$. The resulting total mass agrees very well with the net mass ejection $7.2{\rm M}_{\odot}$ inferred from the numerical simulation, as indicated by figure 1b. The spatial and temporal distribution of the ejected mass also follows a similarity form, given in terms of the dimensionless similarity variable $x=r/v_{\rm o}t$ by

[TABLE]

where the latter expression gives the relationship to the standard mass coordinate $M_{r}$. This represents a cumulative distribution function for the fraction of mass ejected above a scaled velocity $x=v/v_{\rm o}$.

If we next apply the similarity form for the velocity $v=r/t$, we see that $m(v/v_{\rm o})$ also gives the mass dependence on ejecta velocity. Defining then the inverse function $\chi(m)\equiv m^{-1}(\chi)$, we can write a similarity form for the mass distribution of ejecta velocity,

[TABLE]

Figure 6 compares plots for the velocity vs. mass for the indicated time snapshots of the our full hydrodynamical simulation; the dashed curve represents the associated analytic form $v(M_{r})$ from eqn.  (14). For large times and masses above the marginally unbound mass with $M_{r}=M_{\ast}-\Delta M=92.8{\rm M}_{\odot}$, there is again remarkably close agreement between the analytic similarity relation and the numerical simulation results. Note that only a small amount of mass reaches the highest speeds $v>3000$ km s*-1*.

4.4 Total kinetic energy of ejecta

We can also use these similarity relations to obtain an analytic scaling for the total kinetic energy of the ejected material,

[TABLE]

where the final numerical evaluation is again in quite good agreement with what’s found in the full numerical simulations. It is also comparable to the kinetic energy inferred for the ejecta in $\eta$ Carina (Smith et al., 2003).

The scaling in the last line of (15) shows that the characteristic speed $v_{\rm o}$ from the exponential truncation of density is the key parameter setting the specific kinetic energy of the ejecta, $\Delta K/\Delta M$. The associated averaged speed is

[TABLE]

We show in section 6.1 (see figure 9b) that the connection in the last equality between the average ejecta speed and the star’s surface escape speed is a quite general property for eruption models with moderate energy factors $f<1$.

Since the ejecta is gravitationally unbound, and generally also highly supersonic, the kinetic energy essentially represents the total energy of the ejecta. From total energy conservation between the initial and final state, we thus see that the final binding energy of the remaining star just depends on the kinetic energy and the net initial energy,

[TABLE]

As this is less negative than the binding energy ${E_{\ast}}=-6.65\times 10^{50}$ erg of the original, unperturbed star, it indicates that the post-eruption stellar envelope remains somewhat inflated compared to a star in full thermal equilibrium. If this gives a sufficiently enhanced luminosity, it could lead to an extended epoch of further mass loss through a quasi-steady, super-Eddington, continuum-driven wind (Quataert et al., 2016; Owocki et al., 2017). Detailed consideration of this will be left to future work.

5 Energy Addition to Core

5.1 Adding half the star’s binding energy to its core

To explore how the ejection of mass depends on the location of the added energy, let us next consider a case with the same addition factor $f=0.5$, but now with the added energy distributed evenly in mass over the stellar core, with $M_{r}<76{\rm M}_{\odot}$. Figure 8 shows key results for such a model.

The top panel plots contours of the mass coordinate $M_{r}$ vs. log radius and time over the full range out to the final simulation time, here just $600$ ks. Note that all the mass below $M_{r}<98{\rm M}_{\odot}$ now remains bound, but is excited into a pulsation with a gradual subsequent damping.

The middle panel shows that the Bernoulli energy undergoes a sharp jump from a shock that reaches from the core to near-surface after just $\sim 13$ ks. The red contour for $B_{\rm n}=0$ delineates the unbound mass above from the bound (oscillating) mass below. Extensions to longer time shows that only about $\Delta M\approx 1.7{\rm M}_{\odot}$ attains the positive net energy needed to escape, which is substantially lower than the $7.2{\rm M}_{\odot}$ lost in the standard $f=0.5$ case of energy addition to the stellar envelope.

The bottom panel plots the scaled density vs. similarity speed $r/t$. Note that this now extends out to much higher speeds ($>$5000 km s*-1*) than in the fiducial case of envelope heating. Indeed, at speeds above a “breakout” speed $v_{\rm b}=3500$ km/s, the density no longer fits the usual exponential decline (shown by the dashed line in this semi-log plot), but rather follows a steep power-law (dotted curve) defined by

[TABLE]

where the visual best-fit gives $d_{\rm p}\approx 2.5\times 10^{9}$ g cm*-3* s3. Integrating this similarity form from the breakout speed $v_{\rm b}$ shows that the associated mass in this power-law tail is actually very small,

[TABLE]

This is only about $0.08$% of the total ejecta mass $\Delta M\approx 1.7{\rm M}_{\odot}$ (as estimated from the mass with a positive total energy at the final time $t_{\rm fin}=600$ ks). Analogous integration for the kinetic energy shows that this power-law tail also represents only a small fraction, $\sim$2%, of the total kinetic energy of the ejecta, $\Delta K\approx 1.9\times 10^{49}$ erg. The mean speed of the ejecta is ${\bar{v}}=1055$ km s*-1*, again very near the stellar escape speed, $v_{\rm esc}=1044$ km s*-1*.

The cumulative distribution of ejected mass fraction vs. scaled speed, defined as $m(x)$ in eqn. (13) for the case of a purely exponential similarity, can be readily generalized to take account of this power-law tail for $x>x_{\rm b}\equiv v_{\rm b}/v_{\rm o}$,

[TABLE]

In summary, such a core heating model leads to initial, high-speed ejecta that has only a small fraction of the total ejected mass, which itself is much less than found in the standard envelope heating model with the same energy addition factor $f=0.5$.

5.2 Adding twice the star’s binding energy to its core

To explore potential links to core-collapse SNe, let us next consider a model in which the energy added to the core is twice the binding energy (i.e., $f=2.0$), and so could in principle completely disrupt the entire star.

Figure 8 shows results. As in figure 8, the middle and bottom panels again show respectively the evolution of the Bernoulli energy, and the variation of scaled density with the similarity speed $r/t$. In the former, the red contour for $B_{\rm n}=0$ shows that the parcels deep into the stellar core are developing a positive energy, which eventually leads to ejection of more than $87{\rm M}_{\odot}$. The latter again shows a strong power-law tail (dotted curve) at high similarity speed (now even exceeding 10,000 km s*-1*), with the same form as eqn. (18), but now with a lower transition breakout speed, $v_{\rm b}=2750$ km s*-1*. This and the higher best-fit normalization ($d_{\rm p}\approx 1.34\times 10^{11}$ g cm*-3* s3) gives a higher tail mass $\Delta M_{\rm tail}\approx 0.33{\rm M}_{\odot}$ and kinetic energy $\Delta K_{\rm tail}\approx 7.4\times 10^{49}$ erg. The corresponding fractions of total ejecta mass ($\sim$0.4%) and kinetic energy ($\sim$11%) are thus greater than for the $f=0.5$ case, but still quite modest.

Indeed, it seems that the bulk of the mass ejection follows closely the similarity solution based on the exponential density decline in similarity speed. The top panel of figure 8 shows the mass distribution of the speed $v$ at the cited times. There is quite remarkably good agreement with the analytic form (14) (dashed curve), apart from the slight kink in the simulation results at the mass $M_{r}\approx 76{\rm M}_{\odot}$ near the core-envelope boundary (marked by arrows in figure 8). This suggests that even the marked differences in stellar structure between core and envelope have only a relatively minor effect on the density and speed of the mass ejecta. Thus the similarity forms based on the exponential scaling $\exp(-r/v_{\rm o}t)$ for density may apply quite robustly to eruptive mass loss from both core and envelope heating with energy addition factors of order unity. For such cases, gravity still plays a key role in regulating the mass ejecta properties, for example setting the mean speed to be of order the surface escape speed.

But these models with energy added in the core also lead to a prompt, high-speed ejecta from the sudden shock breakout through the stellar surface, with density following a power-law instead of an exponential decline in similarity speed. Similar shock breakouts occur in core-collapse SNe (Dessart et al., 2010; Ro & Matzner, 2013, 2017), with the resulting ejecta also following the simple similarity form for the velocity, $v=r/t$. Analysis of the mass and density distribution of such SN ejecta by Matzner & McKee (1999) shows a rather complex dependence on the structure of the pre-SN star; but to the extent that the ejecta density follows a similarity relation, Matzner & McKee (1999) fit it with a broken power law in $v=r/t$.

The appearance of a similar power-law tail in velocity for the core-heating cases here provides a potential bridge toward such scalings for SNe. One may conjecture that further increases in the energy addition factor to large $f\gg 1$, as occur in SNe, will lead to greater prominence of such a scale-free power-law. As gravitational binding becomes increasingly overwhelmed, there should be a concomitant diminution in the exponential similarity form, since this enforces the role gravity through the characteristic speed $v_{\rm o}$, which sets the mean speed to near the stellar escape speed. Detailed exploration of such transition from eruptive to explosive mass ejection is left for future work.

6 Parameter study

6.1 Dependence of ejecta mass and energy on level and location of energy addition

Let us next examine how the total mass and kinetic energy of the ejecta depend more broadly on the energy addition factor $f$, as well as on the location of this added energy. For this we have computed a grid of models with energy factors ranging from $f=0.1$ up to $f=2.0$ and with three distinct locations for the energy addition, namely the core ($M_{r}<76{\rm M}_{\odot}$), envelope ($M_{r}>76{\rm M}_{\odot}$), and near-surface ($M_{r}>89{\rm M}_{\odot}$). Each model is computed to a final time (typically $\sim 1$ Ms) that is sufficiently long to show convergence for the mass coordinate with zero Bernoulli energy, $B_{\rm n}(M_{e},t)\rightarrow 0$, giving then for the total ejected mass $\Delta M=M_{\ast}-M_{e}$.

As in the fiducial model (cf. eqn. 17), the total kinetic energy of the ejecta can be estimated from the drop in the energy of the star from the initial to final state,

[TABLE]

where the latter can be readily determined from the simulation. From eqn. (16), the total mass and kinetic energy of the ejecta can then be used to compute its energy-averaged speed ${\bar{v}}$.

The top panel of figure 9 shows that, for all three locations of the energy deposition, the ejected mass increases with $f$; but at the largest $f\rightarrow 2.0$ they each approach a distinct asymptote, set by the total mass above the onset of heating, as indicated by the horizontal dotted lines. In effect, even a large energy addition only ejects mass that lies above the deepest energy deposition. For the core heating model, this means all the mass can be ejected, but for the envelope and surface heating cases, this is much more limited, with asymptotic ejecta values of respectively $24{\rm M}_{\odot}$ and $11{\rm M}_{\odot}$, representing the entire mass of the heating region.

For most cases, the unbound mass is just a small fraction of the total mass, e.g. $<$10% for the f$<0.5$ models, which is the relevant range for direct energy input by stellar mergers. This is consistent with the $<$8% mass loss found in simulations of head-on collisions of high-mass stars (Glebbeek et al., 2013), though the details of the energy budget in those simulations are complex and thus still unclear. Simulations of the dynamical phase of common envelope evolution (e.g., Taam & Sandquist, 2000; Passy et al., 2012; Ricker & Taam, 2012; Nandez & Ivanova, 2016; Ohlmann et al., 2016; Iaconi et al., 2018) also indicate inefficient mass ejection, but until very recently (see, e.g., Ricker et al., 2018), these have generally concentrated on low-mass stars, which have significantly different stellar structures and energy budgets. In multi-D simulations the energy from mergers can take different forms like rotation or complex circulations not accounted for in the present idealized 1D simulations. But the breadth of this 1D parameter study can, with appropriate caution, provide a helpful general guide for interpreting the mass and energetics of ejecta in more realistic multi-D models.s

The middle panel of figure 9 shows the associated mean ejecta speeds ${\bar{v}}$, defined in eqn. (16). For small energy addition factor $f$, the mean speeds for all cases are comparable to the surface escape speed of the unperturbed star, ${\bar{v}}\approx v_{\rm esc}$, as shown by the horizontal dashed line. This is reminiscent of the tendency of a steady stellar wind to reach a terminal speed that scales with the surface escape speed, $v_{\infty}\sim v_{\rm esc}$ (Puls et al., 2008; Owocki, 2010), reflecting the fact that the driving needed to lift material from the star tends also to do a comparable amount of work against inertia to propel the wind’s kinetic energy. In the present eruptive context, it means that as long as the mass of the heated region does not limit the ejected mass, there is a rough equipartition between gravitational and kinetic energy in the ejecta. However, for envelope and near-surface energy deposition models, the saturation of the ejecta mass to that of the deposition region means that for higher energy cases $f>1$, the increased energy deposition leads to a higher specific kinetic energy, and thus higher ejecta speeds, ${\bar{v}}>v_{\rm esc}$. An analogous increase in stellar wind speed occurs when the mass loss from the surface is “choked”, e.g. by using a reduced density at the lower boundary (Owocki, 2010).

The bottom panel of figure 9 shows that, for small $f$, only a small fraction of sthe added energy is imparted to the ejecta as kinetic energy, but for larger $f\rightarrow 2.0$, the ejecta kinetic energy increases in direct proportion to the added energy, with an offset that increases with the depth of the energy addition. The difference $-f{E_{\ast}}-\Delta K$ represents energy left in the stellar envelope, leading to final net binding energy $E_{\rm fin}=(1-f){E_{\ast}}-\Delta K$ (cf. eqn. (21)). For the case of core energy deposition, we find $E_{\rm fin}\rightarrow 0$ as $f\rightarrow 2.0$, showing that the entire star becomes disrupted.

For the cases with energy deposition in the envelope or near surface, there always remains a negative net binding energy for the retained mass below the heating region. To the extent that this is above the binding energy of the original, unperturbed star (${E_{\ast}}=-6.65\times 10^{50}$ erg), it indicates the remaining star is likely to be left in a thermally excited state; if the associated luminosity exceeds the Eddington limit, this could drive a strong, super-Eddington wind for up to a thermal relaxation timescale, resulting in a post-eruption mass loss that could cumulatively rival that from the initial ejecta. This is an important topic for future investigation.

In summary, for modest energy factors $f<1$ the ejecta mass increases with increasing $f$, with ejecta speed comparable to the stellar escape speed. For large factors $f>1$, the ejection mass for the envelope and near-surface heating model is limited to the mass in those regions, with the speed and kinetic energy then increasing with increasing energy factor. By contrast, adding such large energy in the stellar core can completely disrupt the star, ejecting up to the entire stellar mass with average speed near the surface escape speed.

Appendix B develops a simple fitting function for the $f$-dependence of the mass ejecta (scaled by the mass above the bottom of the heating region, i.e. $\Delta M/M_{\rm bot}$), while figure 12 compares the visual best fits with simulation results for each of the three locations for the energy addition.

6.2 Applicability of similarity forms and relation to SNe

As noted at the end of section 5.2, the extension of eruption simulations to higher energy factors above unity, and the inclusion of models with the energy addition in the stellar core, provides a potential bridge to the properties of core-collapse supernovae (SNe), for which $f\gg 1$, with the energy concentrated near the center of the enriched core. A full discussion is beyond the scope of the current work aimed at modeling eruptive mass loss, for which the assumed pre-eruptive interior and core structure is much different than for SNe; but in this context, let us briefly consider the nature and applicability of similarity relations in SNe vs. the eruption models here.

While SN ejecta are commonly modeled in terms of a similarity form for the velocity $v=r/t$, the analysis by Matzner & McKee (1999) of their mass and density distribution shows a rather complex dependence on the structure of the pre-SN star. But to the extent that ejecta density follows a similarity relation, Matzner & McKee (1999) fit it with a broken power law in $v=r/t$. This is contrasts with the exponential density decline found for our eruptive models, $\rho\sim e^{-r/v_{\rm o}t}$, with the characteristic velocity $v_{\rm o}$ playing a key role in the resulting scalings of the ejecta mass and kinetic energy, as given by eqns. (12) and (15).

Analysis of our extended parameter study indicates that for energy addition factors $f\lesssim 1$, both the model series with surface and envelope energy addition also show this exponential decline in density, with the characteristic velocity following roughly the $v_{\rm o}\approx{\bar{v}}/\sqrt{12}$ scaling given in eqn. (16). The mass distribution of velocity likewise follows the similarity form given by eqn. (14).

For such cases with higher energy addition, $f>1$, the time evolution can be more complex, and not always well fit with this simple similarity form for density. The strong energy deposition in the envelope or near the surface excites oscillations in the underlying core, leading to repeated outward shock waves (similar to the reflected shock in the $f=0.5$ envelope model; see figure 2), and so an extended time interval before the dynamical evolution settles into a simple similarity form.

For the models with core energy addition there are no such reflected shocks, so we do generally find the evolution approaches a similarity form for the full range of energy factors $0.1<f<2.0$. Specifically, at moderate $v=r/t\lesssim 3000$ km s*-1*, the density agains shows an exponential decline $e^{-r/v_{\rm o}t}$, but then shifts to more of a power-law for large speeds $v\gtrsim 3000$ km s*-1*. We interpret this as resulting from the greater strengthening of the shock over the envelope propagation from its origins in the heated core, with then the surface layers expelled at high speed in a way similar to shock breakout in SNe (Dessart et al., 2010; Ro & Matzner, 2013, 2017) But for the greater mass in the subsurface layers, the non-negligible retarding effect of the gravity leads to an exponential density cutoff, with the characteristic velocity $v_{\rm o}$ reflecting the competition between gravity and energy addition that sets the bifurcation between bound and ejected mass. In contrast, the models with energy added to the envelope or near-surface don’t have this extended shock growth, and so have little or no material ejected with the high speeds associated with this power-law tail in similarity variable $v=r/t$.

There is also a potential analogy between such eruptions and possible failed explosions from core-collapse (Dessart et al., 2010; Lovegrove & Woosley, 2013; Coughlin et al., 2018). In the latter, the energy from collapse escapes by neutrino emission, with insufficient coupling to the stellar envelope to induce an explosion; but the mass lost in neutrinos and the associated reduced gravitational binding lead to a weak, low-speed ejection of some of the envelope mass (Lovegrove & Woosley, 2013) in a way that is in some respects similar to the eruptions from the mass-proportional addition of energy in the models here.

We thus see that the core heating models here do indeed suggest some potential links between such eruptions and both SNe explosions and failed explosions, but we leave a more complete examination of this to future work, which should include examination of the potential role of fundamental differences in the pre-explosion stellar structure.

7 Concluding summary and future outlook

Motivated by the eruptive mass loss inferred from LBV stars, particularly the 1840’s great eruption of $\eta$ Carinae, this paper uses time-dependent hydrodynamical simulations to examine the mass ejection that results from impulsive energy addition to the interior of a very massive star. This contrasts with previous models of steady-state energy addition that results in a continuum-driven super-Eddington stellar wind (Quataert et al., 2016; Owocki et al., 2017). Moreover, unlike supernova explosion models – for which the extreme energy addition in the core leads to complete, high-speed ejection of the overlying stellar material –, the energy addition here is at an order-unity factor of the stellar binding energy, and so results in ejection of only some outer fraction of the stellar mass, with the rest falling back onto an expanded, though still bound, surviving star.

As in SNe, the velocity of the expanding ejecta closely follows a similarity relation in radius and time, $v\approx r/t$. But a further, key result is that the density of this expansion closely fits a declining exponential $\rho\sim e^{-r/v_{\rm o}t}$ (see figure 5 and eqn. (10)), with now a characteristic velocity scale $v_{\rm o}$ that enforces the effect of gravitational binding in limiting the ejected mass to the outermost material. This contrasts with the scale-free, power-law similarity scalings inferred for SNe (Matzner & McKee, 1999), and to our knowledge is a distinctly novel property for such eruptive mass loss with a bifurcation between ejected and bound mass (see, however Dwarkadas & Chevalier (1998)).

Integration of this exponential similarity form for density yields analytic scalings for the ejected mass and kinetic energy (eqns. (12) and (15)), with an associated analytic form for the mass distribution of the velocity (see figure 6 and eqn. (14)). An associated scaling for the temperature (eqn. (26)) should prove useful for future analysis of the observational signatures of such eruption models.

Our broader parameter study shows how the ejecta mass and kinetic energy depend on the level and location of the energy addition222It likely also will depend on the stellar structure, but since we here have fixed the stellar model, quantifying this will require future investigation; see below.. For energy addition less than the stellar binding energy ($f<1$), the ejected mass increases with increasing $f$, with an average flow speed generally near the surface escape speed, ${\bar{v}}\approx v_{\rm esc}$. For energy addition exceeding the binding energy by up to a factor $f=2$, results depend on the location of the energy deposition. For addition to the envelope or near-surface, the ejected mass is limited to the mass of the heated region, with the excess energy thus increasing the ejecta’s specific energy and thus flow speed. For addition in the stellar core, such a large energy factor can lead to ejection of nearly all the stellar mass, but with slower speed. This provides a potential bridge between the present eruption models and full SN explosions that warrants further investigation.

In the present context of exploring the nature of eLBV stars, the dynamical simulations here have assumed a very massive star at an intermediate main-sequence evolutionary phase. Future work will be needed to explore how the specifics of mass ejection depend on stellar mass and evolutionary state. For example, in their exploration of a binary merger scenario for SN1987A, Morris & Podsiadlowski (2007, 2009) analyzed how this could lead to anisotropic ejection of the extended envelope of a red supergiant. The low binding energy of this envelope facilitates its ejection, but the ejected mass is much lower than in the eLBV scenarios explored here. Overall, the form, speed and mass of the ejecta should depend on both the details of the energy and angular momentum addition, and the mass and evolutionary state of the merging stars. It will be of interest to explore how some of the key findings of the present study, e.g. the exponential similarity of the ejecta density, may extend to such a broader range of mass eruptions.

In general, the hydrodynamical simulations and similarity relations presented here provide important new insights into the properties of the eruptive mass loss that can arise from impulsive energy addition to various regions of a massive-star interior. The scales of the ejected mass, over a range from $<0.1{\rm M}_{\odot}$ to $>10{\rm M}_{\odot}$, and of the ejecta speeds, from a few hundred to several thousand km s*-1*, are comparable to what’s inferred from eLBV’s (Smith et al., 2011). Moreover, a small amount of material is ejected at a speed much higher than the average, providing a possible rationale for the high-speed ejecta inferred by Smith et al. (2018) for $\eta$ Carinae (see also Smith, 2006, 2008, 2013).

But clearly the 1D, adiabatic models here – with an assumed impulsive energy addition of unspecified physical origin – represent only one further step toward understanding eruptive mass loss from massive stars. Future improvements should include a treatment of radiative transport, particularly for the near-surface layers that will set the observational signatures of the eruptions. Any more detailed comparison with specific eLBV’s, most notably $\eta$ Carinae, will also require significant extensions, to include a specific mechanism for the energy addition, and to identify its level and location, and the degree to which it can be characterized as impulsive, steady, or, most likely, something intermediate between these idealizations. In particular, in the case of a binary merger model for $\eta$ Carinae, there is a need for 2-D and 3-D dynamical models that account for the role of the binary orbital angular momentum in spinning up the merged star and centrifugally ejecting mass, while also accounting for the energy deposition from merger of the separate stellar cores.

Acknowledgments

We acknowledge Nathan Smith for helpful discussions and many useful comments and suggestions on an early draft. We thank Chris Matzner for pointing out to us previous findings of exponential similarity scalings for density in SNIa explosions. We also thank Lorne Nelson for sharing computing facilities, and the referee, Alex de Koter, for constructive criticisms that helped improve the final manuscript. The computations here were partially carried out on facilities managed by Calcul Québec and Compute Canada. RH was supported by the JSPS Overseas Research Fellowship No.29-514 and a grant from the Hayakawa Satio Fund awarded by the Astronomical Society of Japan. SPO acknowledges support of a visiting scholar grant from the Royal Society. This work has also been supported by a Humboldt Research Award to Ph.P. at the University of Bonn, and by the Oxford Hintze Centre for Astrophysical Surveys, which is funded through generous support from the Hintze Family Charitable Foundation. Many of the ideas for this study grew out of discussions during a massive-star research program at the Kavli Institute for Theoretical Physics, and as such was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915.

Appendix A Similarity relations for internal energy and temperature

Let us also examine the time and space variation of the temperature for ejected material within this idealized model of adiabatic expansion. To do this, we first consider the variation of total internal energy $E=3P_{\rm rad}+(3/2)P_{\rm gas}$, as given by eqn. (5). Applying a similarity form for the velocity $v=r/t$, this can be written as

[TABLE]

Defining $\beta\equiv P_{\rm gas}/{P}$, we can recast this in the form,

[TABLE]

This implies that, if $\beta$ were constant, the internal energy of mass parcels with fixed $M_{r}$ would decline as a simple power law in time, $E\sim 1/t^{\delta}$. Note in particular that the power index ranges from $\delta=4$ for the radiation-dominated ($\beta=0$) case, to $\delta=5$ in the gas-dominated ($\beta=1$) limit, with $\delta=4.33$ for the special case $\beta=1/2$ with equal gas and radiation pressure.

For the fiducial model with energy-addition factor $f=0.5$ applied in the stellar envelope, figure 10 shows a contour plot of $\beta$ as a function of mass $M_{r}$ and time $t$. The heavy red line shows the equipartition value $\beta=1/2$, while the other contours show that for most of the ejecta, $\beta$ remains within about $\pm 0.2$ of this equipartion value.

Figure 11 plots time snapshots of the internal energy scaled by time as $Et^{4.33}$, i.e. with the exponent $\delta=4.33$ appropriate for $\beta\approx 0.5$. The dashed line represents the best-fit linear trend, corresponding again in this semi-log plot to an exponential decline $e^{-r/v_{\rm o}t}$, with $v_{\rm o}\approx$ 325 km s*-1*. The total energy thus follows a similarity relation of the form,

[TABLE]

where the proportionality constant has the CGS value $E_{\rm o}t_{\rm o}^{4.33}\approx 6.\times 10^{31}$ erg cm*-3* s4.33.

If we also apply the approximation $\beta\approx 1/2$ in the expression for total energy, we find

[TABLE]

The overall similarity fit for the temperature thus takes the approximate form,

[TABLE]

where the scaling constant

[TABLE]

These similarity relations for temperature, density, and velocity can be used, together with a prescription for the opacity, to model observational properties of the eruption, such as the light curve, or even the evolving spectrum. But reliable results will also require inclusion of a radiation transport and the associated radiative losses near the expanding photosphere. This will be a topic for future work.

Appendix B Fits to ejecta mass vs. energy fraction

The parameter study in section 6 shows that the mass ejected at high energy-addition factor $f>1$ is limited to the mass above the bottom of the heating region, $M_{\rm bot}$. But in the opposite limit of small energy factor, $f\ll 1,$ the mass ejection scales roughly as a power in $f$. Let us thus define a general fitting function for the scaled mass ejection that bridges these two limits,

[TABLE]

where the parameter $f_{\rm b}$ represents the transition from the power-law form $\Delta M\sim f^{n}$ at low $f\ll f_{\rm b}$, to the saturation $\Delta M\rightarrow M_{\rm bot}$ at large $f>f_{\rm b}$. For the parameter values in Table 2, figure 12 compares this fitting function to the mass ejection results for each of the three locations of energy addition.