# Wave-Driven Shocks in Stellar Outbursts: Dynamics, Envelope Heating, and
Nascent Blastwaves

Christopher D. Matzner1 & Stephen Ro2<EMAIL_ADDRESS>1Department of
Astronomy & Astrophysics, University of Toronto, 50 St. George St., Toronto,
ON M5S 3H4, Canada  2Astronomy Department and Theoretical Astrophysics Center,
University of California, Berkeley, Berkeley, CA 94720

(Received September 25, 2020)

###### Abstract

We address the shocks from acoustic pulses and wave trains in general one-
dimensional flows, with an emphasis on the application to super-Eddington
outbursts in massive stars. Using approximate adiabatic invariants, we
generalize the classical equal-area technique in its integral and differential
forms. We predict shock evolution for the case of an initially sinusoidal but
finite wave train, with separate solutions for internal shocks and head or
tail shocks, and demonstrate detailed agreement with numerical simulations.
Our internal shock solution motivates improved expressions for the shock
heating rate. Our solution for head and tail shocks demonstrates that these
preserve dramatically more wave energy to large radii and have a greater
potential for the direct ejection of matter. This difference highlights the
importance of the waveform for shock dynamics. Our weak-shock analysis
predicts when shocks will become strong and provides a basis from which this
transition can be addressed. We use it to estimate the mass ejected by sudden
sound pulses and weak central explosions.

††software:  FLASH4.6 (Fryxell et al. 2000)

##  1\. Introduction

A fraction of massive stars undergo periods of intense, episodic mass loss as
they approach core collapse. Observations of eruptions in luminous blue
variables, supernova (SN) impostors, and pre-SN outbursts provide direct
evidence of these events (Foley et al., 2007; Pastorello et al., 2007;
Margutti et al., 2014; Smith et al., 2014), while dense and massive zones of
circumstellar media (CSM) provide indirect evidence. Dense CSM is revealed by
its interaction with supernova (SN) ejecta (Pastorello et al., 2008a; Smith et
al., 2011; Kiewe et al., 2012; Moriya et al., 2014; Pastorello et al., 2016;
Shivvers et al., 2016; Hosseinzadeh et al., 2017; Shivvers et al., 2017), or
by recombination after ‘flash ionization’ from SN shock breakout (Gal-Yam et
al., 2014; Groh, 2014; Khazov et al., 2016). SN types IIn and Ibn require
especially dense CSM from previous eruptions (Chugai et al., 2004; Pastorello
et al., 2008b), and SN spectral evolution from Type Ib to IIn (Milisavljevic
et al., 2015; Margutti et al., 2017) and the transitional Type IIn/Ibn SNe
(Pastorello et al., 2008c; Smith et al., 2012) show strong interactions with
dense, hydrogen-rich CSM.

Importantly, ejecta speeds in some outbursts are high enough to require
acceleration by shocks, as opposed to winds. The pre-explosion outburst of SN
2009ip emitted material at 13,000 km s-1(Pastorello et al., 2013), and there
is tentative evidence of 15,000 km s-1ejecta in the similar source SN 2015bh
(Ofek et al., 2016). The famous Galactic example $\eta$ Car has survived
multiple eruptions; its last Giant Eruption ejected traces material at SN
velocities (${\sim}20\,000$ km s-1; Smith et al. 2018) along with about
${\sim}15$ $\mathrm{M}_{\odot}$ of slower (${\sim}600$ km s-1) material over
two decades (Smith, 2008).

Such phenomena require that energy escapes the stellar interior by a conduit
capable of transporting super-Eddington luminosities and creating shock
motions along the way. Sound waves are a prime suspect for this conduit, for
which the energy sources may include unstable phases of oxygen-silicon shell
burning in low-mass pre-SN stars (Woosley & Heger, 2015) or pulsational pair
instability in high-mass pre-SN stars (Woosley et al., 2007). More commonly,
waves from vigorous convection in advanced burning regions (Quataert & Shiode,
2012; Shiode & Quataert, 2014; Smith & Arnett, 2014) can transport energy to
the outer envelope at super-Eddington rates. Wave dissipation is then thought
to trigger outflows with mass loss rates and velocities characteristic of Type
IIn SNe CSM (Quataert et al., 2016; Fuller, 2017; Fuller & Ro, 2018).

Thanks to the pulsational pair instability, wave-driven mass loss plays a role
in determining the lower edge of an upper black hole mass gap (Woosley, 2017;
Leung et al., 2019; Farmer et al., 2019), a phenomenon that is currently of
interest due to the high inferred initial masses in the black hole coalescence
GW190521 (Abbott et al., 2020).

We seek a way to accurately predict the formation and evolution of wave-driven
shocks and the distribution of shock-deposited heat. Critically, a complete
theory must also evaluate the conditions under which some of the wave energy
persists all the way to the stellar surface, driving an eruption. Closed-form
solutions are especially useful: while waves and shocks can be calculated
using hydrodynamical codes, the very wide range of physical scales makes
direct simulation very costly in practice. Thousands of cells are needed per
acoustical wavelength even with high-order codes, and the wavelength may be
very small compared to the stellar radius.

Fortunately, very accurate predictions are possible using an analytical
technique, which in fact is most accurate in the same (low-amplitude and high-
frequency) limit in which direct numerical simulations are the least feasible.
It is also a classic technique, as it builds on a proposal by Lighthill (1978)
for the generalization of a method due to Whitham (1974), which in turn
follows a suggestion of Landau (1945) on how to extend the seminal theory by
Riemann (1860) past the point of shock formation. We start here by tracing
these contributions, and reviewing our own work on the validity of Lighthill’s
proposal. We then extend the theory somewhat by re-expressing it in
differential form, before considering its implications for shock dissipation
within stars. Our results improve upon, and generalize, the shock heating
prescription developed by Fuller & Ro (2018). We shall see that the potential
for shock emergence depends critically on waveform of the driving source, as
well as its amplitude and period.

Notation: Our space and time coordinates are $r$ and $t$; $r$ will be the
stellar radius when we consider spherical symmetry. We use $\rho$, $c$, $p$,
$\gamma$, $v$, and $A$ for the density, sound speed, pressure, adiabatic
index, fluid velocity, and cross-sectional area. For simplicity we consider
only ideal gases, in which $\gamma$ is fixed in each fluid element (but may
vary across the structure); and we define the auxiliary quantity
$q=(\gamma-1)/2$. The prefix $\delta$ indicates perturbations from the
background state, which is indicated by subscript $0$. The jump of a quantity
across a shock is denoted by prefix $\Delta$.

##  2\. The equal-area method and its generalization

The simplest problem in nonlinear fluid dynamics involves a constant-entropy
fluid moving within a channel of constant cross-sectional area (or in planar
symmetry), in the absence of viscosity or any body forces. This problem was
solved exactly by Riemann (1860), who re-expressed the fluid equations in
terms of quantities $J_{\pm}=v\pm 2c/(\gamma-1)$ that are conserved along
sound fronts: outward or inward trajectories defined by $\dot{r}=v\pm c$, also
known as $C^{\pm}$ characteristics. Riemann’s solutions apply only while the
flow is continuous; however, they generically produce discontinuities after
some finite time. This is most obvious for ‘simple’ waves, in which one of the
two wave families is constant: for instance, if $J_{-}=2c_{0}/(\gamma-1)$ for
some constant $c_{0}$, then $v=2(c-c_{0})/(\gamma-1)=J_{+}/2$ and the wave
speed $v+c=(\gamma+1)v/2+c_{0}$ is a pure function of $J_{+}$ or of $v$. This
implies that the velocity waveform $v(r-c_{0}t)$ will shear at a strictly
uniform rate, as shown in Figure 1, with wavefronts moving horizontally (at
constant $v$). Within this shearing waveform, the inverse velocity gradient
$(\partial v/\partial x)^{-1}$ increases linearly with time, at the uniform
rate $(\gamma+1)/2$, along each wavefront. Therefore, compressive regions of
the waveform – those in which $\partial v/\partial r$ is initially negative –
will develop infinite slope in a finite time. This happens first on whichever
sound front or fronts carry the greatest compression, at a time
$t_{s}=[2/(\gamma-1)]\min(-\partial v/\partial x)_{t=0}^{-1}$. After
$t=t_{s}$, Riemann’s solution requires modification because it predicts a
shearing waveform that has ‘broken’ and become multi-valued, whereas in
reality the flow is single-valued and contains shock discontinuities.

The ‘equal-area’ method (Whitham, 1974, following Landau 1945; see also Landau
& Lifshitz 1959 § 102) provides the modified waveform, including the location
and strength of the shocks that form. One assumes that Riemann’s solution
remains valid away from the shock discontinuities; shocks simply remove a
portion of the waveform. Mass and momentum conservation require that the
removed sections of $v(r)$ must be of equal area, and this determines the
shock’s location as a function of time. Because the waveform shears in a
simple and uniform way, as do the omitted areas, one can identify the shock
transition at time $t$ directly from the initial waveform (at $t=0$), simply
by determining a line of appropriate slope that cuts away sections of equal
area. Figure 1 provides an example.

Once the strength of each shock is known, there are two ways to obtain the
rate at which it transfers wave energy ($E_{w}$) into heat ($E_{\rm diss}$)
within the stellar envelope. One is to employ the local formula

$\frac{dE_{\rm
diss}}{dr}=-\frac{dE_{w}}{dr}=\frac{\gamma+1}{12}A\rho_{0}\frac{(\Delta
v)^{3}}{c_{0}},$ (1)

where $\Delta v$ is the velocity jump across the shock. The other involves
inspecting the waveform to determine how much of its total energy has been
removed by the excision of sections around the shock discontinuities.

While Riemann’s solution is exact and nonlinear, this cannot be said for the
equal-area method. However, its approximations are highly accurate so long as
the wave amplitude is small ($v\ll c_{0}$) and the shocks are correspondingly
weak. A sufficiently weak shock adds negligible entropy to the fluid, as the
entropy perturbation is of order the shock strength cubed; nor does it create
a significant reflected wave (change in $J_{-}$), either of which would ruin
Riemann’s solution for the fluid downstream of the shock. In the small-
amplitude limit, moreover, the perturbations of the conserved mass density
$\rho$ and momentum density $\rho v$ are proportional to that of $v$, so the
equal-area constraint applies to any of these waveforms.

Figure 1.—  Whitham’s equal-area technique applied to a two-cycle outward-
traveling sinusoidal wavetrain of amplitude $v_{p}$, in the simple problem of
an isentropic fluid within a constant-area channel. The invariance of $J_{+}$
along $C^{+}$ characteristics leads the waveform to shear at a uniform rate,
so shocks form within the waveform at the locations of maximal $-\partial
v/\partial x$ at time $t=t_{s}$ (red dots). Shock transitions at time $t$
shear back to lines in the initial conditions (black dotted lines,
corresponding to shocks at $t=3t_{s}$) that cut away equal positive and
negative areas of the initial waveform (dashed blue regions). Internal shocks
evolve more rapidly than head and tail shocks, with implications for energy
deposition and shock ejection.

###  2.1. Lighthill’s generalization

In general, unlike the simple problem, the cross-sectional area, fluid
density, and sound speed may all vary with $r$. A further modification is thus
required to adapt the equal-area method to more complicated environments like
stellar envelopes, with their spherical geometries and strongly stratified
density profiles. For this, Lighthill (Waves in Fluids, 1978, chapter 2)
proposed to generalize Whitham’s method by replacing $J_{\pm}$ with an
approximate invariant (a linearized adiabatic invariant, to be precise) of
acoustical motions in the more complicated problem. Lighthill chose the
approximate invariant $\delta p/\sqrt{Z}$, where $\delta p$ is the pressure
perturbation and $Z(r)=\rho_{0}c_{0}/A$ is the acoustic impedance. However we
prefer to use the ‘wavefront amplitude’,

$\mathcal{R}_{\pm}=\sqrt{\rho_{0}c_{0}A}~{}~{}\delta J_{\pm}$ (2)

(Matzner & Ro, 2020, hereafter MR20). Here, $\delta J_{\pm}$ is the acoustical
perturbation of $J_{\pm}$ from the background state. Our quantity
$\mathcal{R}_{\pm}$ is also an approximate invariant. Its advantage is that it
provides a more direct connection to $J_{\pm}$ and also to the acoustical wave
luminosity, which for waves of linear amplitude is

$L_{w}=\frac{1}{4}(\mathcal{R}_{+}^{2}-\mathcal{R}_{-}^{2}).$ (3)

Conservation of $\mathcal{R}_{+}$ along characteristics in a purely outward-
traveling wave ($\mathcal{R}_{-}=0$) is therefore equivalent to the
conservation of $L_{w}$ in a wave that neither reflects nor dissipates.

In the generalized method, the wavefront shears at a changing rate that
depends on the local conditions. This can be accommodated by defining a new
evolution coordinate $y(r)$ to replace $t$, as we describe in more detail in §
3. Before we do, we wish to review a couple results from MR20 that help to
clarify the validity of the generalized method.

First, it is important to note that the quantity $\mathcal{R}_{\pm}$ is not
strictly conserved along characteristics, but rather varies in linear theory
according to

$\frac{d}{dr^{\pm}}\mathcal{R}_{\pm}=\mathcal{R}_{\mp}\frac{\partial}{\partial
r}\ln\sqrt{Z}$ (4)

where $d/dr^{\pm}$ is the derivative along $C^{\pm}$ with respect to $r$. The
right-hand side causes reflection whenever the impedance $Z$ varies
significantly on scales of order the wavelength, but is completely negligible
when it does not (the high-frequency limit). Therefore, Lighthill’s
generalized equal-area method requires that the wave not only be linear in
amplitude, but high in frequency relative to its surroundings. This is not a
strong constraint in the process of shock formation, because the dynamical
timescales are much shorter in the core than in the envelope of an evolved
star, and because a wave must be low-amplitude to travel many wavelengths
before it forms a shock. However, if the wave were to meet a rapid change in
$Z$ (at a composition boundary, for instance), then one would need to evaluate
what fraction is transmitted and apply our technique to the transmitted wave.
Reflection may also affect shocks that accelerate near the stellar surface, as
we discuss in § 6.3.

Second, equation (4) shows that $\mathcal{R}_{\pm}$ is an adiabatic invariant
of the linearized problem, but shock formation is a nonlinear process. In MR20
we show that nonlinear adiabatic invariants, when these can be found, almost
always tend toward $\mathcal{R}_{\pm}$ in the linear limit. More generally, we
find that the nonlinear error one incurs by assuming that $\mathcal{R}_{\pm}$
is constant is second-order in the wave amplitude, and also proportional to
the logarithmic gradient of a combination of environmental quantities. Note
also that the ‘conserved area’ $\int\mathcal{R}_{+}(t_{i})dt_{i}$ has units of
the square root of the action, which should be conserved in adiabatic
evolution. These results give us confidence that Lighthill’s proposed
generalization is sound.

##  3\. The method in differential form

Here we shall present an analysis that is physically equivalent to the equal-
area method, but expressed in differential rather than integral form. As
before, we focus on outward-traveling waves that encounter no inward sound
energy, and that also are sufficiently high-frequency to suffer negligible
reflection; accordingly, $\mathcal{R}_{-}=\delta J_{-}=0$. This implies
$\delta J_{+}=2v$ and $L_{w}=A\rho_{0}c_{0}v^{2}$, so

$\frac{v^{2}}{c_{0}^{2}}=\frac{L_{w}}{L_{{\rm max}}(r)}$

where $L_{{\rm max}}(r)=A\rho_{0}c_{0}^{3}$.

We label the $i$-th characteristic by the time $t_{i}$ at which it crosses a
fiducial radius $r_{1}$, and by its conserved amplitude $\mathcal{R}(t_{i})$
and luminosity, $L_{w}(t_{i})=\mathcal{R}(t_{i})^{2}/4$. It is convenient to
define the normalized waveform $f(t_{i})=\mathcal{R}(t_{i})/\mathcal{R}_{p}$
where $\mathcal{R}_{p}$ is the amplitude at the wave peak (and
correspondingly, $L_{p}=\mathcal{R}_{p}^{2}/4$). Travelling at speed
$(v+c)_{i}=c_{0}+(q+1)v_{i}$, characteristic $i$ arrives at radius $r$ at the
time

$\displaystyle t$ $\displaystyle=$ $\displaystyle
t_{i}+\int_{r_{1}}^{r}\frac{dr^{\prime}}{(v+c)_{i}}$ (5) $\displaystyle=$
$\displaystyle
t_{i}+\int_{r_{1}}^{r}\frac{dr^{\prime}}{c_{0}(r^{\prime})}-f(t_{i})y(r),$

Above, we take the Taylor expansion of the right hand side in terms of
$v_{i}/c_{0}\ll 1$ to generate the second equation, and we define two new
quantities:

$y(r)=L_{p}^{1/2}Y(r)$ (6)

and

$Y(r)=\int_{r_{1}}^{r}\frac{dr^{\prime}}{c_{0}(r^{\prime})}\frac{\gamma(r^{\prime})+1}{2L_{\rm
max}^{1/2}(r^{\prime})}.$ (7)

The evolution coordinate $y(r)$ takes the role played by $t$ in the simple
problem, in that it determines how much the waveform has sheared.

Since $t$ is a function of $r$ and $t_{i}$, the partial derivative with
respect to $t_{i}$ commutes with the integration over $r$. The density of
characteristics is indicated by $\partial t_{i}/\partial t$, which is the
inverse of

${\partial t}{/\partial t_{i}}=1-f^{\prime}(t_{i})y(r).$ (8)

The wave variation evolves according to

$\frac{\partial f}{\partial t}=\frac{df}{dt_{i}}\frac{\partial t_{i}}{\partial
t}=\frac{f^{\prime}(t_{i})}{1-f^{\prime}(t_{i})y(r)}.$ (9)

Characteristics converge (disperse) with time in compressive (rarefactive)
regions of a wave, since $y(r)$ is positive and increasing. Note that
$\partial f/\partial t$ and $\partial v/\partial r$ are opposite in sign:
compressive regions have $\partial f/\partial t>0$.

###  3.1. Shock formation

A shock forms where $\partial f/\partial t$ becomes infinite; this occurs
first on those characteristics, labeled $t_{i,s}$, that are local maxima of
the compression rate and therefore local maxima of $f^{\prime}(t_{i})$. (This
is a necessary but not always sufficient criterion: it is possible for such
characteristics to collide with another shock before forming their own.) The
shock formation radius $R_{s}$ is found by the relation

$y_{s}\equiv y(R_{\rm s})=\frac{1}{f^{\prime}(t_{i,s})}.$ (10)

Nearby characteristics will collide with the shock: these correspond to the
sections of the waveform excised in the equal-area method.

The shock formation criterion, equation (10), was derived by Lighthill (1978,
his equation 254), and re-derived by Lin & Szeri (2001), Tyagi & Sujith
(2003), Tyagi & Sujith (2005, who noted agreement with Lighthill’s criterion)
and us (Ro & Matzner, 2017). As we discuss in MR20, it assumes the wave’s
slope is not altered by reflection at some location where $Z$ changes rapidly;
this is more of a restriction for internal shocks than head shocks.

###  3.2. Weak Shock Propagation

To address the trajectory $R(t)$ of a shock after it forms, we follow the
characteristic-tracking procedure of Friedrichs (1948) and Whitham (1974). The
velocity $V=dR/dt$ of a weak shock is just the average of the characteristic
velocities immediately to each side of it:

$V=\frac{(v+c)_{u}+(v+c)_{d}}{2},$ (11)

to first order in $v/c_{0}$ (Courant & Friedrichs, 1948), where $u$ and $d$
refer to the states immediately upstream and downstream of the shock,
respectively. This condition enforces the equal-area rule we discussed in §
2.1, but in differential form.

In order to evaluate $V$, we require the rate at which characteristics arrive
at the shock front from either side:

$\left.\frac{dt_{i}}{dt}\right|_{s}=\partial_{t}t_{i}+V\partial_{r}t_{i}=\frac{v+c-V}{v+c}\partial_{t}t_{i};$

the second equality uses $\partial_{t}t_{i}=-(v+c)\partial_{r}t_{i}$.
Substituting equation (8),

$\left(1-f^{\prime}(t_{i})y\right)\left.\frac{dt_{i}}{dt}\right|_{s}\simeq\frac{v+c-V}{v+c},$

which can be evaluated on either side of the shock. Combining this with
equation (11) and using $(q+1)\sqrt{L_{p}/L_{{\rm
max}}}\left.d/dt\right|_{s}=\left.d/dy\right|_{s}$ gives coupled differential
equations for the initial times of characteristics as they arrive at the
shock:

$\displaystyle 2\left(1-f^{\prime}(t_{id})y\right)\frac{dt_{id}}{dy}$
$\displaystyle\simeq$ $\displaystyle f(t_{id})-f(t_{iu}),$ (12a)
$\displaystyle 2\left(1-f^{\prime}(t_{iu})y\right)\frac{dt_{iu}}{dy}$
$\displaystyle\simeq$ $\displaystyle f(t_{iu})-f(t_{id}).$ (12b)

The shock trajectory, described by $t_{iu}(y)$ and $t_{id}(y)$, is found by
integrating these equations from the point of shock formation,
$t_{iu}=t_{id}=t_{i,s}$. When it is desirable to work with dimensionless
solutions scaled to the formation of a given shock, we note that equations
(12a) and (12b) remain valid and unchanged with the replacements
$t_{i}\rightarrow\tilde{t}_{i}$ and $y\rightarrow\tilde{y}$, where

$\tilde{t}_{i}\equiv\frac{t_{i}-t_{i,s}}{y_{s}},~{}~{}~{}~{}~{}~{}\tilde{y}(r)\equiv\frac{y(r)}{y_{s}}=\frac{Y(r)}{Y_{s}},$

and $Y_{s}$ means $Y(R_{s})$.

To specify the deposition of wave energy we define a normalized energy
variable $\tilde{E}=E/(L_{p}y_{s})$. Equation (1) becomes

$\frac{d\tilde{E}_{\rm diss}}{dr}=\frac{(\Delta
f)^{3}}{6}\frac{d\tilde{y}}{dr}=\frac{d\tilde{E}_{\rm
diss}}{d\tilde{y}}\frac{d\tilde{y}}{dr}$ (13)

We see from this that shock dissipation can be broken into two parts. The
first, $d\tilde{E}_{\rm diss}/d\tilde{y}=(\Delta f)^{3}/6$, depends only on
the form of the initial pulse and the evolution coordinate $\tilde{y}$. The
second, $d\tilde{y}/dr$, depends only on the wave amplitude and the structure
of the background envelope.

##  4\. Particular Solutions

Equations (12b) and (12a) are both of the form

$\left(f(t_{id,iu})-f(t_{iu,id})\right)y^{\prime}=2\left(1-f^{\prime}(t_{id,iu})y\right),$

which is solvable if there is a definite relationship between the values of
$f$ just upstream and downstream of the shock. Suppose
$f(t_{id,iu})-f(t_{iu,id})=mf(t_{id,iu})$ for some non-zero constant $m$; then

$b\left(1-f^{\prime}(t_{i})y(t_{i})\right)=f(t_{i})y^{\prime}(t_{i}),$ (14)

where $b=2/m$ and $t_{i}$ is shorthand for either $t_{id}$ or $t_{iu}$. The
exact solution is

$y(t_{i})=y_{s}\frac{f(t_{i,s})^{b}}{f(t_{i})^{b}}+\frac{b}{f(t_{i})^{b}}\int_{t_{i,s}}^{t_{i}}f(t_{i}^{\prime})^{b-1}dt_{i}^{\prime}$
(15)

The dimensionless form of this solution is exactly the same, with
$\\{t_{i,s},t_{i},y\\}\rightarrow\\{0,\tilde{t}_{i},\tilde{y}\\}$.

For a shock that forms at a wave node ($f(t_{i,s})=0$), the first term in
equation (15) is zero and can be omitted.

This solution addresses some common cases. Internal shocks are those for which
there is wave power both upstream and downstream of the shock. When $f(t_{i})$
is locally symmetric around $t_{i,s}$, we know that
$t_{iu}-t_{i,s}=t_{i,s}-t_{id}$ and $f(t_{iu})=-f(t_{id})$; these shocks are
described by $m=2$ and $b=1$. Internal shocks travel at the undisturbed sound
speed ($V=c_{0}$), to first order in the wave amplitude. Head and tail shocks
are those in which a wave packet suddenly begins or ends in a compressive
phase. For these, either $f_{id}$ or $f_{iu}$ is zero, so they are described
by $m=1$ and $b=2$. Leading shocks travel somewhat faster than $c_{0}$;
whereas trailing shocks travel slightly slower than $c_{0}$ (and yet are
supersonic relative to the upstream fluid).

###  4.1. Shocks in Sinusoidal Waves

Consider a sinusoidal wave of $N$ cycles, $f=\sin\omega t_{i}$ for $0<\omega
t_{i}<2N\pi$. We assume that it begins and ends in a state of compression, so
that head and tail shocks form in addition to its internal shocks. All these
shocks share the common value $y_{s}=1/\omega$, so $\tilde{y}=\omega y$ and
one can define $\tilde{t}_{i}=\omega(t_{i}-t_{i,s})$ for any shock labeled by
$t_{i,s}$. The example in Figure 1 consists of two cycles.

####  4.1.1 Internal shocks

Internal shocks evolve according to equation (15) with $b=1$ and
$f(t_{i,s})=0$, so
$\tilde{y}(\tilde{t}_{i})=\tilde{t}_{i}/\sin(\tilde{t}_{i})$ and

$\tilde{t}_{i}={\rm sinc}^{-1}\frac{1}{\tilde{y}(r)}.$ (16)

This coincides with what one would obtain by construction for the internal
shock in Figure 1.

Symmetry around the internal shocks implies $u_{u}=-u_{d}$ and so $\Delta
f=2f_{d}$. Evaluating equation (13),

$\frac{d\tilde{E}_{\rm diss}}{d\tilde{y}}=\frac{4}{3}\left[\sin\left({\rm
sinc}^{-1}\frac{1}{\tilde{y}}\right)\right]^{3}$ (17)

and the cumulative dissipated energy is

$\tilde{E}_{\rm
diss}=\tilde{t}_{i}-\sin\tilde{t}_{i}\cos\tilde{t}_{i}-\frac{2}{3}\tilde{t}_{i}(\sin\tilde{t}_{i})^{2},$
(18)

which can be converted to $\tilde{E}_{\rm diss}(\tilde{y})$ using (16). We
obtained this by integrating the wave energy in the equal-area portions of the
initial waveform that will be removed by the shock (see Figure 1); we then
checked for numerical consistency with equation (17).

The total energy deposited by an internal shock is $\tilde{E}_{\rm
diss}(\infty)=\int_{1}^{\infty}(d\tilde{E}_{\rm
diss}/d\tilde{y})d\tilde{y}=\pi$, or in dimensional terms, $E_{\rm
diss}(\infty)=\pi L_{p}/\omega$. This matches the initial wave energy carried
by characteristics that will intersect the shock. We will also be interested
in what fraction of the wave energy survives at late times. At large values of
$\tilde{y}$, $d\tilde{E}_{\rm
diss}/d\tilde{y}\rightarrow(4\pi^{3}/3)\tilde{y}^{-3}$, so the residual wave
energy that will ultimately interact with the shock is

$\tilde{E}_{w}(\tilde{y})\rightarrow\frac{2}{3}\pi^{3}\tilde{y}^{-2}.$ (19)

####  4.1.2 Head and tail shocks

Head and tail shocks evolve in opposite directions but otherwise in exactly
the same way (in linear theory), so we focus on the case of a head shock.
Using $b=2$ and $f(t_{i,s})=0$, equation (15) yields
$\tilde{y}=2/(1+\cos\tilde{t}_{i})$, so

$\tilde{t}_{i}=\cos^{-1}\left(\frac{2}{\tilde{y}}-1\right).$ (20)

In this case $\Delta f=f(t_{i})$; using $\sin[\cos^{-1}(x)]=\sqrt{1-x^{2}}$,

$\Delta f=2\frac{\sqrt{\tilde{y}(r)-1}}{\tilde{y}(r)};$

and by equation (13)),

$\frac{d\tilde{E}_{\rm
diss}}{d\tilde{y}}=\frac{4}{3}\frac{(\tilde{y}-1)^{3/2}}{\tilde{y}^{3}}.$ (21)

The total energy deposited by a head or tail shock shock, $\tilde{E}_{\rm
diss}(\infty)=\pi/2$, is only half of what is deposited at an internal shock,
because external shocks interact with characteristics from only half a wave
period. However, external shocks expend this energy much more slowly. For
large $\tilde{y}$, external shocks obey $d\tilde{E}_{\rm
diss}/d\tilde{y}\rightarrow(4/3)\tilde{y}^{-3/2}$, so the residual wave energy
is

$\tilde{E}_{w}(\tilde{y})\rightarrow\frac{8}{3}\tilde{y}^{-1/2},$ (22)

in sharp contrast to internal shocks (eq. [19]). For this reason, internal
shocks deliver more energy to radii of order the shock formation radius, while
leading and trailing shocks have a much greater potential to create strong
shocks in stellar atmospheres. Figure 2 compares the normalized residual
energy $\tilde{E}_{w}(\tilde{y})$ between internal and head/tail shocks.

Although we concentrate on purely sinusoidal wave trains, it is important to
realize that the asymptotic properties of shock decay are not sensitive to the
details of the waveform. For head and tail shocks, for instance, the limit
$f(t_{id})\rightarrow 2\tilde{y}^{-1/2}$ is a special case of the result

$f(y)^{2}=\frac{2}{y}\int f(t_{i})\,dt_{i},$ (23)

which can be obtained from the equal-area construction (see, for instance,
Whitham), where the range of integration includes all those characteristics
that will interact with the shock. This construction is valid so long as the
driving pulse is finite and significantly shorter than a preceding quiet
period. (For it to hold at $y$ the quiet period must last a duration $yf(y)$.)

Figure 2.— Residual wave energy (normalized by $L_{p}/\omega$) available to
drive internal and external (leading or trailing) shocks in a sinusoidal wave
train, shown as a function of the inverse of the dimensionless evolution
coordinate $\tilde{y}$. Internal shocks, which deposit wave energy from an
entire wave period, have twice the energy budget of an external shock.
However, external shocks preserve much more wave energy for late stages in
their evolution.They therefore have a much greater capacity to drive strong
shocks in stellar atmospheres.

##  5\. Hydrodynamic tests

To check our analytical results, we employ FLASH4.6 (Fryxell et al. 2000)
using its HLLC Riemann solver and 5th-order reconstruction. We resolve a
spherical domain ($A=4\pi r^{2}$ for $r_{a}<r<r_{b}$) with a uniform
resolution and a static mesh, as we have found adaptive refinement boundaries
to affect wave propagation. We initialize either a power-law or polytropic
hydrostatic structure, employing a static gravitational fields to balance
forces in the background state. A train of acoustic waves is introduced into
the simulation domain by varying the fluid velocity, density, and pressure at
the inner boundary. We fix the adiabatic index to $\gamma=5/3$ throughout. The
simulation stops when the waves reach the (reflective) outer boundary.

###  5.1. Power-Law Envelope

For our light, polytropic, power-law atmosphere we consider $\rho_{0}\propto
r^{-n}$ and adopt the Keplerian gravitiy of a central point mass, requiring
$p_{0}\propto r^{-(n+1)}\propto\rho_{0}^{1+1/n}$, and $c_{0}\propto r^{-1/2}$.
We adopt $n=5/2$. Our grid spans a factor of two in radius (5.6 in density),
which we resolve with $1.4\times 10^{5}$ cells. Into this we launch a train of
sinusoidal waves ($v(r_{a})=v_{p}\sin\omega t$ for $t>0$) with
$\omega=490r_{a}/c_{0}(r_{a})$, so the initial wavelength is $r_{a}/78$, or
1840 cells. The initial amplitude is $v_{p}/c_{0}=1/100$. At this amplitude
and wavelength the linear high-frequency limit holds very well, and the
numerical resolution is sufficient. As we show in Figure 3, the simulation
agrees very precisely with our analytical prediction for the radius of shock
formation ($R_{s}=1.15r_{a}$), as well as the development of the head shock
and all subsequent internal shocks.

The comparison is equally successful when we initialize a wave of only one
sinusoidal cycle, which develops head and tail shocks but no internal shocks;
see Figure 4. (We have explored other waveforms, such as sawtooth and reverse-
sawtooth waves. These also show excellent agreement, in the linear regime,
between analytical and numerical solutions.)

Figure 3.— Top: Velocity distribution of an outward traveling sinusoidal wave
in an ambient medium with power-law structure (normalized by sound speed at
the inner edge). The distribution is shown across the top three panels. Red
dashed lines are the prediction for the maximum wave amplitude barring shock
formation. Blue and black dashed lines are the predicted post-shock values for
a forward and internal shock, respectively. Vertical yellow line is the
predicted shock formation radius. Bottom: Distribution of density, pressure,
and relative entropy at several instances in time. The last panel shows a
comparison of the predicted shock luminosity and measured normalized maximum
wave luminosities.

In the bottom panels of Figures 3 and 5 we show the local maximum value of
$L_{w}$ in units of the input peak value $L_{p}$, demonstrating that wave
luminosity is conserved in the acoustic regime and decays as expected during
shock evolution.

Figure 4.— Velocity distribution (normalized by sound speed at the inner edge)
of a single sinusoidal wave in an ambient medium with power-law structure at
several instances in time. The head (postive phase) and tail shocks (negative
phase) evolve identically in linear heory.  Figure 5.— Top and middle: Mach
number distribution across an $n=2.5$ stellar polytrope excited by outward-
traveling, small amplitude sinusoidal waves (i.e., $v/c_{0}=0.01$ at
$r/r_{a}=1$). The distribution is divided into two spatial domains. Red dashed
lines are the prediction for the maximum wave amplitude barring shock
formation. Blue and black dashed lines are the predicted post-shock values
given a forward and internal shock, respectively. Vertical yellow line is the
predicted shock formation radius. Bottom: A comparison of the predicted shock
luminosity and normalized maximum wave luminosities. Also shown are several
instances of the relative entropy profile $s/s_{0}$ where $s=p/\rho^{\gamma}$.

###  5.2. Stellar Polytrope

For a polytropic stellar model we adopt the self-gravitating Lane-Emden
solution of radius $R_{*}$, in which $p_{0}\propto\rho_{0}^{1+1/n}$; again we
choose polytropic index $n=5/2$. We arrange our grid to span radii
$0.1R_{*}-0.999R_{*}$ with $1.5\times 10^{5}$ equally-spaced cells. Over this
range $\rho_{0}$ varies by a factor of $2.6\times 10^{8}$ and $L_{{\rm max}}$
varies by a factor $2.9\times 10^{11}$. At the inner boundary we introduce
sinusoidal waves of low amplitude ($v/c_{0}=10^{-2}$) and high frequency
($R_{*}/200$ in wavelength) so that each wave is initially resolved by
approximately 800 cells. This simulation is displayed in Figure 5. Again, a
head shock and multiple internal shocks form at the predicted location
($R_{s}=0.147R_{*}$) and evolve as predicted in the linear regime.

In contrast to the power-law atmosphere, however, the sharp drop in $L_{{\rm
max}}$ near the stellar surface causes the head shock to become nonlinear. The
departure from linear theory is visible by $0.98R_{*}$, where flow downstream
of the head shock reaches $v/c_{0}\sim 0.5$, and grows from there. Linear
theory remains remarkably accurate for flow Mach numbers of order unity. New
features, such as the presence of significant post-shock velocity, become
visible as the head shock becomes strong.

We note that the internal shocks remain weak as the head shock strengthens.
Indeed, the internal shocks continue to decelerate (because $V\approx c_{0}$
for a weak shock) as the head shock begins to accelerate; for this reason, a
gap opens behind the advancing head shock. (Tail shocks differ from head
shocks in the nonlinear regime: they cannot strengthen and advance to the same
degree.) Although our simulation volume does not capture the region outside
$R_{*}$, we can be confident that the head shock develops sufficient strength
to eject matter from the stellar surface.

##  6\. Discussion

###  6.1. Shock decay and ejection: importance of the waveform

The most striking result of our analysis is the sharp difference in the decay
properties of head and tail shocks, compared to internal shocks, in the simple
model of a sinusoidal wave train composed of a finite number of cycles.
Consider the decay phase in which
$f(t_{i,d})\rightarrow[\pi\tilde{y}^{-1},2\tilde{y}^{-1/2}]$ for [internal,
head/tail] shocks, and let us define a characteristic local propagation time

$t_{\rm char}(r)\equiv Y(r)L_{{\rm max}}^{1/2}$

to eliminate $Y$. Note that $t_{\rm char}$ need not be monotonic in $r$. With
this and our other definitions, the asymptotic shock strength is given by

$\frac{v_{d}^{2}}{c_{0}^{2}}=\frac{L_{w}(t_{i,d})}{L_{{\rm
max}}}\rightarrow\left\\{\begin{array}[]{lc}\frac{\pi^{2}}{\omega^{2}t_{\rm
char}^{2}}&\mbox{internal~{}shock}\\\ \frac{4}{\omega t_{\rm
char}}\frac{L_{p}^{1/2}}{L_{{\rm
max}}(r)^{1/2}}&\mbox{head/tail~{}shock}\end{array}\right.$ (24)

Internal shocks, remarkably, do not depend in this limit on the initial wave
amplitude; this implies a universal asymptotic shock heating behaviour that
depends only on $\omega$ and local conditions, which we shall put to use in §
6.2. Internal shocks dissipate wave energy relatively rapidly, but
nevertheless strengthen in dimensionless terms in regions where $t_{\rm
char}(r)$ declines outward.

Head and tail shocks, in contrast, retain some sensitivity to the input
luminosity, and strengthen once the combination $t_{\rm char}(r)L_{{\rm
max}}(r)^{1/2}$ begins to decline. Head and tail shocks evolve differently in
the nonlinear regime; head shocks accelerate more readily to become nascent
blastwaves.

Modelling the input wave as a pure sinusoidal wave train is highly idealized,
so we consider the internal shocks and head or tail shocks as limiting cases
of what would be observed in a real star. We anticipate that a continuous but
incoherent wave source will tend to resemble the internal shocks studied here,
perhaps with a rare head- or tail-like shock arising by chance. Transient
events, on the other hand, may tend to resemble the head and tail shocks but
also produce internal shocks in their ring-down phases. We defer a full
analysis to later work.

###  6.2. Dissipation formulae for wave luminosity

Fuller & Ro (2018) present a useful differential equation for the dissipation
of wave energy at internal shocks, building on results from Landau & Lifshitz
(1959), Ulmschneider (1970), and Mihalas & Mihalas (1984). We check this
formula by examining the asymptotic behaviour of internal shocks in sinusoidal
waves: $\tilde{y}\rightarrow\pi/f$ for large $\tilde{y}$, from equation (16).
Differentiating this and using definitions for $\tilde{y}$ and $f$ to express
the result in terms of $L_{w}$, and then using $A\rho_{0}dr=dm$ (defining the
mass coordinate $m$), we find

$dL_{w}=-\frac{\gamma+1}{\pi}\omega c_{0}^{2}\left(\frac{L_{w}}{L_{{\rm
max}}}\right)^{3/2}dm.$ (25)

Note, $L_{w}$ here means $L_{w}(t_{i,d})$ – that is, it refers to the the
conditions at the shock. The phase-averaged luminosity $\bar{L}_{w}$ is the
more relevant quantity. Given that $L_{w}\propto f^{2}$,
$\bar{L}_{w}=L_{w}(t_{i,d})/3$ in the asymptotic ‘N-wave’ state; therefore

$d\bar{L}_{w}=-\frac{\gamma+1}{\pi}\sqrt{3}\omega
c_{0}^{2}\left(\frac{\bar{L}_{w}}{L_{{\rm max}}}\right)^{3/2}dm.$ (26)

The dissipation rate presented by Fuller & Ro (2018) agrees with this result,
up to an inconsistency in their equation (17): on its left-hand side, the
symbol $L_{w}$ refers to $\bar{L}_{w}$, while on its right-hand side $L_{w}$
means $L_{w}(t_{id})=\bar{L}_{w}/3$. Note also that equation (26) is accurate
with our definition of $L_{{\rm max}}$, which is double the definition used by
Fuller & Ro.

Expression (26) is convenient in form, as it can be augmented to account for
non-shock dissipation, as Fuller & Ro have done. However, it has a couple
drawbacks. It applies only to the phase of asymptotic decay, so it cannot
accurately describe the bulk of wave heating. It can only be applied in the
region $\tilde{y}>1$ where shocks exist; but even there, the initial value of
$\bar{L}_{w}$ depends on what survives the early phase of shock evolution.

For these reasons we recommend instead using the exact form of
$\bar{L}_{w}(\tilde{y})$ that results from recognizing that $\bar{L}_{w}$ is
proportional to the residual wave energy $\tilde{E}_{w}$ as shown in Figure 2:

$\frac{\bar{L}_{w}(\tilde{y})}{\bar{L}_{w}(0)}=1-\frac{\tilde{E}_{\rm
diss}(\tilde{t}_{i}(\tilde{y}))}{\pi}$ (27)

for $\tilde{y}>1$, where $\tilde{E}_{\rm diss}(\tilde{t}_{i})$ is given by
equation (18) and $\tilde{t}_{i}(\tilde{y})={\rm sinc}^{-1}(1/\tilde{y})$.

Or, for an even more convenient alternative, an approximation based on the
common asymptotic evolution of internal shocks in (eq. 24, using
$\bar{L}_{w}=L_{w}(t_{i,d})/3$):

$\bar{L}_{w}(r)\simeq\left[\bar{L}_{w}(0)^{-k}+\left(\frac{\pi^{2}L_{{\rm
max}}(r)}{3\omega^{2}t_{\rm char}(r)^{2}}\right)^{-k}\right]^{-1/k}$ (28)

which respects both limits, and remains within $10\%$ for $k=0.64$.

Both the exact and approximate expressions (eqs. 27 and 28) assume the input
wave is a pure uninterrupted sinusoid that remains in the linear, high
frequency regime, so its shock dissipation can be predicted by the results of
§ 4.1.1. For a more sophisticated approach one must start from a more
realistic input waveform. These expressions also neglect other types of
dissipation, such as radiation diffusion; we (Ro & Matzner, 2017) have found
this to be valid when $\bar{L}_{w}$ exceeds the radiation diffusion luminosity
of the background state.

###  6.3. Surface ejection by head shocks: weak limit

We wish to make a rough estimate for the mass ejected by a head shock that
strengthens in the outer stellar envelope, in the weak-wave limit for which
this occurs very close to the stellar surface ($z\equiv(R_{*}-r_{0})/R_{*}\ll
1$, where $R_{*}$ is the stellar radius). We begin with ejection by a weak
acoustic pulse, in which shock formation also occurs in the outer envelope,
and then consider the case of a central point explosion in the next
subsection.

Our analysis here is only asymptotic, however, as we neglect radiation
diffusion, which would truncate the acceleration toward $v_{\rm esc*}$, the
stellar escape velocity, for sufficiently low wave energies (Linial et al.
2020, in prep.). We also require that the wave be sufficiently radial that
lateral pressure gradients do not alter the character of the flow (Matzner et
al., 2013; Salbi et al., 2014; Linial & Sari, 2019).

The outer polytropic scalings are $A=4\pi R_{*}^{2}$,
$\rho_{0}=\rho_{0h}z^{n}$, and $c_{0}=c_{0h}z^{1/2}$, where the ‘h’ subscript
indicates the outer normalization and $z=(R_{*}-r_{0})/R_{*}$ is the relative
depth. From these we derive the external mass $m_{h}z^{n+1}$ (where
$m_{h}=4\pi R_{*}^{3}\rho_{0h}/(n+1)$), maximum luminosity $L_{{\rm
max}}=L_{{\rm max},h}z^{n+3/2}$ (where $L_{{\rm max},h}=4\pi
R_{*}^{2}\rho_{0,h}c_{0,h}^{3}$), and characteristic time $t_{\rm
char}=t_{{\rm char},h}z^{1/2}$ (where $t_{{\rm
char},h}=(\gamma+1)(2n+1)R_{*}/8c_{0,h})$).

For a reasonable approximation to the ejected mass, we first identify the
value of $z=z_{S}$ at which the shock becomes strong, in the sense that
$v_{d}(z_{S})=c_{0}(z_{S})$, using equation (24) to extrapolate the head shock
out of the weak shock regime. We determine the corresponding shock velocity
$v_{s}(z_{S})=\left(\gamma+1+\sqrt{17+\gamma(2+\gamma})\right)c_{0}(z_{S})/4$.
We then apply the self-similar scalings for a strong planar shock,
$v_{s}\propto\rho_{0}^{-\beta}$ (Sakurai, 1960) until the ejection condition
$f_{\rm sph}Cv_{s}=v_{\rm esc*}$ is met. Here $v_{\rm esc*}$ is the stellar
escape velocity, $\beta\simeq 0.19$ is the shock acceleration index, and
$C\sim 2$ is the post-shock acceleration factor for purely planar, non-
relativistic flow (Sakurai, 1960; Ro & Matzner, 2013). The factor $f_{\rm
sph}$ accounts for the difference between spherical and planar flow, which
matters for finite but small $E_{\rm in}$ because of the slow nature of post-
shock acceleration (Litvinova & Nadezhin, 1990); for non-gravitating flow
$f_{\rm sph}\simeq(1-\\{0.51,0.34\\}z^{1/3})$ for $n=\\{\frac{3}{2},3\\}$
(Matzner & McKee, 1999).

We must also allow for the possibility of wave reflection, in the case that
the high-frequency limit does not hold throughout its transition to a strong,
accelerating shock. (Reflection is not an issue in the strong regime, as
Sakurai’s solution accounts for it self-consistently.) The high-frequency
limit requires $\omega zR_{*}/c_{0}\gg 1$; evaluating this at the weak-to-
strong transition,

$\omega\frac{z_{S}R_{*}}{c_{0}(z_{S})}=\left[\frac{L_{p}}{L_{{\rm
max},h}}\left(\frac{t_{{\rm
char},h}c_{0h}}{4R_{*}}\right)^{2}\left(\frac{c_{0h}}{\omega
R_{*}}\right)^{3+2n}\right]^{\frac{1}{5+2n}}.$ (29)

As an analysis of wave reflection is beyond the scope of this paper, we
capture it using a velocity transmission prefactor $f_{\rm tm}\leq 1$, where
$f_{\rm tm}=1$ in the absence of reflection. Then

$\displaystyle\frac{m_{\rm ej}}{m_{h}}=\left[f_{\rm tm}f_{\rm
sph}K\frac{c_{0h}}{v_{\rm esc*}}\left(\frac{L_{p}}{L_{{\rm
max},h}}\frac{16}{\omega^{2}t_{\rm char,h}^{2}}\right)^{\frac{1+2\beta
n}{5+2n}}\right]^{\gamma_{p}/\beta}$ (30)
$\displaystyle\longrightarrow\left(\frac{3.47f_{\rm tm}f_{\rm
sph}\,c_{0h}}{v_{\rm esc*}}\right)^{7.18}\left(\frac{L_{p}}{L_{{\rm
max},h}}\frac{1}{\omega^{2}t_{{\rm char},h}^{2}}\right)^{1.38}$

The coefficient $K\simeq Cv_{s}(z_{S})/c_{0}(z_{S})$ depends on our
approximate matching between weak and strong limits. The evaluation on the
second line of (30) is for radiation-dominated flow in an $n=3$ polytrope
($\gamma=\gamma_{p}=4/3$).

###  6.4. Mass ejection from a weak central explosion

Our results can also be applied to the similar problem of mass ejection caused
by a explosion of small but finite strength in the stellar centre, which has
been studied by Wyman et al. (2004) and Linial, Sari, & Fuller (2020, in
prep.) for polytropes, and by Dessart et al. (2010), Owocki et al. (2019), and
Kuriyama & Shigeyama (2020) for realistic stellar models. We consider
polytropes for simplicity. The explosion is strong until it decelerates to
about $c_{0}$; for a complete polytrope, this occurs on a time $t_{\rm
dec}\simeq[2/(5c_{0c})]^{5/3}(E_{\rm in}/\rho_{0c})^{1/3}$ for an explosion of
energy $E_{\rm in}$ that is small compared to the binding energy; subscript
$c$ indicates conditions at the center, and ‘dec’ refers to the conditions at
the deceleration radius.

While there is no shock-free phase in this problem, we can nevertheless, for a
good approximation, adopt a head shock model in which the shock forms at the
deceleration radius ($r_{s}\simeq r_{\rm dec}$); the energy in the wave pulse
is $\sim E_{\rm in}/2$ (accounting for heat deposited within $r_{\rm dec}$);
and the pulse duration, or wave half-period, is $\pi/\omega\simeq t_{\rm
dec}$. Correspondingly, $L_{p}\simeq E_{\rm in}/t_{\rm dec}$.

We require the normalized coordinate $\tilde{y}(r)=Y(r)/Y_{s}=Y(r)/Y(r_{\rm
dec})$. For the denominator, we note that $Y(r)$ is almost constant in a
uniform medium: $Y(r_{\rm dec})=\Lambda(\gamma+1)(4\pi
c_{0c}\rho_{0c}^{5})^{1/2}/2$; we drop $\Lambda$, which is a logarithm of
$r_{\rm dec}$. For the numerator,
$Y(z)=2[(\gamma+1)/(2n+1)]zR_{*}/[c_{0}(z)L_{{\rm max}}(z)^{1/2}]$ in the sub-
surface region ($z\ll 1)$.

We now apply the asymptotic head shock formula
$v_{d}^{2}/c_{0}^{2}=(4/\tilde{y})L_{p}/L_{{\rm max}}$ to identify $z_{S}$ as
the depth at which the shock regains its strength ($v_{d}=c_{0}$), and proceed
to estimate $m_{\rm ej}$ as before. The full expression is too cumbersome to
write, but we note that the ejected mass scales $\propto(f_{\rm tm}f_{\rm
sph})^{\gamma_{p}/\beta}E_{\rm in}^{\alpha}$, where

$\alpha=\frac{4\gamma_{p}(1+2\beta n)}{3\beta(5+2n)}.$ (31)

The coefficient can be determined from the polytropic structure, which fixes
$(\rho_{0h}/\rho_{0c})^{1/n}=(c_{0h}/c_{0c})^{2}=(\xi\theta^{\prime})_{\xi_{n}}$
where $\theta(\xi)$ is the Lane-Emden function with outer coordinate
$\xi_{n}$. For the case $n=3/2$ and $\gamma=5/3$, we obtain

$\frac{m_{\rm ej}}{M_{*}}\simeq\left(2.23f_{\rm tm}f_{\rm
sph}\right)^{7.52}\left(\frac{E_{\rm in}}{E_{\rm bind}}\right)^{2.09}$ (32)

where $E_{\rm bind}$ is the negative total energy of the polytrope. This is
consistent with the preliminary results of Linial and collaborators (Itai
Linial, private communication, Oct. 2020) if $Kf_{\rm tm}f_{\rm
sph}\propto(E_{\rm in}/E_{\rm bind})^{0.12}$.

We suspect that the threshold energy observed by Wyman et al. (2004) is a
consequence of finite numerical resolution in that work. However, radiation
diffusion sets a true lower limit for $m_{\rm ej}$. For analysis of this
point, see Linial et al. (2020).

###  6.5. Failed Supernova Shocks

A ‘failed supernova’ is a collapsing massive star in which the stalled
protoneutron star bounce shock fails to revive and eject the envelope.
Instead, the protoneutron star is expected to collapse into a black hole since
the total accreted mass will exceed the Tolman-Oppenheimer-Volkoff limit. The
nearly instantaneous liberation of $\sim\\!0.3\,\mathrm{M}_{\odot}$ in
neutrinos from the formation of the protoneutron star causes the over-
pressurized envelope to accelerate outward at all radii, at least initially.
The expansion launches an outward-propagating spherical acoustic pulse that
steepens into a weak shock. In the wake of the shock is an expanding
rarefaction wave that informs the exterior envelope of a collapse in the
stellar core.

Red and yellow supergiant stars are well-approximated by the light polytropic
envelope structure described in § 5.1, with index $n=2.5$, over several
decades of $r$. For all positive values of $n$, there are at least two self-
similar solutions to the flow structure: one in which the weak shock becomes
strong and is described by the Sedov-Taylor solutions, and another where its
strength vanishes and it merges with the leading edge of the rarefaction wave.
For the range $2<n<3.5$, Coughlin et al. (2018) discovered a third set of
intermediate self-similar solutions involving weak shocks of a finite
strength. However, Coughlin et al. (2019) show by linear analysis that these
self-similar weak shock solutions are linearly unstable to perturbations in
the shock strength, and Ro et al. (2019) show these shocks asymptotically
approach the rarefaction and Sedov-Taylor self-similar solutions.

Why do these intermediate solutions only exist for $2<n<3.5$? The lower limit
presumably derives from the fact that, for $n<2$ the Sedov-Taylor blastwave
weakens over time, as it implies $V/c_{0}\propto t^{(n-2)/(5-n)}$. For $n\geq
3.5$, Ro et al. (2019) found that all shocks strengthen, evolving toward the
Sedov-Taylor limit. Using our solutions from §4, we find that a head shock
obeys $v_{d}/c_{0}\propto r\sqrt{A\rho_{0}c_{0}}\propto r^{(n-7/2)/4}$: it
strengthens over time when $n>3.5$. In these conditions, a weak shock tends to
outrun outward characteristics, causally disconnecting itself from the
collapsing interior.

##  7\. Conclusion

A main conclusion of this work is that shocks created by radially-propagating
sound waves – at least in the relevant limit of high frequencies and small
amplitudes – can be captured using a straightforward extension of the
classical equal-area method, either in its integral or differential form. The
future shock evolution can effectively be read from the initial waveform and
the stellar profile of density and sound speed. The method leads to improved
expressions for the shock strength and phase-averaged luminosity that should
find application in the modelling of pre-supernova outbursts.

Any detailed prediction for shock evolution, however, requires some knowledge
of how the wave is driven. This is very clear for the case we study, in which
the input waveform is a sinusoidal function that abruptly starts and stops at
its wave nodes. Such waveforms generate internal shocks, which dissipate
rapidly after shock formation, as well as head and tail shocks (if they begin
or end in a compression) that decay slowly, and that preserve significant
energy even to the stellar surface – as in the example shown in Figure 5. The
distinct potential of head shocks to drive shock-driven outbursts is a novel
aspect of our analysis.

We highlight several avenues for future study:

1\. Examination of strengthening shocks. Weak shocks, especially head shocks,
can evolve to become strong in sufficiently stratified regions like outer
stellar envelopes. There already are approximate theories that handle this
transition (e.g., Brinkley & Kirkwood, 1947; Ro, 2017) but there may exist a
smooth interpolation into the strong-shock regime, like the one used by Tan et
al. (2001) to describe trans-relativistic explosions, which would improve on
the estimate we give in § 6.3. One should also investigate the potential for
some of the wave energy to be reflected before the transition is complete, as
well as spherical effects in post-shock acceleration.

2\. Analysis of more general waveforms, such as Gaussian random noise with a
specified power spectrum. We anticipate that the dissipation profile, as well
as the rate at which head- and tail-like shocks are generated at random, will
depend on the effective bandwidth of the input spectrum.

3\. Adaptation to moving backgrounds. The current analysis is limited to
hydrostatic environments, but significant motions could be present either in
the initial state or because of previous shocks.

4\. Characterization of the driving source. For transient events such as the
pulsational pair instability this would involve a semi-analytical analysis, or
direct numerical simulation until the emerging waves are well into the high-
frequency regime. For vigorous convection this might be accomplished as
discussed above.

5\. Extension to multipole emission and multi-dimensional propagation.
Although emission from the core will tend to be of low or zero multipole
order, some driving sources, like shell convection or compact common-envelope
in-spiral, are intrinsically multipolar. Lighthill (1978) discusses how
geometrical acoustics can be used to adapt one-dimensional analyses to the
case of non-radial propagation.

6\. Extension to moving flows. This will be necessary, for instance, to test
Leung & Fuller (2020)’s scenario in which the shock dissipation of high-
frequency waves creates a low-frequency wave, which subsequently shocks.

7\. Further studies of shock-driven ejection. Models like those of Dessart et
al. (2010), Fuller (2017), and Fuller & Ro (2018) can be connected to the
specifics of the driving waveform using the tools developed here.

## Acknowledgements

We thank Itai Linial, Re’em Sari, and Jim Fuller for sharing their work prior
to publication, and Jim Fuller for clarifying questions. This work was funded
by an NSERC Discovery Grant (CDM) and by the Gordon and Betty Moore Foundation
through Grant GBMF5076 (SR).

Greenhouse gas emissions: We estimate $100$ kg CO2 equivalent from simulations
presented here, on the basis of $\sim 400$ kWh of Berkeley, CA electricity.

## References

  * Abbott et al. (2020) Abbott, R., Abbott, T. D., Abraham, S., & Acernese, F. 2020, arXiv e-prints, arXiv:2009.01190
  * Brinkley & Kirkwood (1947) Brinkley, S. R., & Kirkwood, J. G. 1947, Physical Review, 71, 606
  * Chugai et al. (2004) Chugai, N. N., Blinnikov, S. I., Cumming, R. J., et al. 2004, Monthly Notices of the Royal Astronomical Society, 352, 1213
  * Coughlin et al. (2018) Coughlin, E. R., Quataert, E., & Ro, S. 2018, ApJ, 863, 158
  * Coughlin et al. (2019) Coughlin, E. R., Ro, S., & Quataert, E. 2019, arXiv e-prints, arXiv:1901.04487
  * Courant & Friedrichs (1948) Courant, R., & Friedrichs, K. O. 1948, Supersonic flow and shock waves
  * Dessart et al. (2010) Dessart, L., Livne, E., & Waldman, R. 2010, MNRAS, 405, 2113
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53
  * Foley et al. (2007) Foley, R. J., Smith, N., Ganeshalingam, M., et al. 2007, ApJ, 657, L105
  * Friedrichs (1948) Friedrichs, K. O. 1948, Communications on Pure and Applied Mathematics, 1, 211
  * Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273
  * Fuller (2017) Fuller, J. 2017, MNRAS, 470, 1642
  * Fuller & Ro (2018) Fuller, J., & Ro, S. 2018, MNRAS, 476, 1853
  * Gal-Yam et al. (2014) Gal-Yam, A., Arcavi, I., Ofek, E. O., et al. 2014, Nature, 509, 471
  * Groh (2014) Groh, J. H. 2014, A&A, 572, L11
  * Hosseinzadeh et al. (2017) Hosseinzadeh, G., Arcavi, I., Valenti, S., et al. 2017, ApJ, 836, 158
  * Khazov et al. (2016) Khazov, D., Yaron, O., Gal-Yam, A., et al. 2016, ApJ, 818, 3
  * Kiewe et al. (2012) Kiewe, M., Gal-Yam, A., Arcavi, I., et al. 2012, ApJ, 744, 10
  * Kuriyama & Shigeyama (2020) Kuriyama, N., & Shigeyama, T. 2020, A&A, 635, A127
  * Landau & Lifshitz (1959) Landau, L., & Lifshitz, E. 1959, Fluid mechanics, A-W series in advanced physics (Pergamon Press)
  * Landau (1945) Landau, L. D. 1945, J. Phys. USSR, 9, 496
  * Leung & Fuller (2020) Leung, S.-C., & Fuller, J. 2020, ApJ, 900, 99
  * Leung et al. (2019) Leung, S.-C., Nomoto, K., & Blinnikov, S. 2019, arXiv e-prints, arXiv:1901.11136
  * Lighthill (1978) Lighthill, J. 1978, Waves in fluids (Cambridge University Press)
  * Lin & Szeri (2001) Lin, H., & Szeri, A. J. 2001, Journal of Fluid Mechanics, 431, 161
  * Linial & Sari (2019) Linial, I., & Sari, R. 2019, Physics of Fluids, 31, 097102
  * Linial et al. (2020) Linial, I., Sari, R., & Fuller, J. 2020, in prep.
  * Litvinova & Nadezhin (1990) Litvinova, I. Y., & Nadezhin, D. K. 1990, Soviet Astronomy Letters, 16, 29
  * Margutti et al. (2014) Margutti, R., Milisavljevic, D., Soderberg, A. M., et al. 2014, ApJ, 780, 21
  * Margutti et al. (2017) Margutti, R., Kamble, A., Milisavljevic, D., et al. 2017, ApJ, 835, 140
  * Matzner et al. (2013) Matzner, C. D., Levin, Y., & Ro, S. 2013, ApJ, 779, 60
  * Matzner & McKee (1999) Matzner, C. D., & McKee, C. F. 1999, ApJ, 510, 379
  * Matzner & Ro (2020) Matzner, C. D., & Ro, S. 2020, Journal of Fluid Mechanics, submitted, arXiv:2011.07180
  * Mihalas & Mihalas (1984) Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics
  * Milisavljevic et al. (2015) Milisavljevic, D., Margutti, R., Kamble, A., et al. 2015, ApJ, 815, 120
  * Moriya et al. (2014) Moriya, T. J., Maeda, K., Taddia, F., et al. 2014, MNRAS, 439, 2917
  * Ofek et al. (2016) Ofek, E. O., Cenko, S. B., Shaviv, N. J., et al. 2016, ApJ, 824, 6
  * Owocki et al. (2019) Owocki, S. P., Hirai, R., Podsiadlowski, P., & Schneider, F. R. N. 2019, MNRAS, 485, 988
  * Pastorello et al. (2007) Pastorello, A., Smartt, S. J., Mattila, S., et al. 2007, Nature, 447, 829
  * Pastorello et al. (2008a) Pastorello, A., Mattila, S., Zampieri, L., et al. 2008a, MNRAS, 389, 113
  * Pastorello et al. (2008b) —. 2008b, MNRAS, 389, 113
  * Pastorello et al. (2008c) Pastorello, A., Quimby, R. M., Smartt, S. J., et al. 2008c, MNRAS, 389, 131
  * Pastorello et al. (2013) Pastorello, A., Cappellaro, E., Inserra, C., et al. 2013, ApJ, 767, 1
  * Pastorello et al. (2016) Pastorello, A., Wang, X.-F., Ciabattari, F., et al. 2016, MNRAS, 456, 853
  * Quataert et al. (2016) Quataert, E., Fernández, R., Kasen, D., Klion, H., & Paxton, B. 2016, MNRAS, 458, 1214
  * Quataert & Shiode (2012) Quataert, E., & Shiode, J. 2012, MNRAS, 423, L92
  * Riemann (1860) Riemann, B. 1860, Abhandlungen der Königlichen Gesellschaft der Wissenschaften in Göttingen, 8, 43
  * Ro et al. (2019) Ro, S., Coughlin, E. R., & Quataert, E. 2019, The Astrophysical Journal, 878, 150
  * Ro & Matzner (2013) Ro, S., & Matzner, C. D. 2013, ApJ, 773, 79
  * Ro & Matzner (2017) —. 2017, ApJ, 841, 9
  * Ro (2017) Ro, S. S. 2017, PhD thesis, University of California, Berkeley
  * Sakurai (1960) Sakurai, A. 1960, Comm. Pure Appl. Math, 13
  * Salbi et al. (2014) Salbi, P., Matzner, C. D., Ro, S., & Levin, Y. 2014, ApJ, 790, 71
  * Shiode & Quataert (2014) Shiode, J. H., & Quataert, E. 2014, ApJ, 780, 96
  * Shivvers et al. (2016) Shivvers, I., Zheng, W. K., Mauerhan, J., et al. 2016, MNRAS, 461, 3057
  * Shivvers et al. (2017) Shivvers, I., Zheng, W., Van Dyk, S. D., et al. 2017, MNRAS, 471, 4381
  * Smith (2008) Smith, N. 2008, Nature, 455, 201
  * Smith & Arnett (2014) Smith, N., & Arnett, W. D. 2014, ApJ, 785, 82
  * Smith et al. (2011) Smith, N., Li, W., Filippenko, A. V., & Chornock, R. 2011, MNRAS, 412, 1522
  * Smith et al. (2014) Smith, N., Mauerhan, J. C., & Prieto, J. L. 2014, MNRAS, 438, 1191
  * Smith et al. (2012) Smith, N., Mauerhan, J. C., Silverman, J. M., et al. 2012, MNRAS, 426, 1905
  * Smith et al. (2018) Smith, N., Rest, A., Andrews, J. E., et al. 2018, MNRAS, 480, 1457
  * Tan et al. (2001) Tan, J. C., Matzner, C. D., & McKee, C. F. 2001, ApJ, 551, 946
  * Tyagi & Sujith (2005) Tyagi, M., & Sujith, R. 2005, Physica D: Nonlinear Phenomena, 211, 139
  * Tyagi & Sujith (2003) Tyagi, M., & Sujith, R. I. 2003, Journal of Fluid Mechanics, 492, 1
  * Ulmschneider (1970) Ulmschneider, P. 1970, Sol. Phys., 12, 403
  * Whitham (1974) Whitham, G. 1974, Linear and Nonlinear Waves, A Wiley-Interscience publication (Wiley)
  * Woosley (2017) Woosley, S. 2017, The Astrophysical Journal, 836, 244
  * Woosley et al. (2007) Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390
  * Woosley & Heger (2015) Woosley, S. E., & Heger, A. 2015, ApJ, 810, 34
  * Wyman et al. (2004) Wyman, M. C., Chernoff, D. F., & Wasserman, I. 2004, MNRAS, 354, 1053


# Eulerian formulation of the constitutive relation for an electro-magneto-
elastic material class

Deepak Kumar
Department of Mechanical Engineering
Maulana Azad National Institute of Technology Bhopal
Madhya Pradesh, India 462003
<EMAIL_ADDRESS>

###### Abstract

A novel class of electro-magneto-elastic (EME) materials comprise electro-
active and magneto-active particles in the polymer matrix that change their
elastic behavior with an applied electromagnetic field. The material response
for such a material class is usually formulated in terms of Lagrangian strain
tensor along with Lagrangian electromagnetic field vectors as “pushed forward”
to the current configuration. This letter article presents a novel formulation
of an electro-magneto-elasticity in terms of an Eulerian strain tensor and
Eulerian electromagnetic field vectors referring to the current configuration.
Such an Eulerian formulation is often favorable from both theoretical and
computational standpoints, which avoids the “pushed forward” operation to get
the current configuration. Additionally, an exercise to deduce the
constitutive relation for an EME material class available in the existing
literature from the newly proposed relation is also illustrated.

From experience, we know that two bodies of the same ‘shape’ and dimensions,
with the same distribution of external load but made of different materials,
will behave differently [1, 2, 3]. Connecting the ‘deformation’ with the
stresses became necessary in such conditions [4, 5]. Generally, such
connections are the simplified mathematical models relating the ‘deformation’
with the ‘external loads’, known as the constitutive relations for any
material class [6, 7, 8]. Such relations primarily aim to represent
approximately the material behavior without intending to model this behavior
in all the possible external conditions. In this regard, to track the changes
in elastic behavior of a smart elastomer responding to an applied
electromagnetic field, another set of constitutive relations must be derived
from the theory of electro-magneto-elasticity [9, 10]. In such relations, the
material response is usually formulated using two commonly known Lagrangean
and Eulerian strain tensors. Lagrangean and Eulerian strain tensors are
generally related through the “pushed forward” and “pushed backward”
operations. The “pushed forward” operation is a process of taking a Lagrangian
quantity and writing its Eularian analog and vice-versa for reverse. The
necessity of using Eulerian-based relations is often favorable from both
theoretical and computational standpoints, which avoid the “pushed forward”
operation to get the current configuration. In the reference configuration of
an EME material class, the electro-active and magneto-active particles are not
aligned. Rather, they are random. So, upon applying an electric field,
magnetic field, or both together, the individual particles are now aligned in
the current configuration, and the defined stress tensor will have a physical
sense. But, when we convert all the variables to the reference configuration,
that configuration becomes a pseudo configuration, not the actual one. So, to
avoid such confusion, we propose an Eulerian approach to get the direct
physical sense and the correct measure of the stresses, including all the
kinematically dependent variables.

In the current scenario, various researchers, especially the material
experimentalist and material modelers community, are very much interested in
characterizing and segregating the various newly developed smart material
classes, which fall under the category of non-elastic material polymers [11,
12]. Such non-elastic material polymers are smart polymers in which different
coupled phenomena like electromagnetostriction [9] and electro-magneto-
viscoelastic deformation [5, 13] can be observed with an electromagnetic field
of application. In view of the works on the electro-magneto-elasticity, a
limited number of works, especially classical works [14, 15, 16], exist in the
literature. Later, few authors [17, 18] had explored the advancement in
electro-magneto-elasticity by developing various micromechanical methods
investigating the deformation response of the smart elastomers. In the same
connection, the recent works [9, 10, 12, 13, 19, 20] also extended the novel
advancements in electro-magneto-elasticity in obtaining the constitutive
universal relations [9, 10] and analysing the smart actuator response with
filler particles [12, 19, 20] including the damage-induced stress softening
effect [13].

The existing works on electro-magneto-elasticity discussed above motivated us
to develop novel constitutive equations investigating the deformation response
of a novel electro-magneto-elastic (EME) material class. Inspired by them, the
current study aims to present a direct Eulerian formulation of the
constitutive relation for a novel EME material class. To the best of our
knowledge, such an explicit Eulerian formulation of the constitutive relation
in electro-magneto-elasticity is a first of its kind.

In the next, a significant contribution of this paper is formulated by
presenting a theoretical modeling framework to develop a direct Eulerian
formulation of the constitutive relation for a novel electro-magneto-elastic
(EME) material class. Consider an EME solid continua occupying an initial
stress-free configuration $\beta_{0}$ at time $t_{0}$ in the absence of any
electric E, magnetic H, or mechanical loading. Let F denote the deformation
gradient of the initial stress-free configuration $\beta_{0}$ transformation
to the current configuration $\beta$ at time $t$ subjected to any electric,
magnetic, or mechanical loading condition. At the same time, the electro-
magnetic field variables are governed by the following equations referred to,
for example, [16] as

$\displaystyle\begin{split}\mathrm{curl}\ \textbf{E}=0,\quad\mathrm{div}\
\textbf{D}=0,\quad\mathrm{curl}\ \textbf{H}=0,\quad\mathrm{div}\
\textbf{B}=0,\end{split}$ (1)

where new variables D and B denote the electric induction (or electric
displacement) and magnetic induction (or magnetic displacement) vectors,
respectively. In addition, curl and div are the curl and divergence operators
with respect to the current spacial vector. For a given continua, these new
electromagnetic field variables are defined as [21, 22]

$\displaystyle\begin{split}\textbf{D}=\epsilon_{0}\textbf{E}+\textbf{P},\quad\textbf{B}=\mu_{0}(\textbf{H}+\textbf{M}),\end{split}$
(2)

where P and M represent the polarization and magnetization densities of the
material having $\epsilon_{0}$ and $\mu_{0}$ electric permittivity and
magnetic permeability of free space, respectively. If E and H are taken as the
independent field variable then the above equation (2) gives the other
electric variable D (respectively P) when P (respectively D) is prescribed by
a constitutive equation, along with the magnetic variable B (respectively M)
when M (respectively B) is prescribed. For a given EME continua subjected to
an electromagnetic field deformation case, the Lagrangian (or “pull back” from
$\beta$ to $\beta_{0}$) version of characteristic applied electric and
magnetic field directions in the reference configuration are given by

$\displaystyle\begin{split}\textbf{E}^{l}=\textbf{F}^{T}\textbf{E},\quad\textbf{D}^{l}=\textbf{F}^{-1}\textbf{D},\quad\textbf{P}^{l}=\textbf{F}^{-1}\textbf{P},\quad\textbf{H}^{l}=\textbf{F}^{T}\textbf{H},\quad\textbf{B}^{l}=\textbf{F}^{-1}\textbf{B},\quad\textbf{M}^{l}=\textbf{F}^{-1}\textbf{M}.\end{split}$
(3)

The Eulerian (or “pull forward” from $\beta_{0}$ to $\beta$) equivalent
version of electric and magnetic field directions in the current configuration
are now obtained by differentiation of equations (3) with respect to time as

$\displaystyle\begin{split}\dot{\textbf{E}}^{l}=\textbf{L}^{T}\textbf{E}^{l},\quad\dot{\textbf{D}}^{l}=\textbf{L}^{T}\textbf{D}^{l},\quad\dot{\textbf{P}}^{l}=\textbf{L}^{T}\textbf{P}^{l},\quad\dot{\textbf{H}}^{l}=\textbf{L}^{T}\textbf{H}^{l},\quad\dot{\textbf{B}}^{l}=\textbf{L}^{T}\textbf{B}^{l},\quad\dot{\textbf{M}}^{l}=\textbf{L}^{T}\textbf{M}^{l},\end{split}$
(4)

where $\textbf{L}=\mathrm{grad}\textbf{v}=\dot{\textbf{F}}\textbf{F}^{-1}$
represents the velocity gradient tensor. Note that the electric
$\textbf{E}^{l}$ and magnetic $\textbf{H}^{l}$ may be replaced by
$-\textbf{E}^{l}$ and $-\textbf{H}^{l}$ for reverse direction, which is
equivalently defined the characteristic anisotropy. To get rid of such non-
uniqueness, the electric e and magnetic h tensors may be defined as

$\displaystyle\begin{split}\textbf{e}=\textbf{E}^{l}\otimes\textbf{E}^{l},\quad\textbf{h}=\textbf{H}^{l}\otimes\textbf{H}^{l}.\end{split}$
(5)

Using equations (3) and (4), evolution equations for the above electric e and
magnetic h tensors are given by

$\displaystyle\begin{split}\dot{\textbf{e}}=\textbf{L}\textbf{e}+\textbf{e}\textbf{L}^{T},\quad\dot{\textbf{h}}=\textbf{L}\textbf{h}+\textbf{h}\textbf{L}^{T}.\end{split}$
(6)

At the same time, mechanical deformation in the Eulerian description that is
an evolution equation similar to (6) for a defined left Cauchy–Green
deformation $\textbf{b}=\textbf{F}\textbf{F}^{T}$ in the Lagrangian
description is obtained as

$\displaystyle\begin{split}\dot{\textbf{b}}=\dot{\textbf{F}}\textbf{F}^{T}+\textbf{F}\dot{\textbf{F}}^{T}=\dot{\textbf{F}}\left(\textbf{F}^{-1}\textbf{F}\right)\textbf{F}^{T}+\textbf{F}\left(\textbf{F}^{-1}\textbf{F}\right)^{T}\dot{\textbf{F}}^{T}\quad=(\dot{\textbf{F}}\textbf{F}^{-1})(\textbf{F}\textbf{F}^{T})+(\textbf{F}\textbf{F}^{T})(\textbf{F}^{-T}\dot{\textbf{F}}^{T})=\textbf{L}\textbf{b}+\textbf{b}\textbf{L}^{T}.\end{split}$
(7)

Hence, a kinematic tensor $\dot{\textbf{b}}$ for the mechanical deformation
and the kinematic electric tensor $\dot{\textbf{e}}$ along with a magnetic
part $\dot{\textbf{h}}$ for an electromagnetic deformation are described in
the Eulerian description via evolution equations (6) and (7). Such a set of
primarily defined kinematic field variables
$(\dot{\textbf{b}},\dot{\textbf{e}},\dot{\textbf{h}})$ in the Eulerian
description are now get ready to deduce the constitutive relations for an EME
material class.

To develop the constitutive relations for an EME material class, a purely EME
deformation is assumed here by ignoring an internal dissipation of the
material. The corresponding thermodynamically consistent Clausius–Duhem
inequality-based internal dissipation $D_{int}$ vanishing for the defined EME
system of continua with the defined effective conduction current
$\textbf{J}_{e}$, effective electric current
$\textbf{E}_{e}=\textbf{E}+\textbf{v}\times\textbf{B}$, and effective
magnetization
$\textbf{M}_{e}=\mu_{0}^{-1}\textbf{M}+\textbf{v}\times\textbf{P}$ (see [23,
24, 25], for details) is given by

$\displaystyle\begin{split}D_{int}=\textbf{T}:\textbf{d}-\varrho\dot{w}+\textbf{J}_{e}.\textbf{E}_{e}-\dot{\textbf{E}}_{e}\textbf{P}-\textbf{M}_{e}.\dot{\textbf{B}}=0,\end{split}$
(8)

where T represents the Cauchy stress tensor,
$\textbf{d}=(\textbf{L}+\textbf{L}^{T})/2$ denotes the symmetric part of the
deformation rate tensor, and $\varrho$ and $w$ are the current mass density
along with the specific free energy per unit mass. By targeting the Eulerian
description of the constitutive relations for an EME material class, the
specific free energy $w(\dot{\textbf{b}},\dot{\textbf{e}},\dot{\textbf{h}})$
is assumed as a function of the defined kinematic variables. Consequently, the
rate of the free energy is given by

$\displaystyle\begin{split}\dot{w}(\dot{\textbf{b}},\dot{\textbf{e}},\dot{\textbf{h}})=\dfrac{\partial
w}{\partial\textbf{b}}:\dot{\textbf{b}}+\dfrac{\partial
w}{\partial\textbf{e}}:\dot{\textbf{e}}+\dfrac{\partial
w}{\partial\textbf{h}}:\dot{\textbf{h}}.\end{split}$ (9)

Substituting the equations (6) and (7) in the above equation (9), the rate of
the free energy expression is re-written as

$\displaystyle\begin{split}\dot{w}(\dot{\textbf{b}},\dot{\textbf{e}},\dot{\textbf{h}})=2(\textbf{Z}:\textbf{L}),\end{split}$
(10)

where

$\displaystyle\begin{split}\textbf{Z}=\dfrac{\partial
w}{\partial\textbf{b}}\textbf{b}+\dfrac{\partial
w}{\partial\textbf{e}}\textbf{e}+\dfrac{\partial
w}{\partial\textbf{h}}\textbf{h}.\end{split}$ (11)

With the use of equation (11) on (8), the internal dissipation expression
vanishing for the defined EME system of continua is re-written as

$\displaystyle\begin{split}D_{int}=\textbf{T}:\textbf{d}-\varrho\left(\textbf{Z}+\textbf{Z}^{T}\right):\textbf{d}-\varrho\left(\textbf{Z}+\textbf{Z}^{T}\right):\textbf{w}\quad+\textbf{J}_{e}.\textbf{E}_{e}-\dot{\textbf{E}}_{e}\textbf{P}-\textbf{M}_{e}.\dot{\textbf{B}}=0,\end{split}$
(12)

where $\textbf{w}=(\textbf{L}-\textbf{L}^{T})/2$ denotes the anti-symmetric
part of the deformation rate tensor, known as a spin tensor. Assuming the
components of d and w as independent, a given tensor Z is to be needed as
$\textbf{Z}=\textbf{Z}^{T}$. Using the equations (11) in (12) that is expected
to hold for real materials at all times and at every fixed point in space for
a certain class of admissible thermodynamic processes, the constitutive
equations for an EME material class are obtained as

$\displaystyle\begin{split}\textbf{T}=2\varrho\textbf{Z},\quad\textbf{P}=-\varrho\dfrac{\partial
w}{\partial\textbf{E}_{e}},\quad\textbf{M}_{e}=-\varrho\dfrac{\partial
w}{\partial\textbf{B}}.\end{split}$ (13)

The above equation represents an alternative form of the constitutive
relations formulated directly with an Eulerian description for an EME material
class in line with the literature [9, 10, 11, 12]. Utilizing the above newly
proposed relation (13), the upcoming section aimed to deduce the original form
of the relations available in the existing literature.

At last, the original form of the constitutive relations existing in the
literature is deduced from the newly developed constitutive relation (13) in
the previous section. By way of example, considering an electro-magneto-
elasticity case with an applied electric E and magnetic H fields, wherein the
free energy $w(I_{1},I_{2},...I_{9})$ depending on nine invariants given by
[9, 10]

$\displaystyle\begin{split}{I_{1}=\mathrm{tr}\textbf{b}},\quad{I_{2}=\dfrac{1}{2}[(\mathrm{tr}\textbf{b})^{2}-\mathrm{tr}(\textbf{b}^{2})]},\quad{I_{3}=\mathrm{det}\textbf{b}}\quad{I_{4}=\textbf{e}:\textbf{I}},\quad{I_{5}=\textbf{e}:\textbf{b}^{-1}},\quad{I_{6}=\textbf{e}:\textbf{b}^{-2}},\\\
{I_{7}=\textbf{h}:\textbf{I}},\quad{I_{8}=\textbf{h}:\textbf{b}^{-1}},\quad{I_{9}=\textbf{h}:\textbf{b}^{-2}}.\end{split}$
(14)

The first three invariants are the principal ones that correspond to the
mechanical deformation, while the last six invariants correspond to the
electromagnetic field effect. It is to be noted that in the Eulerian
description, the invariants $I_{4}$ to $I_{9}$ linked with electromagnetic
effects are generally functions of both deformation tensor b and the
corresponding field tensor either e or h. The corresponding nontrivial
derivatives of such invariants are obtained as

$\displaystyle\begin{split}\dfrac{\partial
I_{1}}{\partial\textbf{b}}=\textbf{I},\quad\dfrac{\partial
I_{2}}{\partial\textbf{b}}=I_{1}\textbf{I}-\textbf{b},\quad\dfrac{\partial
I_{3}}{\partial\textbf{b}}=I_{3}\textbf{b}^{-1},\quad\dfrac{\partial
I_{4}}{\partial\textbf{e}}=\textbf{I},\quad\dfrac{\partial
I_{5}}{\partial\textbf{e}}=\textbf{b}^{-1},\quad\dfrac{\partial
I_{6}}{\partial\textbf{e}}=\textbf{b}^{-2},\\\ \dfrac{\partial
I_{7}}{\partial\textbf{h}}=\textbf{I},\quad\dfrac{\partial
I_{8}}{\partial\textbf{h}}=\textbf{b}^{-1},\quad\dfrac{\partial
I_{9}}{\partial\textbf{h}}=\textbf{b}^{-2}.\end{split}$ (15)

For a given invariants-dependent free energy function
$w(I_{1},I_{2},...I_{9})$ along with the equation (10) yields the following
expression

$\displaystyle\begin{split}\dfrac{\partial
w}{\partial\textbf{b}}=\sum_{n=1}^{9}w_{n}\dfrac{\partial
I_{n}}{\partial\textbf{b}},\quad\dfrac{\partial
w}{\partial\textbf{e}}=\sum_{n=4}^{6}w_{n}\dfrac{\partial
I_{n}}{\partial\textbf{e}},\quad\dfrac{\partial
w}{\partial\textbf{h}}=\sum_{n=7}^{9}w_{n}\dfrac{\partial
I_{n}}{\partial\textbf{h}},\end{split}$ (16)

where $w_{n}=\dfrac{\partial w}{\partial I_{n}}$. Substituting the
corresponding nontrivial derivatives expressions from (15) in (16), the above
expressions (16) may be re-written as

$\displaystyle\begin{split}\dfrac{\partial
w}{\partial\textbf{b}}=(w_{1}+w_{2}I_{1})\textbf{I}-w_{2}\textbf{b}+w_{3}I_{3}\textbf{b}^{-1}-w_{5}\textbf{b}^{-2}\textbf{e}-w_{6}\textbf{b}^{-3}\textbf{e}-w_{8}\textbf{b}^{-2}\textbf{h}-w_{9}\textbf{b}^{-3}\textbf{h},\\\
\dfrac{\partial
w}{\partial\textbf{e}}=w_{4}\textbf{I}+w_{5}\textbf{b}^{-1}+w_{6}\textbf{b}^{-2},\quad\dfrac{\partial
w}{\partial\textbf{h}}=w_{7}\textbf{I}+w_{8}\textbf{b}^{-1}+w_{9}\textbf{b}^{-2},\end{split}$
(17)

Further, substituting the equations (11), (16), and (17) in a newly developed
form of the constitutive relations (13) formulated directly with an Eulerian
description, we get

$\displaystyle\begin{split}\textbf{T}=2\varrho[(w_{1}+w_{2}I_{1})\textbf{b}-w_{2}\textbf{b}^{2}+w_{3}I_{3}\textbf{I}+w_{4}\textbf{e}+w_{5}(\textbf{b}^{-1}\textbf{e}+\textbf{e}\textbf{b}^{-1})+w_{6}(\textbf{b}^{-2}\textbf{e}+\textbf{e}\textbf{b}^{-2})\\\
+w_{8}(\textbf{b}^{-1}\textbf{h}+\textbf{h}\textbf{b}^{-1})+w_{9}(\textbf{b}^{-2}\textbf{h}+\textbf{h}\textbf{b}^{-2})].\end{split}$
(18)

Finally, the above-obtained relation is an equivalent original form of the
constitutive relation to evaluate the stress response compared with the
existing literature [9, 10, 11, 12].

In the summary, a novel formulation of an electro-magneto-elasticity in terms
of an Eulerian strain tensor and Eulerian electromagnetic field vectors
referring to the current configuration is developed. To develop the same, a
thermodynamically consistent continuum mechanics-based modeling approach is
adopted. A newly proposed relation (13) is further exercised in obtaining the
original form of the constitutive relation (18) available in the existing
literature. In short, it is emphasized in the current work that Eulerian-based
relations are often favorable compared to Lagrangian from both theoretical and
computational standpoints, which avoid the “pushed forward” operation to get
the current configuration.

## Acknowledgement

The author is grateful to Dr. Somnath Sarangi, Associate Professor in the
Department of Mechanical Engineering, IIT Patna, India, for the fruitful
discussions on the topic.

## References

  * [1] D. Kumar, S. Ghosh, S. Roy, S. Santapuri, Modeling and analysis of an electro-pneumatic braided muscle actuator, Journal of Intelligent Material Systems and Structures 32 (4) (2021) 399–409.
  * [2] A. Khurana, A. Kumar, S. K. Raut, A. K. Sharma, M. M. Joglekar, Effect of viscoelasticity on the nonlinear dynamic behavior of dielectric elastomer minimum energy structures, International Journal of Solids and Structures 208 (2021) 141–153.
  * [3] A. K. Sharma, A. Khurana, M. M. Joglekar, A finite element model for investigating the thermo-electro-mechanical response of inhomogeneously deforming dielectric elastomer actuators, European Journal of Computational Mechanics (2021) 387–408.
  * [4] D. Kumar, S. Sarangi, Constitutive modeling of an electro-magneto-rheological fluid, Scientific Reports 12 (1) (2022) 1–12.
  * [5] D. Kumar, S. Sarangi, A novel class of universal relation for incompressible isotropic electro-viscoelastic materials, Mechanics Research Communications 117 (2021) 103784.
  * [6] S. Gour, D. Kumar, A. Khurana, Constitutive modeling for the tear fracture of artificial tissues in human-like soft robots, European Journal of Mechanics-A/Solids (2022) 104672.
  * [7] S. Gour, D. Kumar, V. Yadav, Constitutive modeling for the tear fracture of rubber with filler particles, Sādhanā 47 (2) (2022) 1–9.
  * [8] S. K. Behera, D. Kumar, S. Sarangi, Modeling of electro–viscoelastic dielectric elastomer: A continuum mechanics approach, European Journal of Mechanics-A/Solids 90 (2021) 104369.
  * [9] D. Kumar, S. Sarangi, Electro-magnetostriction under large deformation: modeling with experimental validation, Mechanics of Materials 128 (2019) 1–10.
  * [10] D. Kumar, S. Sarangi, P. Saxena, Universal relations in coupled electro-magneto-elasticity, Mechanics of Materials 143 (2020) 103308.
  * [11] D. Kumar, S. Sarangi, Instability analysis of an electro-magneto-elastic actuator: A continuum mechanics approach, AIP Advances 8 (11) (2018) 115314.
  * [12] D. Kumar, V. Yadav, S. Sarangi, Modeling and analysis of an electro-magneto-elastic rotating cylindrical tube actuator, Journal of Intelligent Material Systems and Structures (2022) 1045389X211072188.
  * [13] S. K. Behera, D. Kumar, S. Sarangi, Constitutive modeling of damage-induced stress softening in electro-magneto-viscoelastic materials, Mechanics of Materials (2022) 104348.
  * [14] J. Alblas, Electro-magneto-elasticity, in: Topics in applied continuum mechanics, Springer, 1974, pp. 71–114.
  * [15] L. Bracke, R. Van Vliet, A broadband magneto-electric transducer using a composite material, International Journal of Electronics Theoretical and Experimental 51 (3) (1981) 255–262.
  * [16] A. Kovetz, Electromagnetic theory, Vol. 975, Oxford University Press Oxford, 2000\.
  * [17] G. P. Carman, K.-S. Cheung, D. Wang, Micro-mechanical model of a composite containing a conservative nonlinear electro-magneto-thermo-mechanical material, Journal of intelligent material systems and structures 6 (5) (1995) 691–699.
  * [18] J. Aboudi, Micromechanical analysis of fully coupled electro-magneto-thermo-elastic multiphase composites, Smart materials and structures 10 (5) (2001) 867.
  * [19] A. Khurana, D. Kumar, A. K. Sharma, M. M. Joglekar, Nonlinear oscillations of particle-reinforced electro-magneto-viscoelastomer actuators, Journal of Applied Mechanics 88 (12) (2021).
  * [20] A. Khurana, D. Kumar, A. K. Sharma, M. Joglekar, Static and dynamic instability modeling of electro-magneto-active polymers with various entanglements and crosslinks, International Journal of Non-Linear Mechanics 139 (2022) 103865.
  * [21] A. Chełkowski, Dielectric physics, Vol. 9, Elsevier Science & Technology, 1980\.
  * [22] E. Callen, Magnetostriction, Journal of Applied Physics 39 (2) (1968) 519–527.
  * [23] R. Bustamante, Mathematical modelling of non-linear magneto-and electro-active rubber-like materials, Ph.D. thesis, University of Glasgow (2007).
  * [24] D. Kumar, Mathematical modeling of smart materials-a continuum mechanics approach, Ph.D. thesis, IIT Patna (2019).
  * [25] D. Kumar, S. Sarangi, R. Bhattacharyya, Universal relations in nonlinear electro-magneto-elasticity, Archive of Applied Mechanics 90 (7) (2020) 1643–1657.