# On the bipolarity of Wolf-Rayet nebulae

D. M.-A. Meyer
Universität Potsdam, Institut für Physik und Astronomie, Karl-Liebknecht-
Strasse 24/25, 14476 Potsdam, Germany
E-mail<EMAIL_ADDRESS>
(Received; accepted)

###### Abstract

Wolf-Rayet stars are amongst the rarest but also most intriguing massive
stars. Their extreme stellar winds induce famous multi-wavelength
circumstellar gas nebulae of various morphologies, spanning from circles and
rings to bipolar shapes. This study is devoted to the investigation of the
formation of young, asymmetric Wolf-Rayet gas nebulae and we present a
2.5-dimensional magneto-hydrodynamical toy model for the simulation of Wolf-
Rayet gas nebulae generated by wind-wind interaction. Our method accounts for
stellar wind asymmetries, rotation, magnetisation, evolution and mixing of
materials. It is found that the morphology of the Wolf-Rayet nebulae of blue
supergiant ancestors is tightly related to the wind geometry and to the
stellar phase transition time interval, generating either a broadened peanut-
like or a collimated jet-like gas nebula. Radiative transfer calculations of
our Wolf-Rayet nebulae for dust infrared emission at $24\,\mu\rm m$ show that
the projected diffuse emission can appear as oblate, bipolar, ellipsoidal or
ring structures. Important projection effects are at work in shaping observed
Wolf-Rayet nebulae. This might call a revision of the various classifications
of Wolf-Rayet shells, which are mostly based on their observed shape.
Particularly, our models question the possibility of producing pre-Wolf-Rayet
wind asymmetries, responsible for bipolar nebulae like NGC 6888, within the
single red supergiant evolution channel scenario. We propose that bipolar
Wolf-Rayet nebulae can only be formed within the red supergiant scenario by
multiple/merged massive stellar systems, or by single high-mass stars
undergoing additional, e.g. blue supergiant, evolutionary stages prior to the
Wolf-Rayet phase.

###### keywords:

methods: MHD – radiative transfer – stars: massive – stars: circumstellar
matter.

##  1 Introduction

Formed out of the gravitational collapse of opaque pre-stellar cores, high-
mass stellar objects, i.e. with mass $\geq 8\,\rm M_{\odot}$, are cosmic
regulators of the cycle of matter at work in the interstellar medium (ISM) of
galaxies (Maeder, 2009; Langer, 2012). One of their characteristic features is
the fast and strong wind blown out from their surface and powered by the
complex nuclear reactions at work in their interior. Throughout their lives,
massive stars first experience a rather long, main-sequence phase
corresponding to the burning of hydrogen in their cores (Ekström et al.,
2012). It is characterised by the release of diluted, hot and supersonic winds
from the stellar surface, and, once their hydrogen is exhausted, stellar
evolution is triggered. The entire evolution of massive stars is fixed by
their zero-age main-sequence properties, such as their mass, angular velocity
and chemical composition, which uniquely determines their evolution and
ultimate fate (Woosley et al., 2002; Vink, 2006; Brott et al., 2011; Szécsi et
al., 2020). These initial stellar characteristics rule the time-dependence of
the stellar surface properties, as well as the number and duration of the
various evolutionary phases, up to their death, either as a core-collapse
supernova leaving behind a plerionic supernova remnant, or directly collapsing
as a black hole (Chevalier, 1977; Weiler & Sramek, 1988; Woosley et al., 2002;
Woosley & Bloom, 2006; Smartt, 2009; Müller & Janka, 2015; Gabler et al.,
2021).

Amongst the many possible evolutionary histories massive stars can undergo,
high-mass stellar objects can evolve to the so-called Wolf-Rayet stage. This
phase is characterised by fast ($\sim 1000$-$5000\,\rm km\,\rm s^{-1}$) and
dense stellar winds expelled at mass-loss rates $\geq 10^{-5}\,\rm
M_{\odot}\,\rm yr^{-1}$ and enriched in C, N and O elements (Hamann et al.,
2006; Sander et al., 2012; Bestenlehner et al., 2014). The precise manner a
massive stellar object acquires the spectroscopic properties of a Wolf-Rayet
star is not completely understood, however, a few typical evolutionary paths
have been identified (Crowther, 2007). The most common channel reposes in an
evolution from the main-sequence to the Wolf-Rayet phase through a red
supergiant phase, and it mostly concerns Galactic stars with masses $\geq
20\,\rm M_{\odot}$. Heavier zero-age main-sequence stars (up to $\sim 60\,\rm
M_{\odot}$) have been shown theoretically to experience a much more complex
evolution, including several successive luminous blue variable phases, during
which the star blows, sometimes eruptively, blue supergiant material (Groh et
al., 2014).

The interaction of the stellar wind of massive stars with the ISM produces gas
nebulae witnessing the release of energy, momentum and heavy chemical elements
which enriched their local ambient medium (Langer, 2012). The classical
picture for stellar wind-ISM interaction is that of Weaver et al. (1977) and
it consists of a several parsec-scale structured circumstellar bubbles
surrounding massive stars. The circumstellar medium of Wolf-Rayet stars is
different, as many of them display unusual smaller-scale ring-like gas nebulae
detected in H$\alpha$ and infrared, both in the Milky Way, see Chu & Treffers
(1981) and Treffers & Chu (1982), and in the Magellanic Clouds (Chu, 1982b,
1983; Dopita et al., 1994; Weis et al., 1997; Hung et al., 2021),
respectively. A classification of circumstellar structures around Wolf-Rayet
stars arose from that studies (Chu, 1981, 1982b, 1982a; Treffers & Chu, 1982;
Chu et al., 1983; Toalá et al., 2015). Several follow-up surveys deepened and
completed these observations in the southern hemisphere, see Marston et al.
(1994) and Marston et al. (1994). Wolf-Rayet nebulae entered the X-rays domain
with Asca (Wrigge et al., 2005) Rosat and Chandra (Guerrero & Chu, 2008), XMM-
Newton (Toalá et al., 2012) observations from their (shocked) stellar winds
witnessing emission of optically-thin high temperature plasma. More recently,
a radio non-thermal synchrotron counterpart of the stellar wind bubble of a
Wolf-Rayet ring nebula indicated that these circumstellar nebulae are the site
of particles accelelration (Prajapati et al., 2019). This accumulation of
observations motivates the present numerical efforts aiming at understanding
the peculiar vicinity of young Wolf-Rayet stars.

Numerical simulations rapidly became an efficient tool to understand the
different physical processes at work in the circumstellar medium of massive
stars (Comerón & Kaper, 1998; Freyer et al., 2003; Dwarkadas, 2005; Freyer et
al., 2006; Toalá & Arthur, 2011; van Marle et al., 2011; van Marle et al.,
2015), to constrain stellar evolution models/probe local ISM conditions
(Mackey et al., 2012; Mackey et al., 2014), and to understand the pre-
supernova surroundings of dying massive stars, with which supernovae shock
waves subsequently interact (Ciotti & D’Ercole, 1989; Chevalier & Liang, 1989;
van Marle et al., 2005, 2006; Meyer et al., 2021). Particularly, the gas
nebulae formed around Wolf-Rayet stars have been studied numerically as
evolved winds, colliding with the material of a larger-scale main-sequence
wind-blown bubble (Brighenti & D’Ercole, 1997; Dwarkadas, 2007; van Marle et
al., 2007; van Marle et al., 2011).

It has recently been shown in Meyer et al. (2020), on the basis of magneto-
hydrodynamical (MHD) simulations, that (i) the circumstellar rings appearing
to be comoving with the fastest-moving Wolf-Rayet stars are a trace of their
very high initial masses ($\geq 35$ to $\sim 60\,\rm M_{\odot}$) inducing
complex evolutionary histories, and that (ii) their large-scale surroundings
are main-sequence stellar wind bow shocks fainter than the inner rings which
became unobservable as the driving star runs into diluted media (Gvaramadze et
al., 2010). This solved the apparent missing bow shock problem around some
high-latitude, Galactic, very fast-moving Wolf-Rayet stars like WR124,
surrounded by its compact gas nebula M1-67 (Chu, 1981; van der Sluys & Lamers,
2003; Marchenko et al., 2010; Toalá et al., 2018) or WR71 (Faherty et al.,
2014; Flagey et al., 2014). If we can now qualitatively explain the formation
of quasi-spherical rings around very fast-moving Wolf-Rayet stars, detailed
observations of these gas nebulae reveal much more complex structures whose
formation are far from being understood. Which mechanisms distinguish between
the production of ring-like shells from ovoidal and/or other complex shapes ?
What can we learn regarding the past evolution of massive stars with nebulae
of such morphologies ?

In this study, we continue our investigations of the shaping of circumstellar
Wolf-Rayet nebulae started in Meyer et al. (2020). We concentrate on the
production of bipolar asymmetries and non-spherical morphologies observed in
the vicinity of some young Wolf-Rayet stars (Toalá et al., 2015). The starting
point of this study is the work of Brighenti & D’Ercole (1997), in which the
development of young aspherical Wolf-Rayet nebulae is examined using an
hydrodynamical wind-wind interaction model. Their simulations focus on the
particular scenario of a red supergiant star, with asymmetric stellar wind
properties, evolving to a Wolf-Rayet stage of symmetric stellar wind
properties. In the same spirit, we build a magneto-hydrodynamical toy model
for rotating blue supergiant stars with asymmetric stellar winds and evolving
to the Wolf-Rayet phase. The effects of asymmetries in the blue supergiant
wind are investigated using the recipe of Raga et al. (2008). We perform near-
infrared synthetic observables of the corresponding bipolar stellar wind
nebulae, in order to discuss their respective emission properties and
facilitate qualitative comparison to observations available in the literature.

Our study is organised as follows. Firstly, we present the numerical methods
utilised for the MHD simulations of Wolf-Rayet gas nebulae evolving after
having undergone a previous blue supergiant phase in Section 2. We show our
results regarding to the dynamical evolution of bipolar MHD Wolf-Rayet nebulae
generated by wind-wind interaction in Section 3. The results are further
discussed in Section 4, and finally, we present our conclusions in Section 5.

##  2 Method

This section presents the numerical methods used in this project, together
with the initial conditions and parametrised boundary conditions of the
simulations. Last, we introduce the reader to the models performed in this
study.

###  2.1 Governing equations

We set our problem in the frame of non-ideal magneto-hydrodynamics, described
by the following series of equations,

$\frac{\partial\rho}{\partial
t}+\mn@boldsymbol{\nabla}\cdot\big{(}\rho\mn@boldsymbol{v})=0,$ (1)
$\frac{\partial\mn@boldsymbol{m}}{\partial
t}+\mn@boldsymbol{\nabla}\cdot\Big{(}\mn@boldsymbol{m}{\color[rgb]{0,0,0}\otimes}\mn@boldsymbol{v}{\color[rgb]{0,0,0}-}\mn@boldsymbol{B}{\color[rgb]{0,0,0}\otimes}\mn@boldsymbol{B}+\mn@boldsymbol{\hat{\mathbf{I}}}p_{\rm
t}\Big{)}=\mn@boldsymbol{0},$ (2) $\frac{\partial E}{\partial
t}+\mn@boldsymbol{\nabla}\cdot\Big{(}(E+p_{\rm
t})\mn@boldsymbol{v}-\mn@boldsymbol{B}(\mn@boldsymbol{v}\cdot\mn@boldsymbol{B})\Big{)}=\Phi(T,\rho),$
(3)

and,

$\frac{\partial\mn@boldsymbol{B}}{\partial
t}+\mn@boldsymbol{\nabla}\cdot\Big{(}\mn@boldsymbol{v}{\color[rgb]{0,0,0}\otimes}\mn@boldsymbol{B}-\mn@boldsymbol{B}{\color[rgb]{0,0,0}\otimes}\mn@boldsymbol{v}\Big{)}=\mn@boldsymbol{0},$
(4)

where $\mn@boldsymbol{B}$ represents the magnetic field vector, $\rho$ stands
for the mass density,

$\mn@boldsymbol{m}=\rho\mn@boldsymbol{v},$ (5)

is the linear momentum vector, $\mn@boldsymbol{\hat{\mathbf{I}}}$ the identity
matrix, $p_{\rm t}$ is the total pressure and $v$ the gas velocity,
respectively. The total energy of the system reads,

$E=\frac{p}{(\gamma-1)}+\frac{\mn@boldsymbol{m}\cdot\mn@boldsymbol{m}}{2\rho}+\frac{\mn@boldsymbol{B}\cdot\mn@boldsymbol{B}}{2},$
(6)

and the sound speed is defined as,

$c_{\rm s}=\sqrt{\frac{\gamma p}{\rho}},$ (7)

which closes the above system of equations, and where $\gamma=5/3$ is the
adiabatic index.

Radiative cooling and heating of the gas,

$\mathit{\Phi}(T,\rho)=n_{\mathrm{H}}\Gamma(T)-n^{2}_{\mathrm{H}}\Lambda(T),$
(8)

are explicitly included via optically-thin processes for loss $\Lambda(T)$ and
gain $\Gamma(T)$ in the the source term, where,

$T=\mu\frac{m_{\mathrm{H}}}{k_{\rm{B}}}\frac{p}{\rho},$ (9)

is the gas temperature, and using the laws of photoionized gas described in
great detail in Meyer et al. (2021). Hence, our method ignores any explicit
treatment of the radiation pressure from the photon field of the central star.
This is correct for cool blue supergiant stars, while we neglect the detailed
position of the ionization front generated by the young Wolf-Rayet stars, as
well as its effects on the shell instabilities (Toalá & Arthur, 2011).

A passive scalar tracer is included in the stellar wind,

$\frac{\partial(\rho Q)}{\partial
t}+\mn@boldsymbol{\nabla}\cdot(\mn@boldsymbol{v}\rho Q)=0,$ (10)

to follow the advection of the new-born Wolf-Rayet circumstellar nebula, as
compared to the previous blue supergiant wind. This system of equations is
solved by use of the eight-wave algorithm (Powell, 1997), within a second-
order Runge-Kutta with parabolic reconstruction of the variables between
neighbouring cells together with the HLL Riemann solver (Harten et al., 1983).
This unsplit scheme ensures that the magnetic field is divergence-free
everywhere in the computational domain. The simulations timesteps are
controlled by the Courant-Friedrich-Levi condition that we set to $C_{\rm
cfl}=0.1$.

Table 1:  Stellar surface parameters used in our two-winds models.  $\mathrm{Quantities}$ | $\mathrm{Symbol}$ | $\mathrm{Supergiant}$ |  $\mathrm{Wolf}$$-$$\mathrm{Rayet}$
---|---|---|---
Accretion rate | $\dot{M}\,(\rm M_{\odot}\,\rm yr^{-1})$ | $10^{-6}$ | $10^{-4.3}$
Wind velocity | $v_{\rm w}\,(\rm km\,\rm s^{-1})$ | 500 | 1900
Angular velocity | $v_{\rm rot}\,(\rm km\,\rm s^{-1})$ | 60 | 10
Surface field | $B_{\star}\,(\rm G)$ | 1 | 100
Stellar radius | $R_{\star}\,(\rm R_{\odot})$ | 20 | 2.3

###  2.2 Initial conditions

We carry out MHD simulations with the pluto code (Mignone et al., 2007, 2012;
Vaidya et al., 2018). The calculations are performed using a 2.5-dimensional
spherical coordinate system ($O$;$r$,$\theta$,$\phi$) of origin $O$, with $r$
the radial direction, $\theta$ the poloidal coordinate and $\phi$ the toroidal
component. The coordinate system $[r_{\rm in},r_{\rm
out}]\times[0,\pi]\times[0,2\pi]$ is mapped with a mesh that is uniform along
the polar and toroidal directions, while it expands logarithmically along the
radial direction. It permits to reach high spatial resolutions close to the
stellar wind boundary $r_{\rm in}=0.02\,\rm pc$, while reducing the total
number of grid zones in the simulation domain and reducing the computational
costs. In total, we use $200\times\,200\times\,1$ grid zones. Outflow boundary
conditions are assigned at the outer boundary $r_{\rm out}=10\,\rm pc$,
reflective boundary conditions are imposed along the symmetry axis $\theta=0$
and $\theta=\pi$, while periodic boundary conditions are used at the borders
$\phi=0$ and $\phi=2\pi$.  Since the inner boundary is much larger than the
radius of the photosphere ($r_{\rm in}\gg R_{\star}$), stellar gravity and
wind acceleration mechanisms such as coronal heating and radiation pressure,
are neglected in the computational domain. Therefore, the acceleration of the
stellar wind to the terminal velocity is not modelled in the simulations.

The simulations are initialised with the asymmetric stellar wind of a blue
supergiant star. The adopted latitude-dependence of the wind density is that
developed in the context of the astrosphere of the runaway asymptotic giant
branch star Mira (Raga et al., 2008). It reads,

$\rho_{w}(r,\theta)=\frac{\dot{M}}{4\pi r^{2}v_{\rm
w}}\rightarrow\frac{A}{r^{2}}f(\theta),$ (11)

with $\dot{M}$ the stellar mass-loss rate, $v_{\rm w}$ the stellar wind
velocity and $r$ the radial distance to the central star. The function,

$f(\theta)=\xi-(\xi-1)|\cos(\theta)|^{p},$ (12)

measures the wind anisotropy, with $p$ a parameter determining the flattening
degree of the density towards the equator and $\xi$ is the equator-to-pole
density ratio, respectively. The scaling factor $A$ ensures that the latitude-
dependence in the stellar wind conserves the amount of mass lost by the star
per unit time. It is evaluated by writing the mass conservation at the stellar
surface,

$\dot{M}=4\pi
r^{2}\int_{0}^{\pi/2}\rho_{w}(r,\theta)v_{w}(r,\theta)\sin(\theta)d\theta,$
(13)

yielding,

$\dot{M}=4\pi Av_{\rm w}\int_{0}^{\pi/2}|f(\theta)|^{p}\sin(\theta)d\theta,$
(14)

where $v_{\rm w}$ is the stellar wind terminal velocity, and, finally, one
obtains,

$A=\frac{\dot{M}}{4\pi v_{\rm w}}\frac{(\xi-1)^{2}}{8\xi^{5/2}-20\xi+12},$
(15)

using $p=1/2$. We chose this prescriptions for stellar wind asymmetries since
it has already been used in the context of massive stars (Fang et al., 2017),
although other recipes exist (Frank et al., 1995; Blondin et al., 1996).

Table 2:  List of models in our study. The table informs on the simulation labels, whether the blue supergiant wind is symmetric or not, the blue supergiant to Wolf-Rayet phase transition timescale $\Delta t$ (in $\rm yr$) and the exponent $\beta$ controlling the variations of the different quantities during the phase transition, the flattening degree of the blue supergiant wind towards the equator $p$ and the equator-to-pole density ratio $\xi$, respectively. The last column gives the general purpose of each simulation models.  ${\rm{Model}}$ | $\mathrm{Asymmetric}\,\mathrm{wind}$ | $\mathrm{\Delta t}\,(\rm yr)$ | $\beta$ | $p$ | $\xi$ | $\mathrm{Purpose}$
---|---|---|---|---|---|---
Run-Base | yes | $0$ | $4$ | $1/4$ | $20$ | Baseline model
Run-Sym | no | $0$ | $-$ | $-$ | $-$ | Model with isotropic progenitor wind
Run-$\Delta$t10kyr | yes | $10^{4}$ | $4$ | $1/4$ | $20$ | Model with $\Delta$t10kyr phase transition
Run-$\xi$20 | yes | $10^{4}$ | $4$ | $1/4$ | $20$ | Model with phase transition and lower pole-to-equator density contrast
Run-$\xi$100 | yes | $10^{4}$ | $4$ | $3/4$ | $100$ | Model with phase transition and higher pole-to-equator density contrast
Run-$\beta$4 | yes | $5\times 10^{4}$ | $4$ | $1/4$ | $20$ | Model with quadratic interpolation between winds during the phase transition
Run-$\beta$2 | yes | $5\times 10^{4}$ | $4$ | $1/4$ | $20$ | Model with linear interpolation between winds during the phase transition

Additionally, Raga et al. (2008) prescribe stellar wind flow assigned at the
wind boundary $r_{\rm in}$ to account for the asymmetry function as,

$v(r_{\rm in},\theta)=\frac{v_{\rm w}}{\sqrt{f(\theta)}},$ (16)

so that the ram pressure of the wind of the inner boundary,

$\rho_{w}(r_{\rm in},\theta)v(r_{\rm in},\theta)^{2}=\frac{A}{r_{\rm
in}^{2}}f(\theta)\Big{(}\frac{v_{\rm
w}}{\sqrt{f(\theta)}}\Big{)}^{2}\propto\Big{(}\frac{v_{\rm w}}{r_{\rm
in}}\Big{)}^{2},$ (17)

remains isotropic. Stellar rotation is included in the model by considering
the following latitude-dependent equatorial rotation velocity,

$v_{\phi}(\theta)=v_{\rm rot}\sin(\theta),$ (18)

where $v_{\rm rot}$ is the star’s equatorial rotation speed. Last, the
magnetisation of the stellar winds is treated as a Parker spiral. It is a
split monopole,

$B_{\rm r}(r)=B_{\star}\Big{(}\frac{R_{\star}}{r}\Big{)}^{2},$ (19)

where $R_{\star}$ is the stellar radius, $B_{\star}$ the stellar surface
magnetic field, plus a toroidal component,

$B_{\phi}(r)=B_{\rm r}(r)\Big{(}\frac{v_{\phi}(\theta)}{v_{\rm
w}}\Big{)}\Big{(}\frac{r}{R_{\star}}-1\Big{)},$ (20)

which relies on both the stellar surface magnetic field strength and on the
stellar rotation velocity $v_{\rm rot}$. This so-called Parker stellar wind is
a parametrisation for magnetised low-mass stars which has been developed and
widely-used in studies devoted to the heliosphere (Parker, 1958; Pogorelov &
Semenov, 1997; Pogorelov & Matsuda, 2000; Pogorelov et al., 2004). It has been
adapted to the intermediate stellar mass regime in works investigating the
shaping of planetary nebula (Chevalier & Luo, 1994; Rozyczka & Franco, 1996;
García-Segura et al., 2018, 2020) and stellar wind bow shocks (Herbst et al.,
2020; Scherer et al., 2020; Meyer et al., 2021), respectively.

###  2.3 Time-dependent stellar boundary conditions

At time $t_{\rm wr}-\Delta t$ the blue supergiant to Wolf-Rayet phase begins
to take place, controlled by several quantities, namely its duration $\Delta
t$ and the manner the surface stellar properties evolve, determined by the
quintuplet ($\dot{M}$,$v_{\rm w}$,$v_{\rm rot}$, $B_{\star}$,$R_{\star}$). The
stellar wind velocity is modulated following (Brighenti & D’Ercole, 1997),

$\displaystyle v_{\rm w}(r_{\rm in},t)=\Bigg{\\{}\begin{array}[]{lll}v_{\rm
w}^{\rm bsg}/\sqrt{f(\theta)}&\mathrm{if}&t\leq t_{\rm wr}-\Delta t,\\\ v_{\rm
w}^{\rm bsg}/\sqrt{f(\theta)}+\Delta vF&\mathrm{if}&t_{\rm wr}-\Delta
t<tandt\leq t_{\rm wr},\\\ v_{\rm w}^{\rm wr}&\mathrm{if}&t_{\rm wr}<t,\\\
\end{array}$ (24)

with $\Delta v=|v_{\star}^{\rm wr}-v_{\star}^{\rm bsg}|$. In the above and
following relations, the superscript ”bsg” and ”wr” stand for the blue
supergiant and Wolf-Rayet winds, respectively. The manner the wind velocity
changes over $\Delta t$ is controlled by,

$F=F(t,t_{\rm wr},\Delta t)=\Big{(}\frac{t-(t_{\rm wr}-\Delta t)}{\Delta
t}\Big{)}^{\beta/2},$ (25)

that is a function invoked when $t_{\rm wr}-\Delta t<t$ and $t\leq t_{\rm
wr}$, defined such that $F(t=t_{\rm wr}-\Delta t)=0$ and $F(t=t_{\rm wr})=1$.
For $\beta=2$ the transition is linear, while it is a power-law for $\beta\neq
2$. The changes in wind density are as follows,

$\displaystyle\dot{M}(r_{\rm in},t)=\Bigg{\\{}\begin{array}[]{lll}\dot{M}^{\rm
bsg}&\mathrm{if}&t\leq t_{\rm wr},\\\ \dot{M}^{\rm wr}&\mathrm{if}&t_{\rm
wr}<t,\\\ \end{array}$ (28)

with a sharp transition at $t_{\rm wr}$ (Brighenti & D’Ercole, 1997). These
relations implies that only the Wolf-Rayet wind is launched isotropically
while the previous winds are affected by the isotropy function $f(\theta)$.

Our MHD setup requires additional, similar prescriptions for the other
quantities necessary to fully fix both the Parker spiral and the passive
tracers tracking the mixing of material into the gas nebula. We introduce, in
the same fashion as done for the hydrodynamical case of Brighenti & D’Ercole
(1997), prescriptions for the toroidal and poloidal components of the stellar
surface magnetic field. They read,

$\displaystyle B_{\star}(r_{\rm
in},t)=\Bigg{\\{}\begin{array}[]{lll}B_{\star}^{\rm bsg}&\mathrm{if}&t\leq
t_{\rm wr}-\Delta t,\\\ B_{\star}^{\rm bsg}+\Delta
B_{\star}F&\mathrm{if}&t_{\rm wr}-\Delta t<tandt\leq t_{\rm wr},\\\
B_{\star}^{\rm wr}&\mathrm{if}&t_{\rm wr}<t,\\\ \end{array}$ (32)

with $\Delta B_{\star}=|B_{\star}^{\rm wr}-B_{\star}^{\rm bsg}|$.
Equivalently, the stellar radius evolves as,

$\displaystyle R_{\star}(r_{\rm
in},t)=\Bigg{\\{}\begin{array}[]{lll}R_{\star}^{\rm bsg}&\mathrm{if}&t\leq
t_{\rm wr}-\Delta t,\\\ \Delta R_{\star}(1-F)+R_{\star}^{\rm
wr}&\mathrm{if}&t_{\rm wr}-\Delta t<tandt\leq t_{\rm wr},\\\ R_{\star}^{\rm
wr}&\mathrm{if}&t_{\rm wr}<t,\\\ \end{array}$ (36)

with $\Delta R_{\star}=|R_{\star}^{\rm bsg}-R_{\star}^{\rm wr}|$. Finally, the
time-dependent evolution of the equatorial angular velocity reads,

$\displaystyle v_{\rm rot}(r_{\rm in},t)=\Bigg{\\{}\begin{array}[]{lll}v_{\rm
rot}^{\rm bsg}&\mathrm{if}&t\leq t_{\rm wr}-\Delta t,\\\ \Delta v_{\rm
rot}(1-F)+v_{\rm rot}^{\rm wr}&\mathrm{if}&t_{\rm wr}-\Delta t<tandt\leq
t_{\rm wr},\\\ v_{\rm rot}^{\rm wr}&\mathrm{if}&t_{\rm wr}<t,\\\ \end{array}$
(40)

with $\Delta v_{\rm rot}=|v_{\rm rot}^{\rm bsg}-v_{\rm rot}^{\rm wr}|$. Last,
the evolution of the passive scalar tracer at the stellar wind boundary obeys
the relation,

$\displaystyle Q(r_{\rm
in},t)=\Bigg{\\{}\begin{array}[]{lll}0&\mathrm{if}&t\leq t_{\rm wr}-\Delta
t,\\\ F&\mathrm{if}&t_{\rm wr}-\Delta t<tandt\leq t_{\rm wr},\\\
1&\mathrm{if}&t_{\rm wr}<t,\\\ \end{array}$ (44)

which permits to separate the cold, dusty, blue supergiant gas ($Q\leq 1/2$)
from the Wolf-Rayet material ($Q>1/2$).

Figure 1:  Number density field in our simulation models Run-Base (a) and Run-
Sym (b). The dashed green contour is the tangential discontinuity of the Wolf-
Rayet nebula, i.e. the location of the nebula made of equal proportion of blue
supergiant and Wolf-Rayet stellar winds, where the passive scalar $Q=0.5$. The
solid black and dashed-dotted red contours on the right-hand part of each
panels represent isovalues of the toroidal component of the velocity
$v_{\phi}$ and magnetic field $B_{\phi}$, respectively. The white cross marks
the position of the star.  Figure 2:  Cross-sections in the density field and
toroidal magnetic field of the evolving Wolf-Rayet nebula. The panels
correspond to different simulations, and on each panel the cross-section are
plotted for both directions along the polar (left) and equatorial (right)
directions of the domain.  Figure 3:  Same as Fig. 1 with our models Run-Base
and Run-$\Delta$t10kyr, exploring the effect of a phase transition timescale
$\Delta t\neq 0$.

###  2.4 Simulation models

The simulations are initialised with stellar wind properties corresponding to
a blue supergiant star. We let the system evolved to a few $10^{4}\,\rm yr$ up
to a time $t_{\rm wr}-\Delta t$. This time interval corresponds to the typical
timescale of a luminous blue variable phase (Smith, 2017) and it permits to
have potential magneto-hydrodynamical boundary effects provoked by the onset
of stellar rotation at $t=0$ to be transported out of the computational domain
under the effects of the constant wind inflow. The blue supergiant stellar
properties are as follows. We use $\dot{M}^{\rm bsg}=10^{-6}\,\rm
M_{\odot}\,\rm yr^{-1}$, which is of the order of magnitude of that of Sher 25
(Smartt et al., 2002), $v_{\rm w}^{\rm bsg}=500\,\rm km\,\rm s^{-1}$, that is
consistent with the values for HD 168625 ($350\pm 100\,\rm km\,\rm s^{-1}$,
Mahy et al., 2016) and for Sk$-$69∘279 ($800\pm 100\,\rm km\,\rm s^{-1}$,
Gvaramadze et al., 2018). The stellar radius is $R_{\star}^{\rm bsg}=20\,\rm
R_{\odot}$, which is that of $Sk-69^{\circ}279$ ($800\pm 100\,\rm km\,\rm
s^{-1}$, Gvaramadze et al., 2018), and $v_{\rm rot}^{\rm bsg}=60\,\rm km\,\rm
s^{-1}$ (Mahy et al., 2016). The surface magnetic field of the blue supergiant
ancestor star is taken to be $B_{\star}^{\rm bsg}=1\,\rm G$, which is of the
order of the magnetic fields observed in other cool stars (Kervella et al.,
2018; Vlemmings et al., 2002, 2005).

The toroidal component of the stellar magnetic field is scaled with the value
for the solar wind measured at $1\,\rm au$ (Herbst et al., 2020; Scherer et
al., 2020). The Wolf-Rayet stellar properties are as follows. We use
$\dot{M}^{\rm wr}=10^{-4.3}\,\rm M_{\odot}\,\rm yr^{-1}$, $v_{\rm w}^{\rm
wr}=1900\,\rm km\,\rm s^{-1}$, $R_{\star}^{\rm wr}=2.3\,\rm R_{\odot}$ and
$v_{\rm rot}^{\rm wr}=10\,\rm km\,\rm s^{-1}$, respectively. These quantities
are the values of the Wolf-Rayet star WR1 (WN4). Our stellar parameters are
summarised in Table 1. The stellar surface magnetic field is rather
unconstrained, and we decided to take $B_{\star}^{\rm wr}=100\,\rm G$.
Similarly, we assume that the toroidal component of the magnetic field scales
with the magnetic field in the solar wind (Herbst et al., 2020; Scherer et
al., 2020).

In this work, we use our toy model to explore the effects of the progenitor
stellar wind and phase transition properties on the shaping of the Wolf-Rayet
nebulae. The simulation models are summarized in Table 2.

Figure 4:  Same as Fig. 1 with our models Run-$\xi$20 and Run-$\xi$100,
exploring the effect of the flattening degree $p$ and density ratio $\xi$ of
the progenitor wind.

###  2.5 Radiative transfer calculations

As most circumstellar nebulae around evolved massive stars have been observed
in the (near)-infrared waveband, we consequently desire, in order to discuss
our models in the context of real data, to know how our simulated objects
would look like if observed at that particular waveband. We perform radiative
transfer calculations of our MHD simulation outputs against dust opacity using
the radmc-3d code (Dullemond, 2012). Dust density fields are constructed by
assuming a standard gas-to-dust ratio of $1/200$ and it is imported into
radmc-3d. The dust temperature is simulated by Monte-Carlo calculation using
the algorithm of Bjorkman & Wood (2001) and photon packages are ray-traced
from the stellar atmosphere to the stellar wind nebulae. The $24\,\mu\rm m$
near-infrared emissivity is estimated using the Silicates dust opacity of Laor
& Draine (1993) and it is projected onto the plane of the sky according to a
particular selected viewing angle. The irradiating central massive star is
assumed to be a spherical black body radiator of effective temperature $T_{\rm
eff}$ that is taken as,

$\displaystyle T_{\rm eff}(t)=\Bigg{\\{}\begin{array}[]{lll}T_{\rm eff}^{\rm
bsg}&\mathrm{if}&t\leq t_{\rm wr}-\Delta t,\\\ \Delta T_{\rm eff}(1-F)+T_{\rm
eff}^{\rm wr}&\mathrm{if}&t_{\rm wr}-\Delta t<tandt\leq t_{\rm wr},\\\ T_{\rm
eff}^{\rm wr}&\mathrm{if}&t_{\rm wr}<t,\\\ \end{array}$ (48)

where $\Delta T_{\rm eff}=|T_{\rm eff}^{\rm wr}-T_{\rm eff}^{\rm bsg}$—. The
stellar radius $R_{\star}(t)$ is determined from Eq. 36. The hot, Wolf-Rayet
material is distinguished from the cold blue supergiant gas using the passive
scalar tracer of Eq. 44. We produce emission maps for a distance to the
stellar source of $1\,\rm kpc$, that is typical for massive star-forming
regions.

##  3 Results

This section presents the results from the modelling of Wolf-Rayet gas nebulae
of blue supergiant ancestors. We compare the effects of the asymmetries and
evolutionary phases of various pre-Wolf-Rayet stellar winds on the geometry of
the gas nebulae.

###  3.1 Wolf-Rayet nebulae of asymmetric versus spherical blue supergiant
wind

In Fig. 1 we plot the number density field in our baseline simulation Run-Base
(left column of panels) and our model Run-Sym (right column of panels), which
assume asymmetric and isotropic blue supergiant progenitor stellar winds,
respectively. The dashed green contour is the tangential discontinuity, where
blue supergiant and Wolf-Rayet stellar winds are in equal proportions
($Q=0.5$). Solid black and dashed-dotted red contours trace $v_{\phi}$ and
$B_{\phi}$, respectively, and the white cross marks the position of the star.
The model Run-Sym obviously generates a spherical nebula (Fig. 1d,f) by fast-
wind-slow-wind collision, as modelled in the work of Meyer et al. (2020). In
model Run-Base the asymmetric blue wind is expelled as a torus of outflowing
material, see the studies of Raga et al. (2008) and Fang et al. (2017),
respectively, onto which stellar rotation and a Parker wind is superposed
(Fig. 1a). Both models have no phase transition, therefore, the Wolf-Rayet
stellar wind directly interacts with the luminous blue supergiant stellar
wind. The initial spherically-symmetric stellar wind of Run-Base is channelled
by its surroundings and adopts a bipolar morphology (Fig. 1c,e). One clearly
sees magnetic dipoles in both the old slow and new fast winds, whose field
lines are compressed in the shell region. Given the size of the nebula ($\leq
10\,\rm pc$), the asymmetric nebulae can develop inside of the region bordered
by the termination shock of the main-sequence stellar wind bubble (Meyer et
al., 2020). These rings should therefore form even in the context of the
supersonic bulk motion of the driving blue supergiant star (Fig. 1b,d,f). This
mechanism is therefore the MHD equivalent of the young aspherical Wolf-Rayet
nebulae of red supergiant stars of Brighenti & D’Ercole (1997), this time in
the context of a blue supergiant ancestor.

Fig. 2 shows cross-sections taken through the number density field of the
Wolf-Rayet nebulae. The figure plots the cross-sections, taken at intervals of
$0.1\,\rm kyr$ through several characteristic simulation snapshots of the blue
supergiant to Wolf-Rayet phase transition event. The left-hand part of each
panel corresponds to the polar direction, whereas the right-hand part of the
panels corresponds to the equatorial plane of the nebulae. The spherical
nebula in Run-Sym (Fig. 2b) is obviously similar along both the polar and
equatorial directions. The shell of swept-up blue wind grows under the effect
of the expelled Wolf-Rayet momentum and the contrast in density along both
direction decreases, i.e. the post-shock density of the Wolf-Rayet nebulae
diminishes with time and so does the corresponding emission (see Section 4.2).
At time $40.1\,\rm kyr$, the cross-sections of the asymmetric model Run-Base
exhibit clear differences between the polar and equatorial directions (Fig.
2a), as a result of the wind asymmetry.  When the Wolf-Rayet material begins
to blow, the shell of swept-up supergiant gas is denser in the equatorial
plane than in the direction perpendicular to the equator, see thick blue and
thin lines of Fig. 2a. This latitude-dependent density contrast is at the
origin of the asymmetric character of some observed Wolf-Rayet nebulae.

Figure 5:  Same as Fig. 1 with our models Run-$\beta$4 and Run-$\beta$2,
exploring the effects of the interpolation exponent $\beta$ during the phase
transition.

###  3.2 Effects of the phase transition timescale

Fig. 3 compares the density fields in our simulation models Run-Base (left
panels) and Run-$\Delta$t10kyr (right panels), which differs by the blue
supergiant to Wolf-Rayet phase transition interval. Such transition period is
non-existent in model Run-Base while it is $\Delta\,t=0.1\,\rm kyr$ in model
Run-$\Delta$t10kyr. Our phase transition timescale is in accordance with both
the typical duration of luminous blue phases $\sim 10^{4}\,\rm yr$ and with
the timescale of blue-to-red phase transition (Mackey et al., 2012; Groh et
al., 2014). The Wolf-Rayet nebula in Run-Base is generated by the simultaneous
changes of the wind velocity and density, engendering a peanut-like
morphology, which shocked layer accumulates compressed magnetic field lines,
as described in Section 3.1. In Run-$\Delta$t10kyr the stellar wind velocity
increases first, provoking a more elongated jet-like nebula, before that the
sudden changes in wind density fills it with denser, isotropic Wolf-Rayet
material (see Eq. 24-28). In the latter case, Rayleight-Taylor instabilities
form in the equatorial plane, as a consequence of the wind-wind collision at
work therein, see the discussion in Brighenti & D’Ercole (1997). This
comparison shows that the properties of the time windows during which the wind
properties evolve has a little influences on the overall appearance of Wolf-
Rayet nebulae. Our simplistic setup uses analytic prescriptions that mimic the
evolution of the central massive star, hence, the detailed shape of the nebula
might be slightly different with more realistic, time-dependent stellar wind
prescriptions, e.g. using pre-computed evolutionary tracks (Brott et al.,
2011).

In model Run-$\Delta$t10kyr, the accelerating stellar wind has already
interacted with the dense equatorial blue supergiant wind, resulting in a
bipolar shell (Fig. 3a,b) that is later filled by the Wolf-Rayet wind. The
wind-wind interaction is stronger for a reduced phase transition interval
$\Delta t$ (Fig. 3a,c), while the density of the shocked material in
Run-$\Delta$t10kyr is reduced, as a consequence of the diminished wind
momentum at the moment just before the sudden changes in mass-loss rate at
$t_{\mathrm{wr}}$. Nevertheless, in both models, the Wolf-Rayet wind
eventually melts with the shell of swept-up blue supergiant material, and,
comparable equatorial densities are reached (Fig. 3a) as it expands outwards
(Fig. 3c). Finally, once the hot gas filled the bipolar cavity, both nebulae
have similar overall dimensions except at the poles, where the magnetic field
lines are collimated (Fig. 3d). This further illustrates that the morphology
of older asymmetric Wolf-Rayet nebulae is governed by the surface
characteristics and stellar wind properties, such as terminal velocity and
mass-loss rate of the ancestor stars, as well as the manner in which the
evolutionary phase transition happens, but not by the ambient medium
properties.

###  3.3 Effects of the blue supergiant stellar wind asymmetries

Fig. 4 compares two simulation models, Run-$\xi$20 and Run-$\xi$100, both with
the same phase transition timescale $\Delta t=10\,\rm kyr$, but different
density flattening degrees towards the equator $p$ and density equator-to-pole
ratios $\xi$, namely, $p=1/4$ and $\xi=20$ (Run-$\xi$20) and $p=3/4$ and
$\xi=100$ (Run-$\xi$100), respectively. The different pre-Wolf-Rayet
morphologies of the circumstellar medium profoundly impact the development of
the shell generated by wind-wind interaction between blue supergiant to the
Wolf-Rayet winds. The denser the material in the equator, the more collimated
the Wolf-Rayet stellar wind (Fig. 4b), as a result of the thicker disc-like
blue supergiant distribution (Fig. 4a). On the other hand, the simulation with
a thinner equatorial plane produces a peanut-like nebula, recalling the
mechanisms for the shaping of the homunculus of $\eta$-Carina (Langer et al.,
1999; González et al., 2004; Hirai et al., 2021). Note the gradual opening of
the magnetic field lines as the wind speed accelerates (Fig. 3a-e,b-f). Again,
this is mostly shaped by the blue supergiant wind whose velocity gradually
increases during the phase transition up to reaching the fast values of the
Wolf-Rayet star (black dotted line in Fig. 4c,d). The discrepancy between
peanut-like and jet-like Wolf-Rayet nebulae persists at later times, when the
stellar phase transition keeps going (Fig. 4e,f).

As soon the Wolf-Rayet stellar wind interacts with the walls of the nebula,
the wind-wind interface becomes naturally much denser in the model with
$\xi=100$ than in the case with $\xi=20$ (Fig. 4c,d). This reflects in the
cross-sections taken through the density fields of the simulations (Fig.
4d,e). Indeed, the model with $\xi=20$ produces a larger shell in the
equatorial plane, whereas that with $\xi=100$ is more compact, thinner and
also much denser. The density jump at the expanding front is larger at time
$40.1\,\rm kyr$ in our model Run-$\xi$100, than in our simulation Run-p1/2 and
it stalls at $\approx 1\,\rm pc$ from the star, whereas the other nebula in
Run-p1/2 is more enlarged in the equatorial plane. The flattening degree of
the density towards the equator and the equator-to-pole density ratio turns to
be a major parameter in the shaping of Wolf-Rayet nebulae, as it controls
their final aspect ratio and therefore their overall morphology. Within the
explored parameter space, we therefore produce either peanut-like or jet-like
shapes.

###  3.4 Effects of phase transition properties

Our last series of comparison tests consist in changing the interpolation
between the blue supergiant and Wolf-Rayet stellar surface properties during
the phase transition. It is governed by the exponent parameter $\beta$ in Eq.
25, set to $\beta=2$ in our model Run-$\beta$2 (linear interpolation) and to
$\beta=4$ in our Run-$\beta$4 (squared polynomial interpolation). The
transition phase timescale is set to $\Delta t=5\,\rm kyr$. In Fig. 5a and
Fig. 5b, the pre-shaped circumstellar medium in which the Wolf-Rayet star
blows its material is similar, and how the blue supergiant stellar wind
gradually adopts a faster speed is the only effect responsible for differences
in the nebula’s morphologies. This keeps on going for another $0.1\,\rm kyr$,
until the tangential discontinuity between the two colliding winds reaches the
cavity. Nevertheless, when the full Wolf-Rayet wind blows through and fully
fills the cavity, morphological differences are rather small (Fig. 5e,f). The
field lines are opened under the effect of the faster equatorial rotation of
the Wolf-Rayet star. Instabilities developing in the equatorial plane, as in
model Run-$\Delta$t10kyr, appear in both runs, despite of the different,
longer phase transition. The apex of the polar lobes remains the only trace of
the different stellar evolution history of the final Wolf-Rayet stars.

###  3.5 Magnetic properties of the Wolf-Rayet nebulae

Figure 6:  Cross-sections in the density field (top panel), toroidal magnetic
field (middle panel) and Alfénic Mach number (bottom panel ) of the evolving
Wolf-Rayet nebula in our simulation Run-Base. On each panel the cross-section
are plotted for both the polar (left) and the equatorial direction (right) of
the computational domain.

In Fig. 6a we plot a time evolution series of cross-sections taken through the
number density field (top panel) in our simulation model assuming an
aspherical blue supergiant stellar wind (Run-Base). The selected times
correspond to the moment the Wolf-Rayet wind begins blowing into the blue
supergiant material. Fig. 6b (middle panel) displays cross-sections taken at
the same time instances and location, through the toroidal density field
$B_{\phi}$ (in $\mu\rm G$) and plotted in the logarithmic scale. One sees that
the magnetic field in the Wolf-Rayet nebula, governed by its toroidal
component is maximum in the region of the expanding shell. The compression
factor of the magnetic stellar wind field reflects that of the expanding shell
density (Fig. 6a,b), and, inversely, the toroidal field $B_{\phi}$ is weaker
along the polar direction. Fig. 6c (bottom panels) plots the Alfénic Mach
number in the circumstellar nebula. The gas is super-Alfénic everywhere, with
lower values $M_{\rm A}\sim 1.5$ in the region of the equatorial plane, in the
compressed shell, and it peaks at $M_{\rm A}\approx 4$ along the polar
direction.

Figure 7:  Correlation between the number density and the velocity of the gas
in the inner $5\,\rm pc$ of the computational domain, for several simulations
of the stellar wind nebulae generated during blue supergiant to Wolf-Rayet
phase transitions, with colours representing the strength of the toroidal of
the magnetic field $B_{\phi}$.

In Fig. 7 we display scatter plots for the distribution of gas in the
computational domain, representing the number density (in $\rm cm^{-3}$) of
the wind as a function of its velocity (in $\rm km\,\rm s^{-1}$), coloured by
the value of the toroidal magnetic field $B_{\phi}$ (in $\mu\rm G$). The
different panels of the figure correspond to several time instances of the
simulations, from the younger (top panels) to the older times (bottom panels),
and each column represent a simulation model. The constant, isotropically
freely-streaming stellar wind of the blue supergiant progenitor is represented
by an horizontal line of decreasing toroidal magnetic field from the high to
the lower densities (Fig. 7a,b). When the star evolves to the Wolf-Rayet
phase, the horizontal is split into two horizontal lines of different
velocities, a fast velocity, high density one (the Wolf-Rayet component), and
a slow, low-density one (the blue component), see Fig. 7c. The model Run-Base
with anisotropic stellar wind vertically spreads the horizontal line of Fig.
7a as a broadened band, whose height is a function of the degree of flattening
of the wind and of the equator-to-pole density ratio, see Fig. 7d,e. The phase
transition, with spherical Wolf-Rayet wind blown into an asymmetric blue
supergiant wind, shifts the high density part of the bar-like distribution to
a higher velocity, higher number density horizontal line (Fig. 7f). The region
between the two winds exhibits variations of the toroidal magnetic field from
the density compression during the wind-wind interaction, which results in
more scattering of the data in the velocity-density plane (Fig. 7f).

The model with both anisotropic pre-Wolf-Rayet stellar wind and long phase
transition timescale, Run-$\Delta$t$10\rm kyr$, exhibits a distribution in the
$\rho$-$v$ diagram both reporting the initial freely-expanding, high-density
stellar wind as a high-density horizontal bar, plus a low-density region
broadened in velocity, as a result of the blue supergiant asymmetries, already
at work at that time (Fig. 7g). At later times, the transition region greatly
affects the gas properties in the equatorial plane and this results in a
larger toroidal magnetic field (Fig. 7h,i). At even later times, the
development of instabilities at the wind-wind interface in the equatorial
plane strongly disperses the high-velocity gas in the $\rho$-$v$ diagram (Fig.
7i). Similarly, the simulation with Run-$\xi$100 exhibits a more important
scattering of the magnetised blue supergiant and expanding Wolf-Rayet stellar
winds. The velocity broadening is strong, as a result of the large equator-to-
pole density ratio $\xi$. The jet-like nebula has a dispersion of the material
therein that contains highly-magnetised, swept-up, slow diluted blue
supergiant stellar wind in the equatorial plane (Fig. 7k,l). The different
velocity distributions in peanut-like ($\xi=20$) and jet-like ($\xi=100$)
Wolf-Rayet nebulae clearly appear as a consequence of the aspect ratio of the
wind bubble, see Fig. 7i,l.

Figure 8:  Same as Fig. 7 for our models Run-$\beta$4 and Run-$\beta$2\.

Fig. 8 shows our models Run-$\beta$4 and Run-$\beta$2, respectively, which
differ by the manner the various quantities of the stellar wind are
interpolated during the phase transition. At time $40\,\rm kyr$ both nebulae
have already grown when the new-born Wolf-Rayet wind is blown into the
accelerated blue supergiant wind, as testifies the S-shape in the density-
velocity plane. The distribution of the gas therein presents both the
broadened low-velocity component from an equatorially-asymmetric blue
supergiant stellar wind and the horizontally-distributed component of a rapid
Wolf-Rayet wind. Both distributions with changing $\beta$ are qualitatively
similar (Fig. 8a,b), while differences at the transition between the two
components appear, as the wind-wind interaction region at the equatorial plane
is not the same. The model with non-linear interpolation ($\beta=4$) reveals
compressed gas with a stronger magnetisation $B_{\mathrm{\phi}}$. The phase
transition continues modifying the curved, S-shaped region of scattered dots
between the two components of the ancestor and evolved stellar winds. Such
curvature is more pronounced in the case with $\beta=4$ because the wind
momentum increases faster with the polynomial interpolation than in the linear
case with $\beta=2$, which results in an accumulation of dense, magnetised and
fast material in the equatorial plane (Fig. 8c,d). The high velocity component
turns into a thin horizontal line once the Wolf-Rayet wind is fully
established. The instabilities developing in those peanut-like Wolf-Rayet
nebulae are responsible for the scattering of points in the high-density,
high-velocity parts of the diagram.

Figure 9:  Emission maps of calculated from the Wolf-Rayet circumstellar
nebulae model Run-Base. The figure plots the $24\,\mu\rm m$ surface brightness
(in $\rm erg\,\rm cm^{2}\,\rm s^{-1}\,\rm arcsec^{-2}$) of our Wolf-Rayet
nebulae. Quantities are calculated excluding the non-dusty Wolf-Rayet material
plotted with the linear scale, under several viewing angles $\theta$. The
white cross marks the position of the star.  Figure 10:  Same as Fig. 9 for
our other models. Images assume an inclination angle of $\theta=45^{\circ}$.

##  4 Discussion

This section discusses the caveats of our method and compares our study to
results in precedent works. We also discuss the infrared emission properties
of the young Wolf-Rayet wind nebulae, and, finally, we compare our results
with observations.

###  4.1 Model limitations

The very first caveat of our simulations is their 2.5-dimensional nature, as
our method naturally imposes axisymmetry to the Wolf-Rayet nebulae (Heiligman,
1980; García-Segura et al., 1999). While this is acceptable in our pilot
study, performing full 3D MHD models better resolving the detailed formation
of, e.g. instabilities during the wind-wind interaction mechanisms (Brighenti
& D’Ercole, 1997; Meyer et al., 2020) should be considered in the future.
Secondly, we neglected several microphysical processes, such as radiation
transport in the close surroundings of the massive star. The high temperatures
of the Wolf-Rayet star ($\sim 10^{5}\,\rm K$) must affect the local
thermodynamical properties and the energy budget of the stellar wind, and, as
radiation pressure is at work therein, this should changes the gas dynamics
(Moens et al., 2021). Any chemical and/or non-ideal MHD mechanisms are also
totally neglected in our simulations. Their consideration, much beyond the
scope of our simplistic work, will be treated in follow-up studies. Last, our
study which continues the work of Brighenti & D’Ercole (1997), utilises
analytic prescriptions for the evolution of the stellar wind properties. It
has the advantage to directly use stellar properties derived from e.g.,
observations of blue supergiant stars (Smartt et al., 2002; Mahy et al., 2016;
Smith, 2017; Gvaramadze et al., 2018), however, at the cost of the
introduction of free parameters such as the timescale $\Delta t$ of the
corresponding phase transition. This can be circumvented by using self-
consistently calculated theoretical stellar evolution models (Brott et al.,
2011; Szécsi et al., 2020).

###  4.2 Infrared emission maps

In Fig. 9 we display $24\,\mu\rm m$ emission maps calculated on the basis of
our Wolf-Rayet nebula model Run-Base, using the RADMC-3D code. The surface
brightness is plotted in the linear scale in $\rm erg\,\rm cm^{2}\,\rm
s^{-1}\,\rm arcsec^{-2}$ for several viewing angles and the images are
calculated excluding the hot dust-free Wolf-Rayet gas, see Meyer et al.
(2020). On each panel the central white cross marks the position of the star.
The first column of panels displays the time-dependent edge-view evolution of
the close circumstellar medium of the massive stars. Before the evolution of
the blue supergiant star, the asymmetric wind generates a bright oblate disc-
like shape extending in the equatorial plane of the star (Fig. 9a). This model
does not have a long phase transition interval $\Delta t$, and the hot Wolf-
Rayet wind directly interacts with the dense equatorial colder gas. This
interaction region emits in the mid-infrared, which produces the bipolar
morphology (Fig. 9d). As the phase transition ends, the nebula further extends
and become fainter at the poles, while conserving its overall bipolar bubbly
shape. The second column of panels plots the same sequence of images
considered with an inclination angle $\theta=45^{\circ}$ with respect to the
plane of the sky. The situation is qualitatively the same as for
$\theta=0^{\circ}$, except that the bright infrared equatorial ring appears as
an ellipse (Fig. 9e,h). The last column of panels displays top-viewed nebulae,
which, as a consequence of the 2.5D nature of the simulations, look like rings
extending away from the evolving central star (Fig. 9f,i). The maximum
emission region is the interaction layer where the Wolf-Rayet wind reaches the
termination shock of the blue supergiant nebula (Fig. 9d,e,f).

Fig. 10 is as Fig. 9 for the other simulations models listed in Table. 1. The
left column of panels shows the nebula when the wind-wind collision begins,
while the right column of panels shows the circumstellar nebula around the
Wolf-Rayet stars at a later time. Both series of images are displayed with a
viewing angle of $\theta=45^{\circ}$ with respect to the plane of sky. The
model Run-$\Delta$t10kyr has a phase transition $\Delta t=10\,\rm kyr$
extending the time interval during which the Wolf-Rayet wind merges with the
previous blue supergiant material. This results in a denser region of wind-
wind interaction and a brighter bipolar nebula (Fig. 10a) materialised as an
ellipsoid in the projected equatorial plane (Fig. 10b). Our model Run-$\xi$20
has a different density flattening degree towards the equator. It emphasizes
the bipolar aspect of the nebula, which is less bubbly and more elongated
(Fig. 10c,d). Our model Run-$\xi$100 with a larger equator-to-pole density
ratio $\xi=100$ generates a tubular nebula as during the phase transition the
wind material is channeled into the direction perpendicular to the equator,
resulting in a less compressed and $24\,\mu\rm m$-fainter expanding shell. The
model Run-$\beta$4 with a non-linear interpolation of the stellar wind
quantities during the phase transition produces a quicker and more violent
wind-wind collision and a brighter interaction region (Fig. 10g,h). Our model
Run-$\beta$2 with a linear interpolation induces a more important compression
of the shell of stellar wind, and also a more pronounced bipolar morphology
(Fig. 10i,j).

###  4.3 Comparison with observations

####  4.3.1 Generalities

The geometry of Wolf-Rayet nebulae is a long-standing problem and still rather
unexplored mystery. All except one simulation performed to date predict their
appearance to be both spherical and affected by hydrodynamical instabilities
originating from the evolved Wolf-Rayet stellar wind, blown into the interior
of its own main-sequence stellar wind bubble and interacting with the shell
generated during its previous, e.g. red supergiant, evolutionary stage (Freyer
et al., 2003, 2006; Toalá & Arthur, 2011; Meyer et al., 2020). All these
models assume a static massive star undergoing a succession of isotropically-
launched stellar winds. As discussed above, a cool, dense stellar wind
launched along the equator of the star can make the Wolf-Rayet nebula strongly
anisotropic, shaping the outflow as a bipolar structure. The physical origin
of the equatorial stellar wind is so-far included as sub-grid models in the
few simulations run so far, see Brighenti & D’Ercole (1997). Our study is
therefore the first work continuing that on asymmetric Wolf-Rayet nebulae of
Brighenti & D’Ercole (1997).

The previous paper of our series explored the effects of the stellar bulk
motion on the formation of Wolf-Rayet stars. Indeed, the fastest runaway
massive stars of the Milky Way are strong-winded and bow-shock-free Wolf-Rayet
stars carrying circular rings, which challenged our understanding of the
formation of Wolf-Rayet nebulae (Meyer et al., 2020). A typical example is the
nebula M1-67 around the star WR124 (van der Sluys & Lamers, 2003). The
proposed simple solution for their formation scenario teaches us that the
stellar motion of the driving star does not influence the growth of Wolf-Rayet
nebulae, since the region of free-streaming between stellar surface and main-
sequence wind bubble’s termination shock is sufficiently extended. This is
largely the case for high-latitude objects, as the diluted medium therein
favors the production of huge and extended bow shocks, making room for wind-
wind shell collision to happen without being affected by the ISM material.
Importantly, the absence of observed ISM material in Wolf-Rayet nebulae
(Fernández-Martín et al., 2012; Esteban et al., 2016) indicates that most
Wolf-Rayet shells, including bipolar Wolf-Rayet nebulae, are probably produced
this way, and, this further implies that the large-scale wind-blown bubble can
be ignored when numerically modelling them.

However, the very complex shape of many observed Wolf-Rayet nebulae challenges
this picture (Toalá et al., 2015). When exactly bipolarity happens and which
factors are in control of this ? Which progenitors are responsible for which
kind of circumstellar nebulae ? The well-studied bipolar Wolf-Rayet nebulae
NGC 6888, of probable red supergiant ancestor, and RCW58, seemingly coming
from a blue supergiant, are the ideal cases to test numerical models like ours
against observations, see the discussion in Garcia-Segura et al. (1996) and
Garcia-Segura et al. (1996). Bipolarity of single blue supergiant origin is
the situation that we simulate in the present study, and it is in accordance
with the interpretation of Garcia-Segura et al. (1996) in the context of
RCW58. However, what can our result teach us regarding the other, red
supergiant-based scenario ?

####  4.3.2 What makes NGC 6888 asymmetric ?

The gaseous nebula NGC 6888, also called the veil nebula, the crescent nebula,
Caldwell 27, Sharpless 105, is an asymmetric wind bubble found around the
evolved massive Wolf-Rayet WN66(h)-type star WR 136. Its particular shape is
made of two bipolar lobes and of an equatorial ring of denser material. It has
long been studied in great detail (Parker, 1978; Treffers & Chu, 1982; Garcia-
Segura et al., 1996; Fernández-Martín et al., 2012; Mesa-Delgado et al., 2014;
Toalá et al., 2014; Esteban et al., 2016; Rubio et al., 2020). Early
observations of the surroundings of NGC 6888 revealed a variation of the level
ionization, traced by the optical line ratios [NII]/H$\alpha$ in the many
clumps along the surface of the shell, witnessing an efficient mixing of
materials preceding the formation of the Wolf-Rayet shell (Parker, 1978;
Fernández-Martín et al., 2012). The presence of dust and enriched elements
around Wolf-Rayet nebulae indicate that (i) they result from wind-wind
interaction of evolved material, (ii) one of the evolved winds is that of a
cold, dust-producing, supergiant phase of stellar evolution and (iii) the kind
of detected dust tells us that the wind from a main-sequence star does not
contribute to the shaping of the nebula. NGC 6888 has been the site of the
first observation of the aftermath of the CNO cycle of nuclear reactions into
the stellar core (Mesa-Delgado et al., 2014), and its progenitor star has been
constrained to be a single, red-supergiant-evolving star of $25-40\,\rm
M_{\odot}$ (Stock & Barlow, 2014; Mesa-Delgado et al., 2014). Hence, these
Wolf-Rayet stars must probably have undergone a red supergiant evolutionary
phase, and, such scenario has been modelled using 1D and 2D spherically-
symmetric hydrodynamical simulations in Garcia-Segura et al. (1996). However,
if it is clear that WR136 has been a cold star prior to the Wolf-Rayet phase
(Stock & Barlow, 2014), the origins of the progenitor’s asymmetries are still
not understood.

Our results show that a non-eruptive single blue supergiant star can not spin
sufficiently fast to produce an asymmetric pre-Wolf-Rayet circumstellar medium
(our Fig. 1). Since red supergiants and other late-type stars do not rotate
fast (Kervella et al., 2018; Vlemmings et al., 2002, 2005), their shells can
not be greatly affected by rotational and/or magnetic effects either.
Therefore, the question is how NGC 6888 can produce a bipolar Wolf-Rayet
nebula with an equatorial ring, obviously witnessing the existence of an
equatorial structure produced during a post-main-sequence evolutionary phase
prior to the Wolf-Rayet phase ? The most accepted picture, supported by CNO
yields measured in the nebula, made of a main-sequence, red supergiant and
Wolf-Rayet single star evolutionary channel, is therefore challenged by our
results, as explanations for the polar asymmetries of the pre-Wolf-Rayet wind
are still missing.

We speculate that a mechanism, similar to that revealed for asymptotic giant
stars, see particularly Decin et al. (2020) and Decin (2020), in which
multiplicity generates a latitude dependence of the stellar wind, should be at
work in the context of red supergiant stars as well, and in particular to the
central star(s) of NGC 6888 before it adopted a Wolf-Rayet spectral type.
Other possibility could be alternative evolution scenario involving an
(eruptive) luminous blue variable event, e.g. within a blue loop (Chita et
al., 2008; Mackey et al., 2012), however less in accordance with the dust size
measures of Rubio et al. (2020) shown to be typical for red supergiants. Their
constraint of the initial mass of WR136 to $\approx 50$-$60\,\rm M_{\odot}$
might reconcile both arguments, dust and asymmetries being generated at
different phases of the complex evolution of a (potentially multiple) stellar
system. Our allegation that the driving stellar object of NGC 6888 is not,
despite the findings of Stock & Barlow (2014); Rubio et al. (2020), a single
star, is supported by the detection of periodical modulation in line emission
from the Wolf-Rayet star WR 134, interpretable as an inhomogeneous outflow
(Morel et al., 1998) or as a trace of a so far undetected close companion
(Meyer et al., 2018).  The close companions to evolved massive stars
Koenigsberger & Schmutz (2020) and the circumstellar nebulae of extragalactic
Wolf-Rayet stars such as GR290 in the Galaxy M33 (Maryeva et al., 2018, 2020)
also support our arguments.  All this motivates further investigations of that
problem.

##  5 Conclusion

In this paper we study the possibility of forming asymmetric Wolf-Rayet
nebulae as a direct result of equatorial anisotropies in the blue supergiant
stellar wind preceding the Wolf-Rayet phase. We present a two-winds toy model
for the investigation of circumstellar medium generated from the evolutionary
phase transition of different latitude-dependent distributions. Our MHD
simulations are performed with a spherically-symmetric grid in the 2.5D
fashion, using the pluto code (Mignone et al., 2007, 2012; Vaidya et al.,
2018). We concentrate on the $\sim 10\,\rm pc$ surrounding the star, and treat
both stellar magnetic field and stellar rotation, using a combination of the
recipe for asymmetric winds of Raga et al. (2008) and magnetised astropsheres
(Parker, 1958). This study updates the pioneering work of Brighenti & D’Ercole
(1997) on the early aspherical evolution of young Wolf-Rayet nebulae from red
supergiant stars. Particularly, our toy model permits to change the rotation,
mass-loss, magnetic properties of the evolving stars, while simultaneously
exploring the phase transition interval between the two evolutionary stages.
It allows us to explore at reduced computational costs, the wind-wind
interaction of typical magnetised stellar outflows responsible for the shaping
of circumstellar nebulae. We also explore, by means of radiative transfer
calculations, the emission properties of those gas nebulae.

The density distribution in the pre-Wolf-Rayet surroundings is directly
responsible for the shaping of our Wolf-Rayet nebulae and the stellar magnetic
field seems to be dynamically unimportant in the process of producing the
bipolarities, in the sense that only evolving MHD flows in rotation produce
spherical nebulae, at least in the corner of the parameter space of blue
supergiant and subsequent Wolf-Rayet stars that we explore. An equatorially
compressed circumstellar distribution gives the Wolf-Rayet stellar wind a
peanut-like morphology, which is further elongated along the polar direction
as a jet-like shape if the polar-to-equatorial density distribution ratio
increases, or in case a substantially long phase transition period ($\sim\rm
kyr$) between the two winds. The topology of the magnetic field lines is
greatly modified when the Wolf-Rayet wind starts blowing into the previous
cold material. The rotating blue supergiant dipole $B_{\phi}$ sees its
magnetic field lines opened and channelled to the poles of the structure.
Interestingly, all our nebulae are shaped by wind-wind interaction taking
place in the free-streaming stellar wind of the ancestor star, i.e. inside of
the termination shock of its main-sequence stellar wind bubble (Weaver et al.,
1977). This implies that the nebulae will appear the same regardless of their
driving star’s bulk motion (Meyer et al., 2020) and that the nebulae are free
of ISM material. We propose to characterise the evolution of MHD asymmetric
Wolf-Rayet nebulae in the $\rho$-$v$-$B_{\phi}$ diagram, which we show to
concisely report the main feature of our MHD Wolf-Rayet nebulae as a
characteristic S-shape.

Radiative transfer calculations against dust opacity of the Wolf-Rayet nebulae
provided us with mid-infrared $24\,\mu\rm m$ synthetic emission maps, showing
that their projected emission reflects the anisotropies of the blue pre-Wolf-
Rayet wind. We demonstrate that projection effects play a non-negligible role
in the shaping of observed Wolf-Rayet nebulae. According to the observer’s
viewing angle, these transcient objects can appear as a oblate form, a bipolar
structure, an ellipse or a ring. Last, we discuss the potential origins of the
pre-Wolf-Rayet circumstellar anisotropies. Two main scenarios have been
identified so far, involving either a red or a blue supergiant star. However,
we show that solely the blue supergiant one can, within the single star
scenario, generate density distributions sufficiently asymmetric to channel
the Wolf-Rayet stellar wind and to produce Wolf-Rayet nebulae’s bipolarity.
This lead us to question the well-accepted past evolution of the famous
asymmetric bubble nebula NGC 6888 around the Wolf-Rayet WR 136, believed to be
the archetypal resulting nebulae of a single red supergiant-to-Wolf-Rayet-
evolving star. We propose that NGC 6888 might be driven by a multiple system,
or underwent additional evolutionary phases.

## Acknowledgements

The author thanks the referee A. Raga from the University of Mexico for his
numerous comments which greatly improved the quality of the paper.  The author
thanks L. Oskinova from the University of Potsdam for discussions regarding
evolved massive stars. The author acknowledges the North-German Supercomputing
Alliance (HLRN) for providing HPC resources that have contributed to the
research results reported in this paper.

## Data availability

This research made use of the pluto code developed at the University of Torino
by A. Mignone and collaborators (http://plutocode.ph.unito.it/) and of the
radmc-3d code developed by C. Dullemond and collaborators at the University of
Heidelberg (https://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/),
respectively. The figures have been produced using the Matplotlib plotting
library for the Python programming language (https://matplotlib.org/). The
data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Bestenlehner et al. (2014) Bestenlehner J. M., Gräfener G., Vink J. S. e. a., 2014, A&A, 570, A38
  * Bjorkman & Wood (2001) Bjorkman J. E., Wood K., 2001, ApJ, 554, 615
  * Blondin et al. (1996) Blondin J. M., Lundqvist P., Chevalier R. A., 1996, ApJ, 472, 257
  * Brighenti & D’Ercole (1997) Brighenti F., D’Ercole A., 1997, MNRAS, 285, 387
  * Brott et al. (2011) Brott I., de Mink S. E., Cantiello M., Langer N., de Koter A., Evans C. J., Hunter I., Trundle C., Vink J. S., 2011, A&A, 530, A115
  * Chevalier (1977) Chevalier R. A., 1977, ARA&A, 15, 175
  * Chevalier & Liang (1989) Chevalier R. A., Liang E. P., 1989, ApJ, 344, 332
  * Chevalier & Luo (1994) Chevalier R. A., Luo D., 1994, ApJ, 421, 225
  * Chita et al. (2008) Chita S. M., Langer N., van Marle A. J., García-Segura G., Heger A., 2008, A&A, 488, L37
  * Chu (1981) Chu Y. H., 1981, ApJ, 249, 195
  * Chu (1982a) Chu Y. H., 1982a, ApJ, 254, 578
  * Chu (1982b) Chu Y. H., 1982b, ApJ, 255, 79
  * Chu (1983) Chu Y. H., 1983, ApJ, 269, 202
  * Chu & Treffers (1981) Chu Y. H., Treffers R. R., 1981, ApJ, 250, 615
  * Chu et al. (1983) Chu Y. H., Treffers R. R., Kwitter K. B., 1983, ApJS, 53, 937
  * Ciotti & D’Ercole (1989) Ciotti L., D’Ercole A., 1989, A&A, 215, 347
  * Comerón & Kaper (1998) Comerón F., Kaper L., 1998, A&A, 338, 273
  * Crowther (2007) Crowther P. A., 2007, ARA&A, 45, 177
  * Decin (2020) Decin L., 2020, arXiv e-prints
  * Decin et al. (2020) Decin L., Montargès M., Richards A. M. S., Gottlieb C. A., Homan W., McDonald I. e. a., 2020, Science, 369, 1497
  * Dopita et al. (1994) Dopita M. A., Bell J. F., Chu Y. H., Lozinskaya T. A., 1994, ApJS, 93, 455
  * Dullemond (2012) Dullemond C. P., , 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library
  * Dwarkadas (2005) Dwarkadas V. V., 2005, ApJ, 630, 892
  * Dwarkadas (2007) Dwarkadas V. V., 2007, ApJ, 667, 226
  * Ekström et al. (2012) Ekström S., Georgy C., Eggenberger P., Meynet G., Mowlavi N., Wyttenbach A., Granada A., Decressin T., Hirschi R., Frischknecht U., Charbonnel C., Maeder A., 2012, A&A, 537, A146
  * Esteban et al. (2016) Esteban C., Mesa-Delgado A., Morisset C., García-Rojas J., 2016, MNRAS, 460, 4038
  * Faherty et al. (2014) Faherty J. K., Shara M. M., Zurek D., Kanarek G., Moffat A. F. J., 2014, AJ, 147, 115
  * Fang et al. (2017) Fang J., Yu H., Zhang L., 2017, MNRAS, 464, 940
  * Fernández-Martín et al. (2012) Fernández-Martín A., Martín-Gordón D., Vílchez J. M., Pérez Montero E., Riera A., Sánchez S. F., 2012, A&A, 541, A119
  * Flagey et al. (2014) Flagey N., Noriega-Crespo A., Petric A., Geballe T. R., 2014, AJ, 148, 34
  * Frank et al. (1995) Frank A., Balick B., Davidson K., 1995, ApJ, 441, L77
  * Freyer et al. (2003) Freyer T., Hensler G., Yorke H. W., 2003, ApJ, 594, 888
  * Freyer et al. (2006) Freyer T., Hensler G., Yorke H. W., 2006, ApJ, 638, 262
  * Gabler et al. (2021) Gabler M., Wongwathanarat A., Janka H.-T., 2021, MNRAS, 502, 3264
  * Garcia-Segura et al. (1996) Garcia-Segura G., Langer N., Mac Low M.-M., 1996, A&A, 316, 133
  * García-Segura et al. (1999) García-Segura G., Langer N., Różyczka M., Franco J., 1999, ApJ, 517, 767
  * Garcia-Segura et al. (1996) Garcia-Segura G., Mac Low M.-M., Langer N., 1996, A&A, 305, 229
  * García-Segura et al. (2018) García-Segura G., Ricker P. M., Taam R. E., 2018, ApJ, 860, 19
  * García-Segura et al. (2020) García-Segura G., Taam R. E., Ricker P. M., 2020, ApJ, 893, 150
  * González et al. (2004) González R. F., de Gouveia Dal Pino E. M., Raga A. C., Velázquez P. F., 2004, ApJ, 616, 976
  * Groh et al. (2014) Groh J. H., Meynet G., Ekström S., Georgy C., 2014, A&A, 564, A30
  * Guerrero & Chu (2008) Guerrero M. A., Chu Y.-H., 2008, ApJS, 177, 216
  * Gvaramadze et al. (2010) Gvaramadze V. V., Kniazev A. Y., Fabrika S., 2010, MNRAS, 405, 1047
  * Gvaramadze et al. (2018) Gvaramadze V. V., Kniazev A. Y., Maryeva O. V., Berdnikov L. N., 2018, MNRAS, 474, 1412
  * Hamann et al. (2006) Hamann W. R., Gräfener G., Liermann A., 2006, A&A, 457, 1015
  * Harten et al. (1983) Harten A., Lax P. D., van Leer B., 1983, SIAM Review, 25, 35
  * Heiligman (1980) Heiligman G. M., 1980, MNRAS, 191, 761
  * Herbst et al. (2020) Herbst K., Scherer K., Ferreira S. E. S., Baalmann L. R., Engelbrecht N. E., Fichtner H., Kleimann J., Strauss R. D. T., Moeketsi D. M., Mohamed S., 2020, ApJ, 897, L27
  * Hirai et al. (2021) Hirai R., Podsiadlowski P., Owocki S. P., Schneider F. R. N., Smith N., 2021, MNRAS, 503, 4276
  * Hung et al. (2021) Hung C. S., Ou P.-S., Chu Y.-H., Gruendl R. A., Li C.-J., 2021, ApJS, 252, 21
  * Kervella et al. (2018) Kervella P., Decin L., Richards A. M. S., Harper G. M., McDonald I., O’Gorman E., Montargès M., Homan W., Ohnaka K., 2018, A&A, 609, A67
  * Koenigsberger & Schmutz (2020) Koenigsberger G., Schmutz W., 2020, A&A, 639, A18
  * Langer (2012) Langer N., 2012, ARA&A, 50, 107
  * Langer et al. (1999) Langer N., García-Segura G., Mac Low M.-M., 1999, ApJ, 520, L49
  * Laor & Draine (1993) Laor A., Draine B. T., 1993, ApJ, 402, 441
  * Mackey et al. (2014) Mackey J., Mohamed S., Gvaramadze V. V., Kotak R., Langer N., Meyer D. M.-A., Moriya T. J., Neilson H. R., 2014, Nature, 512, 282
  * Mackey et al. (2012) Mackey J., Mohamed S., Neilson H. R., Langer N., Meyer D. M.-A., 2012, ApJ, 751, L10
  * Maeder (2009) Maeder A., 2009, Physics, Formation and Evolution of Rotating Stars, springer-verlag edn. Astronomy and Astrophysics Library
  * Mahy et al. (2016) Mahy L., Hutsemékers D., Royer P., Waelkens C., 2016, A&A, 594, A94
  * Marchenko et al. (2010) Marchenko S. V., Moffat A. F. J., Crowther P. A., 2010, ApJ, 724, L90
  * Marston et al. (1994) Marston A. P., Chu Y. H., Garcia-Segura G., 1994, ApJS, 93, 229
  * Marston et al. (1994) Marston A. P., Yocum D. R., Garcia-Segura G., Chu Y. H., 1994, ApJS, 95, 151
  * Maryeva et al. (2018) Maryeva O., Koenigsberger G., Egorov O., Rossi C., Polcaro V. F., Calabresi M., Viotti R. F., 2018, A&A, 617, A51
  * Maryeva et al. (2020) Maryeva O. V., Koenigsberger G., Karpov S. V., Lozinskaya T. A., Egorov O. V., Rossi C., Calabresi M., Viotti R. F., 2020, A&A, 635, A201
  * Mesa-Delgado et al. (2014) Mesa-Delgado A., Esteban C., García-Rojas J., Reyes-Pérez J., Morisset C., Bresolin F., 2014, ApJ, 785, 100
  * Meyer et al. (2018) Meyer D. M.-A., Kuiper R., Kley W., Johnston K. G., Vorobyov E., 2018, MNRAS, 473, 3615
  * Meyer et al. (2021) Meyer D. M. A., Mignone A., Petrov M., Scherer K., Velazquez P. F., Boumis P., 2021, arXiv e-prints, p. arXiv:2107.05513
  * Meyer et al. (2020) Meyer D. M. A., Oskinova L. M., Pohl M., Petrov M., 2020, MNRAS, 496, 3906
  * Meyer et al. (2020) Meyer D. M. A., Petrov M., Pohl M., 2020, MNRAS, 493, 3548
  * Meyer et al. (2021) Meyer D. M. A., Pohl M., Petrov M., Oskinova L., 2021, MNRAS, 502, 5340
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228
  * Mignone et al. (2012) Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G., 2012, ApJS, 198, 7
  * Moens et al. (2021) Moens N., Sundqvist J. O., El Mellah I., Poniatowski L., Teunissen J., Keppens R., 2021, arXiv e-prints, p. arXiv:2104.03968
  * Morel et al. (1998) Morel T., Marchenko S. V., Eenens P. R. J., Moffat A. F. J., Koenigsberger G., Antokhin I. I., Eversberg T., Tovmassian G. H., Hill G. M., Cardona O., St-Louis N., 1998, Ap&SS, 260, 173
  * Müller & Janka (2015) Müller B., Janka H. T., 2015, MNRAS, 448, 2141
  * Parker (1958) Parker E. N., 1958, ApJ, 128, 664
  * Parker (1978) Parker R. A. R., 1978, ApJ, 224, 873
  * Pogorelov & Matsuda (2000) Pogorelov N. V., Matsuda T., 2000, A&A, 354, 697
  * Pogorelov & Semenov (1997) Pogorelov N. V., Semenov A. Y., 1997, A&A, 321, 330
  * Pogorelov et al. (2004) Pogorelov N. V., Zank G. P., Ogino T., 2004, ApJ, 614, 1007
  * Powell (1997) Powell K. G., 1997, An Approximate Riemann Solver for Magnetohydrodynamics. Springer Berlin Heidelberg, Berlin, Heidelberg, pp 570–583
  * Prajapati et al. (2019) Prajapati P., Tej A., del Palacio S., Benaglia P., CH I.-C., Vig S., Mandal S., Kanti Ghosh S., 2019, ApJ, 884, L49
  * Raga et al. (2008) Raga A. C., Cantó J., De Colle F., Esquivel A., Kajdic P., Rodríguez-González A., Velázquez P. F., 2008, ApJ, 680, L45
  * Rozyczka & Franco (1996) Rozyczka M., Franco J., 1996, ApJ, 469, L127
  * Rubio et al. (2020) Rubio G., Toalá J. A., Jiménez-Hernández P., Ramos-Larios G., Guerrero M. A., Gómez-González V. M. A., Santamaría E., Quino-Mendoza J. A., 2020, MNRAS, 499, 415
  * Sander et al. (2012) Sander A., Hamann W. R., Todt H., 2012, A&A, 540, A144
  * Scherer et al. (2020) Scherer K., Baalmann L. R., Fichtner H., Kleimann J., Bomans D. J., Weis K., Ferreira S. E. S., Herbst K., 2020, MNRAS, 493, 4172
  * Smartt (2009) Smartt S. J., 2009, ARA&A, 47, 63
  * Smartt et al. (2002) Smartt S. J., Lennon D. J., Kudritzki R. P., Rosales F., Ryans R. S. I., Wright N., 2002, A&A, 391, 979
  * Smith (2017) Smith N., 2017, Philosophical Transactions of the Royal Society of London Series A, 375, 20160268
  * Stock & Barlow (2014) Stock D. J., Barlow M. J., 2014, MNRAS, 441, 3065
  * Szécsi et al. (2020) Szécsi D., Wünsch R., Agrawal P., Langer N., 2020, arXiv e-prints, p. arXiv:2004.08203
  * Toalá & Arthur (2011) Toalá J. A., Arthur S. J., 2011, ApJ, 737, 100
  * Toalá et al. (2012) Toalá J. A., Guerrero M. A., Chu Y. H., Gruendl R. A., Arthur S. J., Smith R. C., Snowden S. L., 2012, ApJ, 755, 77
  * Toalá et al. (2014) Toalá J. A., Guerrero M. A., Gruendl R. A., Chu Y. H., 2014, AJ, 147, 30
  * Toalá et al. (2015) Toalá J. A., Guerrero M. A., Ramos-Larios G., Guzmán V., 2015, A&A, 578, A66
  * Toalá et al. (2018) Toalá J. A., Oskinova L. M., Hamann W. R., Ignace R. e. a., 2018, ApJ, 869, L11
  * Treffers & Chu (1982) Treffers R. R., Chu Y. H., 1982, ApJ, 254, 569
  * Vaidya et al. (2018) Vaidya B., Mignone A., Bodo G., Rossi P., Massaglia S., 2018, ApJ, 865, 144
  * van der Sluys & Lamers (2003) van der Sluys M. V., Lamers H. J. G. L. M., 2003, A&A, 398, 181
  * van Marle et al. (2011) van Marle A. J., Keppens R., Meliani Z., 2011, Bulletin de la Societe Royale des Sciences de Liege, 80, 310
  * van Marle et al. (2006) van Marle A. J., Langer N., Achterberg A., García-Segura G., 2006, A&A, 460, 105
  * van Marle et al. (2005) van Marle A. J., Langer N., García-Segura G., 2005, A&A, 444, 837
  * van Marle et al. (2007) van Marle A. J., Langer N., García-Segura G., 2007, A&A, 469, 941
  * van Marle et al. (2011) van Marle A. J., Meliani Z., Keppens R., Decin L., 2011, ApJ, 734, L26
  * van Marle et al. (2015) van Marle A. J., Meliani Z., Marcowith A., 2015, A&A, 584, A49
  * Vink (2006) Vink J. S., 2006, in Lamers H. J. G. L. M., Langer N., Nugis T., Annuk K., eds, Stellar Evolution at Low Metallicity: Mass Loss, Explosions, Cosmology Vol. 353 of Astronomical Society of the Pacific Conference Series, Massive star feedback – from the first stars to the present. p. 113
  * Vlemmings et al. (2002) Vlemmings W. H. T., Diamond P. J., van Langevelde H. J., 2002, A&A, 394, 589
  * Vlemmings et al. (2005) Vlemmings W. H. T., van Langevelde H. J., Diamond P. J., 2005, A&A, 434, 1029
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Weiler & Sramek (1988) Weiler K. W., Sramek R. A., 1988, ARA&A, 26, 295
  * Weis et al. (1997) Weis K., Chu Y. H., Duschl W. J., Bomans D. J., 1997, A&A, 325, 1157
  * Woosley & Bloom (2006) Woosley S. E., Bloom J. S., 2006, ARA&A, 44, 507
  * Woosley et al. (2002) Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015
  * Wrigge et al. (2005) Wrigge M., Chu Y.-H., Magnier E. A., Wendker H. J., 2005, ApJ, 633, 248