# Mixing of materials in magnetised core-collapse supernova remnants

D. M.-A. Meyer1, M. Pohl1,2, M. Petrov3 and K. Egberts1
1 Universität Potsdam, Institut für Physik und Astronomie, Karl-Liebknecht-
Strasse 24/25, 14476 Potsdam, Germany
2 Deutsches Elektronen-Synchrotron (DESY), Platanenallee 6, 15738 Zeuthen,
Germany
3 Max Planck Computing and Data Facility (MPCDF), Gießenbachstrasse 2, D-85748
Garching, Germany
E-mail<EMAIL_ADDRESS>
(Last updated 2020 June 10; in original form 2013 September 5)

###### Abstract

Core-collapse supernova remnants are structures of the interstellar medium
(ISM) left behind the explosive death of most massive stars ($\lesssim 40\,\rm
M_{\odot}$). Since they result in the expansion of the supernova shock wave
into the gaseous environment shaped by the star’s wind history, their
morphology constitutes an insight into the past evolution of their progenitor
star. Particularly, fast-moving massive stars can produce asymmetric core-
collapse supernova remnants. We investigate the mixing of materials in core-
collapse supernova remnants generated by a moving massive $35\,\rm M_{\odot}$
star, in a magnetised ISM. Stellar rotation and the wind magnetic field are
time-dependently included into the models which follow the entire evolution of
the stellar surroundings from the zero-age main-sequence to $80\,\rm kyr$
after the supernova explosion. It is found that very little main-sequence
material is present in remnants from moving stars, that the Wolf-Rayet wind
mixes very efficiently within the $10\,\rm kyr$ after the explosion, while the
red supergiant material is still unmixed by $30\%$ within $50\,\rm kyr$ after
the supernova. black Our results indicate that the faster the stellar motion,
the more complex the internal organisation of the supernova remnant and the
more effective the mixing of ejecta therein. blackIn contrast, the mixing of
stellar wind material is only weakly affected by progenitor motion, if at all.

###### keywords:

MHD – stars: evolution – stars: massive – ISM: supernova remnants.

††pubyear: 2023††pagerange: Mixing of materials in magnetised core-collapse
supernova remnants–Mixing of materials in magnetised core-collapse supernova
remnants

##  1 Introduction

High-mass stars are defined as stellar objects of mass $M_{\star}$ heavier
than $8\,\rm M_{\odot}$ (Salpeter, 1955; Kroupa, 2001), and they represent
$\approx 1\%$ of all stars in the Milky Way. They form in dense and highly
opaque giant molecular clouds, as a result of the local free-fall
gravitational collapse of dense cold gas (Shu, 1992). They follow a formation
scenario that is a scaled-up version of that of low-mass stars (Vorobyov &
Basu, 2006; Vorobyov, 2009, 2010), and which is well-described by the so-
called burst mode of accretion in star formation (Meyer et al., 2017; Meyer et
al., 2021a). When the infalling molecular material lands onto a centrifugally-
balanced accretion disc, it undergoes gravitational instability and gaseous
clumps can form and inward-migrate onto stellar surface, adding mass to the
young massive star (Meyer et al., 2019). Once massive protostars are hot
enough, their ionizing feedback dissipate their surrounding molecular
environment and they release strong supersonic winds from the stellar surface
(Vink et al., 2000) to interact directly with the interstellar medium (ISM),
giving birth to stellar wind bubbles (Weaver et al., 1977).

As a function of multiplicity (Sana et al., 2012), metallicity (Brott et al.,
2011a, b) and rotation properties (Maeder & Meynet, 2000; Szécsi et al.,
2022), massive stars evolve from a long initial main-sequence phase by
undergoing several evolutionary phases, characterised by abrupt surface
properties changes, alternating between hot, fast-winded with cold, slow-
winded phases (Langer, 2012). Although exotic evolutionary phases exist
(Smith, 2014; Groh et al., 2014), a high-mass star typically goes through a
red supergiant phase (Levesque, 2010; Davies, 2017), followed, in the case
$M_{\star}\geq 30\,\rm M_{\odot}$, by a Wolf-Rayet phase (Underhill, 1968;
Hamann, 1992; Crowther, 2007), before exploding as a supernova (Woosley &
Weaver, 1986; Weiler & Sramek, 1988; Smartt, 2009). The circumstellar medium
of massive stars is the region of the ISM directly surrounding high-mass
stellar objects and it adopts morphologies and emission properties reflecting
both the stellar feedback and the ambient medium properties.

The surroundings of young massive stellar objects are massive accretion discs
of radius $\sim 1000\,\rm au$ (Johnston et al., 2015; Ilee et al., 2016;
Ahmadi et al., 2018) that are evaporated and blown by the radiation field and
the wind of zero-age main-sequence high-mass stars. Their environments adopt
the theoretical shape of the so-called stellar wind bubble described in Weaver
et al. (1977), which is distorted into a stellar wind bow shock (Gull & Sofia,
1979; Wilkin, 1996) in the case of the $30\%$ of runaway massive stars which
move supersonically through the ISM (Blaauw, 1961; Gies, 1987; Gvaramadze et
al., 2012; Moffat et al., 1998). Such bow shock can be observed throughout the
whole electromagnetic spectrum, from the ultraviolet (Kaper et al., 1997) and
the infrared waveband (Peri et al., 2012; Peri et al., 2015; Kobulnicky et
al., 2016, 2017) to the high-energies (Benaglia et al., 2010; del Valle et
al., 2013; De Becker et al., 2017). When the star evolves to the red
supergiant phase, the wind interacts either with the interior of the main-
sequence wind bubble (Garcia-Segura et al., 1996b, a; Freyer et al., 2006), or
directly with the ISM if the star moves fast enough (Brighenti & D’Ercole,
1995b; Chita et al., 2008; Herbst et al., 2022). The circumstellar medium of
Wolf-Rayet stars takes the form of circular rings in the context of fast-
moving, high-latitude objects, bipolar nebulae in the case of rotating massive
stars or more complex morphologies (Brighenti & D’Ercole, 1995a; Meyer et al.,
2020b; Meyer, 2021). The surroundings of evolved massive stars is therefore
the fingerprint of its past life onto the ISM, and, it constitutes the
location of the supernova explosion. The supernova remnants left behind the
death of high-mass stars are shaped by the interaction of the supernova shock
wave with the circumstellar medium (Weiler & Sramek, 1988). Unfolding them
permit to reconstruct the past evolution and to constrain the properties of
its progenitor star.

Supernova remnants of massive progenitor stars, or core-collapse supernova
remnants, are characterised by the presence of circumstellar material enriched
in the heavy elements synthetised by the progenitor released as stellar winds.
They reveal the presence of elements such as Fe (Uchida et al., 2009), that
have been produced in the core of the progenitor and by the possible presence
of a neutron star left inside of the supernova remnant (Zavlin et al., 1998;
Gvaramadze, 2006), ultimately leading to the formation of pulsar wind nebulae
(Olmi & Bucciantini, 2023, 2019). The supernova explosion is closely connected
to non-thermal physical mechanisms such as the anisotropic liberation of
neutrinos (Kifonidis et al., 2006; Müller et al., 2012; Müller & Janka, 2015)
and the production of gravitational waves (Jardine et al., 2022). Models from
the early instants (Aloy & Obergaulinger, 2021) up to the first year of core-
collapse event has been subject to numerical modelling in Gabler et al.
(2021). The interaction of the liberated supernova shock wave with the last
wind of the progenitor generates characteristic lightcurve functions of the
wind properties, which have been observed and modelled (Chevalier, 1982;
Chevalier & Liang, 1989; Truelove & McKee, 1999). The propagating shock wave
interacts with the circumstellar material before further expanding into the
interstellar medium (Tutone et al., 2020), leading to the production of non-
thermal emission (Orlando et al., 2019; Das et al., 2022). A well-documented
example of such mechanism is that of Cassiopeia A (Finn et al., 2016; Holland-
Ashford et al., 2020).

Once the shock wave further expands, its propagation is constrained by the
distribution of circumstellar material shaped by the supernova progenitor, see
e.g. the case of the historical core-collapse supernova remnant Kes 73
(Katsuda et al., 2018). The reflection towards the center of the explosion of
the shock wave against its own circumstellar matter produces a central region
emitting thermal X-rays that is surrounded by a radio shell. This class of
supernova remnants is called mixed-morphology supernova remnant (Yusef-Zadeh
et al., 2003; Arias et al., 2019; Chiotellis et al., 2019; Derlopa et al.,
2019). Most evolved supernova remnants exhibit morphologies deviating from
sphericity. These asymmetries might have many origins, such as the intrinsic
asymmetry of the explosion (Orlando et al., 2020), but also the aspherical
distribution of circumstellar material like an asymmetric (Blondin et al.,
1996) wind or a bow shocks generated by runaway massive stars. The organised
magnetic field of the ISM imposes stellar wind bubbles of massive stars an
elongated morphology (van Marle et al., 2015), which, in its turn, imposes
rectangular shapes to the supernova remnants (Meyer et al., 2022). On the long
term, they can reflect in the subsequent pulsar wind nebulae of plerionic
core-collapse supernova remnants (Meyer & Meliani, 2022; Bandiera et al.,
2023).

In this study, we explore by means of 2.5D MHD simulations the morphology and
emission properties of core-collapse supernova remnants generated by runaway
massive stars, including the rotation and magnetisation of the stellar winds.
Furthermore, our approach made of high-resolution simulations using a uniform
cylindrical grid (Meyer et al., 2015; Meyer et al., 2020a; Meyer et al.,
2021b) is updated with the inclusion of a system of passive scalar tracers
which permit us to time-dependently follow the materials successfully expelled
from the stellar surface throughout the many evolutionary phases undergone by
the high-mass star. The circumstellar medium of the supernova progenitor is
first simulated from the zero-age main-sequence of the star until the
supernova time using methods similar to that of Comerón & Kaper (1998); van
Marle et al. (2005, 2007, 2011); Mackey et al. (2012); van Marle et al.
(2014); Meyer et al. (2014, 2016, 2017, 2020b); Green et al. (2019). Then, the
early supernova-wind interaction is calculated in a one-dimensional fashion
(Chevalier, 1982; Truelove & McKee, 1999; Dwarkadas, 2005, 2007; Whalen et
al., 2008; van Marle et al., 2010; Telezhinsky et al., 2012a, b, 2013;
Chiotellis et al., 2013; Broersen et al., 2014) and mapped onto the
circumstellar medium in order to model the evolution of the supernova remnant
up to a time $\sim 100\,\rm kyr$ after the explosion, see also Tenorio-Tagle
et al. (1990, 1991); Velázquez et al. (2006); van Veelen et al. (2009); Vigh
et al. (2011); Chiotellis et al. (2012); Toledo-Roy et al. (2014); van Marle
et al. (2015); Fang et al. (2017); Chiotellis et al. (2020, 2021). By
performing a resolution study, we discuss the internal organisation of the
supernova remnants, such as the mixing of material at work therein, and, by
radiative transfer calculation, we predict their non-thermal appearance by
synchrotron emission.

This paper is organised as follows. First, the numerical methods used to model
magnetised core-collapse supernova remnants of runaway, rotating, magnetised
Wolf-Rayet-evolving massive stars are presented in Section 2. We detail the
internal properties of the supernova remnants of a runaway $35\,\rm M_{\odot}$
progenitor in Section 3. We further discuss our results in Section 4 and
finally draw our conclusions in Section 5.

##  2 Method

In this section we present the methods used to simulate the core-collapse
supernova remnants of a runaway rotating magnetised $35\,\rm M_{\odot}$
massive star.

###  2.1 Governing equations

Our modelling of core-collapse supernova remnants is conducted within the
frame of non-ideal magneto-hydrodynamics (Meyer et al., 2020b). The time
evolution of a plasma volume element of the system is described by the mass
continuity equation,

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

the equation for the conservation of linear momentum,

$\frac{\partial\mn@boldsymbol{m}}{\partial
t}+\mn@boldsymbol{\nabla}\cdot\Big{(}\mn@boldsymbol{m}\textcolor{black}{\otimes}\mn@boldsymbol{v}\textcolor{black}{-}\mn@boldsymbol{B}\textcolor{black}{\otimes}\mn@boldsymbol{B}+\mn@boldsymbol{\hat{I}}p_{\rm
t}\Big{)}=\mn@boldsymbol{0},$ (2)

the total energy conservation equation,

$\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 the equation for the evolution of the magnetic field,

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

where $\mn@boldsymbol{m}=\rho\mn@boldsymbol{v}$ is the vector momentum, $\rho$
the gas density, $v$ the flow velocity, $p_{t}$ the total pressure and
$\vec{B}$ the magnetic field vector, 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},$
(5)

and the sound speed closing the system of equations Eq. 1-4 is defined as

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

with $\gamma=5/3$ the adiabatic index.

Optically-thin radiative cooling and heating is included into the equations
using the laws presented in Meyer et al. (2014). They are represented by the
function,

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

where

$n_{\mathrm{H}}=\frac{\rho}{\mu(1+\chi_{\rm He,Z})m_{\mathrm{H}}},$ (8)

is the hydrogen number density, with $\mu$ the adiabatic index, $\chi_{\rm
He,Z}$ the solar helium and metal abundance, and where the quantity

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

is the temperature of the gas, respectively. The functions $\Gamma$ and
$\Lambda$ in Eq. 7 stand for the heating and cooling rates for photoionised
gas per unit time. They account for cooling from H, He and metals Z (Wiersma
et al., 2009) within the solar abundance (Asplund et al., 2009), which largely
dominates the cooling at high T. It additionally includes components for
hydrogen recombination (Hummer, 1994), collisionally excited forbidden line
emission for the metals of atomic number $>2$ such as O, C and N (Raga et al.,
1997). The heating rate represents the photon emitted by the photosphere of
the massive star, which ionize the recombining hydrogenoic ions and provoke
the release of energetic electrons (Osterbrock & Bochkarev, 1989).

The equations are integrated using the pluto
code111http://plutocode.ph.unito.it/ (Mignone et al., 2007, 2012; Vaidya et
al., 2018) using a Godunov-type solver composed of the shock-capturing HLL
Rieman solver (Harten et al., 1983) and the eight-wave finite-volume algorithm
(Powell, 1997) for magnetised flows, which makes sure that the magnetic field
remains divergence-free,

$\vec{\nabla}\cdot\vec{B}=0,$ (10)

throughout the whole computational domain and over the entire simulation time,
from the onset of the stellar wind to the middle-age evolution phase of the
supernova remnant. It has been utilised in the precedent studies of the same
authors devoted to the circumstellar medium of evolving and dying massive
stars (Meyer et al., 2021c). Furthermore, the numerical scheme uses a third-
order Runge-Kutta time integrator, the minmod flux limiter and the WENO3
interpolation scheme for the reconstruction of the primitive variables of the
system. The time-marching algorithm is controlled by the Courant-Friedrich-
Levy rule that is initialised to $C_{\rm cfl}=0.1$ (Meyer, 2021).

Figure 1:  Surface and stellar wind properties of a $35\,\rm M_{\odot}$
supernova progenitor. The left-hand figures display the time evolution of the
stellar mass $M_{\star}$ (in $\rm M_{\odot}$, panel a), the stellar radius
$R_{\star}$ (in $\rm R_{\odot}$, panel b), the effective temperature $T_{\rm
eff}$ (in $\rm K$, panel c). The right-hand figures show the mass-loss rate
$\dot{M}$ (in $\rm M_{\odot}\,\rm yr^{-1}$, panel d), the terminal wind
velocity $v_{\rm w}$ (in $\rm km\,\rm s^{-1}$, panel e) and the equatorial
rotation profile of the star $v_{\rm rot}$ (in $\rm km\,\rm s^{-1}$, panel f),
respectively.

###  2.2 Modelling the circumstellar medium

####  2.2.1 Stellar wind density

The circumstellar medium of massive stars is modelled using a 2.5D cylindrical
coordinate system (R,z) and a $[0;100]\times[-50;50]\,\rm pc^{2}$
computational domain of $8000\times 8000$ grid zones. The stellar wind of the
massive star is injected at the surface of a sphere of radius $r_{\rm in}20$
cells centered onto the origin $O(0;0;0)$. The density is taken to be the
profile

$\rho_{w}(r=r_{\rm in},t)=\frac{\dot{M}(t)}{4\pi r^{2}v_{\rm w}(t)},$ (11)

where $\dot{M}$ is the mass-loss rate of the massive star, $v_{\rm w}$ the
terminal wind velocity of the star and $r$ the radial direction. As in the
previous papers of this series, the stellar wind nebula developing around the
star is calculated in its reference frame. The stellar motion is mimiced by
imposing a gas flow $v=-v_{\star}$ at the outer boundary of the computational
domain at $z=50\,\rm pc$, with $v_{\star}=20\,\rm pc$ the stellar speed
through the ISM, see also Comerón & Kaper (1998); Green et al. (2019); Meyer
et al. (2020b). The uniform ambient medium is in the warm phase of the ISM,
with a temperature of $T_{\rm ISM}\approx 8000\,\rm K$ and a magnetisation of
$B_{\rm ISM}=7\,\mu\,\rm G$ taken to be parallel to the $Oz$ symmetry axis of
the computational domain (Meyer et al., 2017).

####  2.2.2 Stellar wind velocity

The stellar wind properties of the massive star that we consider are shown in
Fig. 1. The figure plots the time evolution of various stellar global and
surface quantities such as the stellar mass $M_{\star}$ (Fig. 1a), the stellar
radius $R_{\star}$ (Fig. 1b), the effective temperature $T_{\rm eff}$ (Fig.
1c), as well as the mass-loss rate (Fig. 1d), that are directly taken from the
stellar evolution model of the geneva
library222https://www.unige.ch/sciences/astro/evolution/en/database/syclist/
(Ekström et al., 2012) with an initial mass of $35\,\rm M_{\odot}$ and a
rotation velocity of

$\frac{\Omega_{\star}(t=0)}{\Omega_{\rm K}}=0.1,$ (12)

and with

$\Omega_{\star}(t)=\frac{v_{\rm rot}(t)}{R_{\star}(t)}$ (13)

the angular frequency of the star and $\Omega_{\rm K}$ the equatorial
Keplerian velocity. The star undergoes three distinct evolutionary phase, i.e.
a hot main-sequence with a mass-loss rate of $\dot{M}\approx 10^{-6.3}\,\rm
M_{\odot}\,\rm yr^{-1}$ lasting at time $\approx 5.1\,\rm Myr$, a short,
$\approx 0.25\,\rm Myr$-long cold red supergiant phase with a dense wind of
$\dot{M}\approx 10^{-4.0}\,\rm M_{\odot}\,\rm yr^{-1}$ and a final Wolf-Rayet
phase with $\dot{M}\approx 10^{-5.0}\,\rm M_{\odot}\,\rm yr^{-1}$. The star
loses most of its mass during the two last phases (Fig. 1a) and the supernova
time is at $\approx 0.63\,\rm Myr$. The wind terminal velocity is calculated
as

$\frac{v_{\rm w}(t)}{v_{\rm esc}(t)}=\sqrt{\beta(T)},$ (14)

or equivalently,

$v_{\rm w}(t)=\sqrt{\beta(T)\frac{2GM_{\star}(t)}{R_{\star}(t)}},$ (15)

with the escape speed of the massive star $v_{\rm esc}(t)$, its stellar radius
$R_{\star}(t)$ (Fig. 1b), and

$\beta_{\rm w}(T)=\begin{cases}1.0&\text{if }T_{\rm eff}\leq 10000\,\rm K,\\\
1.4&\text{if }T_{\rm eff}\leq 21000\,\rm K,\\\ 2.65&\text{if }T_{\rm
eff}>21000\,\rm K,\\\ \end{cases}$ (16)

a corrective factor that is function of the effective temperature (Eldridge et
al., 2006), see Fig. 1c. The evolution of $v_{\rm w}$ decreases as a function
of time, as the main-sequence wind is supersonic ($\approx 1800\,\rm km\,\rm
s^{-1}$), the red supergiant wind velocity is very slow ($\approx 30\,\rm
km\,\rm s^{-1}$) and the final Wolf-Rayet wind velocity increases up to
($\approx 700\,\rm km\,\rm s^{-1}$), see Fig. 1e.

Figure 2:  blackblackTemporal evolution of the surface magnetic field
strength, $B_{\star}$, (panel a), the magnetic field strength in the stellar
wind, $B_{\rm w}$, at $r=0.03\,\rm au$ from the star (panel b), and the values
of the passive scalar tracers (panel c) for the different wind phases of the
stellar evolution.

####  2.2.3 Stellar rotation and magnetic field

The stellar magnetic field is imposed at the inner stellar wind boundary as a
Parker spiral (Parker, 1958; Weber & Davis, 1967; Pogorelov & Semenov, 1997;
Pogorelov & Matsuda, 2000) that is superimposed to the stellar wind (Chevalier
& Luo, 1994; Rozyczka & Franco, 1996; García-Segura et al., 2018, 2020). It is
composed of both a radial

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

and of a toroidal component

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

respectively, with

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

the latitude-corrected azimuthal component of the surface stellar velocity,
where the factor $\sin(\theta)$ account for the latitude-dependence of the
surface rotation (Fig. 1f), The third poloidal component of the magnetic field
is set to $B_{\theta}=0\,\rm G$. As the magnetic field is transported into the
stellar wind of the massive star, we scale the decrease of the strength of the
surface magnetic field to that of the Sun, via the relation

$\frac{B_{\star}(t)}{B_{\odot}}=\frac{B_{\rm w}(t,r=1\,\rm au)}{B_{\rm
w,\odot}(r=1\,\rm au)},$ (20)

where $B_{\rm w}(t,r=1\,\rm au)$ is the star’s magnetic field at a distance of
$r=1\,\rm au$ from the stellar center, $B_{\odot}$ the Sun’s surface magnetic
field and $B_{\rm w,\odot}(r=1\,\rm au)$ the observed strength of the magnetic
field in the solar wind, respectively. This recipe has been used as an
addition to boundary conditions for magnetised stellar winds in the context of
massive stars (Scherer et al., 2020; Herbst et al., 2020; Baalmann et al.,
2020, 2021; Meyer et al., 2021c).

The time-dependence of the stellar surface magnetic field $B_{\star}(t)$ is
estimated by assuming a typical field strength for each individual
evolutionary phase of the massive star. The main-sequence field is taken to
$B_{\star}=500\,\rm G$ (Fossati et al., 2015; Castro et al., 2015; Przybilla
et al., 2016; Castro et al., 2017). Regarding the red supergiant phase, we use
the value derived for Betelgeuse, i.e. $B_{\star}=0.2\,\rm G$ (Vlemmings et
al., 2002, 2005; Kervella et al., 2018). Finally, we assume the ad-hoc value
of $B_{\star}=100\,\rm G$ for the final Wolf-Rayet phase (Meyer, 2021). These
assumed values are reported in Fig. 2a and permit to build a magnetic field
history at any radius of the stellar wind boundary $r=r_{\rm in}$, see Fig.
2b.

####  2.2.4 Mixing of materials

Last, a system of passive scalar tracers permit us to time-dependently follow
the evolution and mixing of stellar wind materials into the circumstellar
medium and the supernova remnant. Since the $35\,\rm M_{\odot}$ of the geneva
library (Ekström et al., 2012) undergoes three typical evolutionary phases,
i.e. an initial long main-sequence phase, a supergiant phase and final Wolf-
Rayet phase, we introduce the corresponding tracers $Q_{\rm MS}$, $Q_{\rm
RSG}$ and $Q_{\rm WR}$, obeying the continuity equations

$\frac{\partial(\rho Q_{\rm MS})}{\partial
t}+\mn@boldsymbol{\nabla}\cdot(\mn@boldsymbol{v}\rho Q_{\rm MS})=0,$ (21)
$\frac{\partial(\rho Q_{\rm RSG})}{\partial
t}+\mn@boldsymbol{\nabla}\cdot(\mn@boldsymbol{v}\rho Q_{\rm RSG})=0,$ (22)

and

$\frac{\partial(\rho Q_{\rm WR})}{\partial
t}+\mn@boldsymbol{\nabla}\cdot(\mn@boldsymbol{v}\rho Q_{\rm WR})=0,$ (23)

respectively. The initial value of the passive scalar is set to $Q_{\rm i}=1$
during the time the star spends in the evolutionary phase $i=\\{\rm MS,\rm
RSG,\rm WR\\}$ and $Q_{\rm i}=0$ otherwise. This implies that $Q_{\rm
MS}+Q_{\rm RSG}+Q_{\rm WR}=1$ at every time of the simulation at the stellar
wind boundary $r=r_{\rm in}$.

###  2.3 Modelling the ejecta-wind interaction

The interaction between the supernova ejecta and the isotropic freely-
expanding stellar wind of the progenitor is calculated using a 1D spherically-
symmetric computational domain [$O$;$r_{\rm out}$]. The supernova ejecta
initially occupy the $0\leq r\leq r_{\rm max}<r_{\rm out}$ region of the
domain. The corresponding density field reads,

$\rho(r)=\begin{cases}\rho_{\rm core}(r)&\text{if $r\leq r_{\rm core}$},\\\
\rho_{\rm max}(r)&\text{if $r_{\rm core}<r<r_{\rm max}$},\\\ \end{cases}$ (24)

where $0<r_{\rm core}$ is filled by the constant high-density plateau

$\rho_{\rm core}(r)=\frac{1}{4\pi n}\frac{(10E_{\rm
ej}^{n-5})^{-3/2}}{(3M_{\rm ej}^{n-3})^{-5/2}}\frac{1}{t_{\rm max}^{3}},$ (25)

while the ejecta velocity display an homologous expansion

$v(r)=\frac{r}{t},$ (26)

which is reaching

$v_{\rm core}=\bigg{(}\frac{10(n-5)E_{\rm ej}}{3(n-3)M_{\rm
ej}}\bigg{)}^{1/2}$ (27)

at a distance $r=r_{\rm core}$ from the center of the explosion, respectively.
The region $r_{\rm core}\leq r\leq r_{\rm max}$ contains the envelope of
ejecta made of material of decreasing density,

$\rho_{\rm max}(r)=\frac{1}{4\pi n}\frac{\left(10E_{\rm
ej}^{n-5}\right)^{(n-3)/2}}{\left(3M_{\rm
ej}^{n-3}\right)^{(n-5)/2}}\frac{1}{t_{\rm max}^{3}}\bigg{(}\frac{r}{t_{\rm
max}}\bigg{)}^{-n},$ (28)

with $n=11$, see Chevalier (1982); Truelove & McKee (1999). At the edge of the
supernova shock wave the gas flow speed is $v_{\rm max}=30000\,\rm km\,\rm
s^{-1}$, defining the time at which we start the simulations,

$t_{\rm max}=\frac{r_{\rm max}}{v_{\rm max}},$ (29)

and $r_{\mathrm{max}}$ is found using the numerical procedure of van Veelen et
al. (2009).

In the above relations, $E_{\rm ej}=5\times 10^{50}\,\rm erg$ is the energy
released by the supernova explosion, and,

$M_{\rm
ej}=M_{\star}-\int_{t_{\mathrm{ZAMS}}}^{t_{\mathrm{SN}}}\dot{M}(t)~{}dt-
M_{\mathrm{NS}}=10.12\,\rm M_{\odot},$ (30)

is the ejecta mass, with $t_{\mathrm{ZAMS}}$ and $t_{\mathrm{SN}}$ the zero-
age main-sequence and supernova times of the evolution of the $35\,\rm
M_{\odot}$ that we consider, and $M_{\mathrm{NS}}=1.4\,\rm M_{\odot}$ the mass
of the neutron star left behind the dead progenitor.

###  2.4 Modelling the supernova remnant

The 1D ejecta-wind interaction is mapped onto a 2.5D computational domain of
size $[0;45]\times[-45;45]\,\rm pc^{2}$ before the supernova shock wave
reaches the termination shock of the Wolf-Rayet stellar wind, which
corresponds to a circle of radius $2.45\,\rm pc$ centered onto the origin. The
mixing of the ejecta into the circumstellar medium at the pre-supernova time
is traced by a passive scalar $Q_{\rm EJ}$ advected via,

$\frac{\partial(\rho Q_{\rm EJ})}{\partial
t}+\mn@boldsymbol{\nabla}\cdot(\mn@boldsymbol{v}\rho Q_{\rm EJ})=0,$ (31)

and whose value is initialy set to $Q_{\rm EJ}=1$ in the ejecta at time
$t_{\rm max}$ and to $Q_{\rm EJ}=0$ in the rest of the circumstellar and
interstellar medium, respectively.

In this study, we focus on the case of a $35\,\rm M_{\odot}$ star that is
static or moving with velocity $v_{\star}=20$$-$$40\,\rm km\,\rm s^{-1}$ in a
uniform medium of number density $0.79\,\rm cm^{-3}$, which implies that the
supernova progenitor moves with a Mach number $M\approx 0$, $1$ and $2$,
respectively. The supernova remnants are calculated with a grid resolution of
$2000\times 4000$ grid zones, while we perform a resolution study with the
model with $v_{\star}=20\,\rm km\,\rm s^{-1}$. black For blackreference, we
compare our results with the supernova remnant model with $v_{\star}=0\,\rm
km\,\rm s^{-1}$ of Meyer et al. (2022). We list the ensemble of performed
models in Table 1.

Table 1:  List of models and spatial resolution adopted in each simulation.  ${\rm{Model}}$ | Grid mesh
---|---
Run-35-MHD-0-SNR | $3000\times 6000^{(a)}$
Run-35-MHD-20-SNR-1000 | $1000\times 2000$
Run-35-MHD-20-SNR-2000 | $2000\times 4000$
Run-35-MHD-20-SNR-4000 | $4000\times 8000$
Run-35-MHD-20-SNR-8000 | $8000\times 16000$
Run-35-MHD-40-SNR-2000 | $2000\times 4000$
(a) Meyer et al. (2022) |

Figure 3:  Number density fields in our magneto-hydrodynamical simulation of
the supernova remnant of the runaway $35\,\rm M_{\odot}$ star rotating with
$\Omega_{\star}/\Omega_{\rm K}=0.1$ and moving with velocity
$v_{\star}=20\,\rm km\,\rm s^{-1}$ (left) and $v_{\star}=40\,\rm km\,\rm
s^{-1}$ (right). The various contours highlight the region with a $50\%$
contribution of supernova ejecta (green), Wolf-Rayet wind (purple), red
supergiant wind (red), main-sequence material (blue), respectively. The black
lines are iso-temperature contours ($T=10^{6}\,\rm K$) and the black arrows
are magnetic field lines.

##  3 Results

This section presents the modelled core-collapse supernova remnant generated
by a Wolf-Rayet-evolving runaway massive progenitor star moving at velocity
$v_{\star}=20\,\rm km\,\rm s^{-1}$ and $v_{\star}=40\,\rm km\,\rm s^{-1}$. We
begin to present the pre-supernova circumstellar medium and the evolution of
the supernova remnants. We focus on the internal properties of the supernova
remnants, their luminosities and the mixing of material happening therein.

###  3.1 The pre-supernova circumstellar medium

Fig. 3 displays circumstellar medium at the supernova time in our MHD
simulations of massive runaway stars moving with velocity $v_{\star}=20\,\rm
km\,\rm s^{-1}$ (panel a) and $v_{\star}=40\,\rm km\,\rm s^{-1}$ (panel b).
The figures represent the number density field (in $\rm cm^{-3}$) in the
models and their right-hand part plots the ISM magnetic field lines (black
arrows) vertical to the symmetry axis and penetrating the region of shocked
material. The density field is structured as a large-scale bow shock generated
during the main-sequence phase of the star’s life. It has the standard
features of the astropshere of moving stellar objects, with a bow shock facing
the direction of stellar motion, separated from the star by a stand-off
distance (Baranov et al., 1971),

$R_{\rm
SO}=\sqrt{\frac{\dot{M}v_{\mathrm{w}}}{4\pi\rho_{\mathrm{ISM}}v_{\star}^{2}}},$
(32)

with $\rho_{\mathrm{ISM}}$ the mass density of the ISM material, and which
overall shape can be approximated by,

$R({\theta})=R_{\rm
SO}\frac{\sqrt{3(1-\theta)\mathrm{cotan}(\theta)}}{\sin(\theta)},$ (33)

where $\theta$ is the angle from the direction of stellar motion (Wilkin,
1996). A tubular low-density region shaped by the winds and a central circular
cavity surrounded by a Wolf-Rayet shell are found at the apex of the bow
shock. A series of contours trace the regions made of $50\%$ in number density
of each kind of material, such as the main-sequence stellar wind (blue), the
red supergiant (red) and the Wolf-Rayet (purple) stellar winds as well as the
supernova ejecta (green).

The cone of the bow shock behind the star is more opened in the case of the
star moving with $v_{\star}=20\,\rm km\,\rm s^{-1}$ (Fig. 3a) than in the case
of the faster-moving star (Fig. 3b). The Wolf-Rayet wind, faster than the
speed of stellar motion $v_{\star}$, forms a shell of evolved stellar wind and
swept-up ISM gas when interacting bow shock (Brighenti & D’Ercole, 1995a, b).
The circumstellar medium of the faster-moving star is affected by Rayleigh-
Taylor instabilities, both ahead of the bow shock where large eddies are seen,
and in the tubular cavity in the wake of the star, which close it and separate
the main-sequence material from the other winds (Fig. 3b). The environment of
the slower-moving stellar object displays a region of red supergiant material
(red) separating the main-sequence wind (blue) from the Wolf-Rayet material
(purple). Hence, the mixing of post-main-sequence materials is also at work
before the supernova explosion when the star moves slowly, because of the size
and unstable character of the low-density channel (Fig. 3a). Conversely, the
supernova ejecta of the fast-moving star are released in a closed bubble of
radius $\approx 10\,\rm pc$ centered onto the location of the defunct star
(Fig. 3a), whereas the direct environment of the supernova of the slowly
moving star is open to the trail of the bow shock (Fig. 3a).

###  3.2 Supernova remnants

The next panels of Fig. 3 show the evolution of the supernova remnants at
times $8\,\rm kyr$ ( middle line) and at time $80\,\rm kyr$ for the model with
$v_{\star}=20\,\rm km\,\rm s^{-1}$ (left) and $v_{\star}=40\,\rm km\,\rm
s^{-1}$ (right). They plot the reverberation of the supernova shock wave
against the walls of the asymmetric Wolf-Rayet cavity at work in the remnants,
together with the subsequent game of reflections and transmissions of waves
happening in it. The imprint of the morphology of that Wolf-Rayet cavity
already strongly modifies the geometry of the expanding forward shock at time
$8\,\rm kyr$. The ejecta filled the region between the center of the explosion
and the former termination shock of the pre-supernova bow shock, and the
reverberated shock towards the center of the explosion carries the imprint of
the distribution of Wolf-Rayet wind, see black contour for the region of gas
hotter than $10^{6}\,\rm K$ in Fig. 3c and Fig. 3d. The red supergiant
material (red) is distributed as a thin layer surrounding the ejecta region
(green), and it has started to efficiently mix with the Wolf-Rayet material
(purple) in the direction opposite of that of the progenitor’s motion. Note
that the red supergiant material distributes, because of the morphology of the
instabilities of the pre-supernova circumstellar medium, as knots in the case
of our faster-moving progenitor (Fig. 3d).

The bottom series of panels in Fig. 3 display the supernova remnants at a
later time of their evolution. The reverberated supernova shock wave
propagated all through the cavity and experienced a second reflection against
the walls of the cavity. The overall remnant has grown in size and a top
region of shocked ejecta and ISM gas expanding into the unperturbed ISM is
clearly visible. It is more pronounced in the model with $v_{\star}=40\,\rm
km\,\rm s^{-1}$ as a result of the smaller amount of mass in the Wolf-Rayet
shell (Meyer et al., 2015). The old supernova remnant is shaped as a Cygnus-
loop-like nebulae in the model with $v_{\star}=20\,\rm km\,\rm s^{-1}$ as a
result of the channelling of the supernova shock wave into the tunnel of low-
density cavity, while that generated by the progenitor moving with
$v_{\star}=40\,\rm km\,\rm s^{-1}$ harbors a more ovoid morphology. The
external envelope of the supernova remnant does not develop large scale
instabilities as a result of the background ISM magnetic field (Meyer et al.,
2021b). The Wolf-Rayet material is entirely mixed with the other kind of
material, except in some very small punctual regions (Fig. 3f), and the red
supergiant material has efficiently mixed with the supernova ejecta. The main-
sequence material is advected downstream the remnant and does not contribute
to its composition.

Figure 4:  Energy evolution in the supernova remnant (thick lines) as a
function of time, for different simulations with velocity $v_{\star}=20$
(left) and $40\,\rm km\,\rm s^{-1}$ (right). The figure displays the kinetic
(top panel), magnetic (middle panel) and potential components (bottom panel)
of the total energy. Each panel distinguish between the contribution from the
stellar wind (red solid line), main-sequence material (dashed blue line), red
supergiant material (dashed dotted line), Wolf-Rayet gas (solid blue line),
supernova ejecta (solid orange line) and shocked ISM gas (dashed dotted black
line). The thin lines give the corresponding percentages of the total energy
at a function of time.

###  3.3 Energetics

Fig. 4 displays the kinetic (top), magnetic (middle) and potential energy
(bottom) energies in the supernova remnants of progenitors moving with
velocity $v_{\star}=20\,\rm km\,\rm s^{-1}$ (left) and $v_{\star}=40\,\rm
km\,\rm s^{-1}$ (right). Each figure distinguishes between the components from
the main-sequence (dashed blue line), red supergiant (dotted dashed blue
line), Wolf-Rayet (solid blue line) stellar winds, the supernova ejecta
(orange), and the ISM material (dashed dotted black line). The kinetic energy
of a given species $i$ is calculated as,

$E_{\rm kin}^{i}=\mu m_{\rm
H}\iint_{\mathrm{SNR}}n|\vec{v}|^{2}Q_{i}\mathrm{dV},$ (34)

by integrating the kinetic energy density over the entire supernova remnant,
with $Q_{i}$ the tracer associated to the considered species and $dV$ the
volume element. Similarly, the magnetic energy of the gas is defined as,

$E_{\rm mag}^{i}=\iint_{\mathrm{SNR}}\frac{|\vec{B}|^{2}}{2}Q_{i}\mathrm{dV},$
(35)

and, finally, the potential energy in the supernova remnant reads,

$E_{\rm int}^{i}=\mu m_{\rm H}\iint_{\mathrm{SNR}}n\epsilon
Q_{i}\mathrm{dV}=\mu m_{\rm
H}\iint_{\mathrm{SNR}}n\frac{p}{(\gamma-1)}Q_{i}\mathrm{dV},$ (36)

where $\epsilon=p/(\gamma-1)$ is the internal energy density. These quantities
are evaluated throughout the supernova remnant simulation, up to a time
$40\,\rm kyr$.

Figure 5:  Soft X-ray (bottom) and hard X-ray (top) lightcurves of the
supernova remnant. Each panel distinguishes the contribution from the stellar
wind (red solid line), main-sequence material (dashed blue line), red
supergiant material (dashed dotted line), Wolf-Rayet gas (solid blue line),
supernova ejecta (solid orange line) and shocked ISM gas (dashed dotted black
line). We also show the time evolution of the relative contributions to the
energy budget for each species in the simulation (thin lines, in per cent).

In both models, the kinetic component of the energy budget of the supernova
remnant ($E_{\rm kin}\sim 10^{50}\,\rm erg$) is more important than that of
the magnetic ($E_{\rm mag}\sim 10^{49}\,\rm erg$) and internal potential
($E_{\rm pot}\sim 10^{26}\,\rm erg$) energy. The decreasing evolution of the
ratio $E_{\rm kin}/E_{\rm mag}$ that in middle-age ($30$$-$$35\,\rm kyr$)
core-collapse supernova remnants $E_{\rm mag}$ accounts for about $10\%$ of
the energy of the supernova remnant and should not be neglected into the
simulations. The kinetic energy of the supernova remnant whose progenitor
moves with velocity $v_{\star}=20\,\rm km\,\rm s^{-1}$ is initially such that
$E_{\rm kin}^{\rm ISM}\sim E_{\rm kin}^{\rm ej}$ and the stellar wind does not
contribute much to the energy budget $E_{\rm kin}^{\rm wind}\ll E_{\rm
kin}^{\rm ej}+E_{\rm kin}^{\rm ISM}$, see Fig. 4a. The wind contribution is
governed by the red supergiant material. Then, the ejecta kinetic energy
diminishes because of the deceleration of the blastwave which enters the
Sedov-Taylor phase once it interacts with the circumstellar medium of the
progenitor and sweeps up more of that material than the mass of the ejecta.
The situation is similar in the model with progenitor’s velocity
$v_{\star}=40\,\rm km\,\rm s^{-1}$, although we note an energetically more
important main-sequence wind in the first $15\,\rm kyr$ of the remnant’s life
(Fig. 4d). This is due to the thinner, more unstable cavity of unshocked
stellar wind left behind the progenitor in the case of a fast-moving star,
which prevents the ejecta to interact directly with the main-sequence
material, see Fig. 3.

The magnetic energy $E_{\rm mag}$ of the gas constituting the supernova
remnant is about two orders of magnitude smaller than the kinetic energy
$E_{\rm kin}$. Within the many materials inside of the supernova remnant, the
ejecta contributes less than one per cent to $E_{\rm mag}$, because in our
models the unshocked supernova material is initially unmagnetised (Section
2.3). The shocked ejecta compress the stellar wind and the circumstellar
magnetic field lines, which mix with the ejecta as the blastwave does through
the material of the wind-blown bubble. The contribution from the stellar wind
is governed by the red supergiant material, while the main-sequence gas, and,
to an even larger extend the Wolf-Rayet material are negligible to the
magnetic energy ($E_{\rm mag}^{\rm RSG}>E_{\rm mag}^{\rm MS}$ and $E_{\rm
mag}^{\rm RSG}>E_{\rm mag}^{\rm WR}$). The denser surroundings generated by
the progenitor moving with velocity $v_{\star}=40\,\rm km\,\rm s^{-1}$ induce
a slightly more magnetised remnant (Fig. 4e). The internal contribution to the
total energy of the supernova remnant is very small (Fig. 4c,f). It is
regulated by the shocked ISM material ($90\%$), the red supergiant material
weights for $1\%$ of the internal energy while the other winds (main-sequence
and Wolf-Rayet) have negligible contributions to $E_{\rm int}$.

###  3.4 X-rays luminosities

Fig. 5 displays the hard (top) and soft (bottom) X-ray lightcurves of the
supernova remnants of the progenitors moving with velocity $v_{\star}=20\,\rm
km\,\rm s^{-1}$ (left) and $v_{\star}=40\,\rm km\,\rm s^{-1}$ (right). As in
Fig. 4, the different stellar winds, the supernova ejecta, and the ISM
material are separated using the system of passive scalar tracers. For a given
species, the X-ray luminosity at energy $E\geq\alpha\,\rm keV$ is estimated
as,

$L_{X_{\textcolor{black}{\mathrm{E}\geq\alpha}}^{i}}=\iint_{\mathrm{SNR}}j_{\mathrm{X}}^{\textcolor{black}{\mathrm{E}\geq\alpha}}(T)n_{\mathrm{H}}^{2}Q_{i}\mathrm{dV},$
(37)

where $j_{\mathrm{X}}$ is the emissivity for diffuse ISM obtained with the
xspec software (Arnaud, 1996). Therefore, the luminosity in a specific energy
band reads,

$L_{X_{\textcolor{black}{\alpha\leq\mathrm{E}\leq\beta}}^{i}}=\iint_{\mathrm{SNR}}\Big{(}j_{\mathrm{X}}^{\mathrm{E}\leq\beta}(T)-j_{\mathrm{X}}^{\alpha\leq\mathrm{E}}(T)\Big{)}n_{\mathrm{H}}^{2}Q_{i}\mathrm{dV},$
(38)

and we perform such integral for the $0.1\leq E\leq 5.0\,\rm keV$ (soft) and
$E\,\geq 5.0\,\rm keV$ (hard).

The overall hard X-ray luminosities, dominated by the ejecta and the shocked
ISM, abruptly increase when the supernova shock wave hits the circumstellar
material, before decreasing with time as the gas temperature lessen (thick
lines of Figs. 5a,b). Further mild variations at time $15\,\rm kyr$ mirror the
reverberation of the ejecta in the center of the remnants. The hard stellar-
wind contribution in the model with $v_{\star}=20\,\rm km\,\rm s^{-1}$
accounts for $<10\%$ of the soft $L_{\rm X}$, while the ISM and the ejecta
represent $\approx 40\%$ each. In the model with $v_{\star}=40\,\rm km\,\rm
s^{-1}$, the wind contribution is more important, especially at the moment of
the reflection of the shock wave inside of the remnant, when it can reach up
to $\approx 40\%$ and exceed that of the ISM plus ejecta (thin lines of Figs.
5a,b). On the contrary, the total soft X-ray luminosity increases and reaches
a plateau largely governed by the radiation from ISM material which exceeds
the luminosity in hard X rays by several orders of magnitude. Note that the
ejecta contribution to the soft X-ray luminosity is about $20\%$ as its
temperature is too hot ($T\geq 10^{6}\,\rm K$) and that the stellar wind
contribution is nebligible once the shock wave has interacted with the
circumstellar medium (at times $5\,\rm kyr$).

Interestingly, the main-sequence wind contribution in the hard X-ray band is
the smallest of all winds, while it is that of the Wolf-Rayet wind which is
the least emitting in the soft X rays. In both wavebands, the red supergiant
material dominates the stellar wind contributions to the total luminosity.
This can be explained since the main-sequence wind is rapidly advected
downstream in the ejecta-wind interaction region by the stellar motion. The
Wolf-Rayet material efficiently mixes inside of the supernova remnant up to
being confined in small regions scattering inside of it (Fig. 3).

Figure 6:  The volume (in $\rm pc^{3}$) occupied by the different materials
inside of the supernova remnant of the slower (top) and the faster progenitor
(bottom), plotted as a function of time. Each panel distinguishes between the
contribution from the stellar wind (red solid line), main-sequence material
(dashed blue line), red supergiant material (dashed dotted line), Wolf-Rayet
gas (solid blue line), supernova ejecta (solid orange line) and shocked ISM
gas (dashed dotted black line). The contributions to the energy budget for
each species in the simulation (thin lines, in per cent) is added to the plot.

###  3.5 Mixing of materials in the supernova remnant

Fig. 6 plots the volume evolution of the remnants whose progenitors are moving
with velocity $v_{\star}=20\,\rm km\,\rm s^{-1}$ (top) and $v_{\star}=40\,\rm
km\,\rm s^{-1}$ (bottom). Each figure distinguishes the components from the
main-sequence (dashed blue line), red supergiant (dotted dashed blue line),
Wolf-Rayet (solid blue line) stellar winds, the supernova ejecta (orange), and
the ISM material (dashed dotted black line). Most of the supernova remnants
are filled with shocked ISM material (see dashed dotted think black lines in
Fig. 6a,b) which occupies about $80$$-$$90\%$ of its volume, all the other
species representing about $10$$-$$20\%$ of the remaining volume. The volume
of the remnant is constant in our simulation because it is calculated taking
into account pre-supernova circumstellar medium, which does not evolve much
after the time of the explosion. However, the volume filled by the ejecta
quickly increases after the explosion during the freely expanding phase of the
blastwave ($\leq 5\,\rm kyr$ after the explosion) and its proportion of
occupied volume reaches $\approx 15\%$ in the model with $v_{\star}=20\,\rm
km\,\rm s^{-1}$ (Fig. 6a) and $\approx 10\%$ $v_{\star}=40\,\rm km\,\rm
s^{-1}$ (Fig. 6b) since the cavity of the slowly-moving progenitor is larger
than that of the fast-moving one (Fig. 3). The other stellar wind materials,
which very initially fill more space than the ejecta, quickly shrink as the
blastwave propagates and mixes with the Wolf-Rayet stellar wind. At time $\leq
10\,\rm kyr$ this induces a transition between the wind material in the
remnant being mostly made of Wolf-Rayet gas to a stage in which the red
supergiant gas dominates.

Figure 7:  black Mixing of the different materials in the supernova remnant
blackduring the first $50\,\rm kyr$ after the supernova explosion. The line
colors distinguish models with stellar velocity $v_{\star}=0\,\rm km\,\rm
s^{-1}$ (solid orange, see Meyer et al. (2022)), $v_{\star}=20\,\rm km\,\rm
s^{-1}$ (dotted red) and $v_{\star}=40\,\rm km\,\rm s^{-1}$ (dotted dashed
blue). The blackpanels are for the main-sequence (a), red supergiant (b),
Wolf-Rayet (c) blackwind material and for the supernova ejecta (d),
respectively.
Figure 8:  black Probability density distribution of the gas number density
(in $\rm cm^{-3}$) in the supernova remnant model with $v_{\star}=20\,\rm
km\,\rm s^{-1}$, plotted for various times from the supernova explosion up to
$44\,\rm kyr$. The figures distinguish the different kind of materials in the
supernova remnants, such as the stellar wind material (red) and the ejecta
(green). The circumstellar material of the progenitor star is further divided
into main-sequence (solid line), red supergiant (dotted line) and Wolf-Rayet
materials (dashed line), respectively.

Fig. 7 displays the mixing factor for the different kind of materials present
in the supernova remnants of the progenitors moving with $v_{\star}=20\,\rm
km\,\rm s^{-1}$ (dashed red line) and with $v_{\star}=40\,\rm km\,\rm s^{-1}$
(dotted dashed blue line). We define the mixing factor of the main sequence,

$f_{\rm MS}=\frac{\mu m_{\rm H}\iint_{Q_{\rm MS}\leq 0.5}nQ_{\rm
MS}\mathrm{dV}}{\mu m_{\rm H}\iint_{\mathrm{SNR}}nQ_{\rm
MS}\mathrm{dV}}=\frac{\iint_{Q_{\rm MS}\leq 0.5}nQ_{\rm
MS}\mathrm{dV}}{\iint_{\mathrm{SNR}}nQ_{\rm MS}\mathrm{dV}},$ (39)

the red supergiant,

$f_{\rm RSG}=\frac{\iint_{Q_{\rm RSG}\leq 0.5}nQ_{\rm
RSG}\mathrm{dV}}{\iint_{\mathrm{SNR}}nQ_{\rm RSG}\mathrm{dV}},$ (40)

and the Wolf-Rayet

$f_{\rm WR}=\frac{\iint_{Q_{\rm WR}\leq 0.5}nQ_{\rm
WR}\mathrm{dV}}{\iint_{\mathrm{SNR}}nQ_{\rm WR}\mathrm{dV}},$ (41)

stellar winds, respectively. The mixing factor of the ejecta is similarly
defined,

$f_{\rm ej}=\frac{\iint_{Q_{\rm ej}\leq 0.5}nQ_{\rm
ej}\mathrm{dV}}{\iint_{\mathrm{SNR}}nQ_{\rm ej}\mathrm{dV}},$ (42)

where the numerator is the mass with subdominant ejecta fraction and where the
denominator is the initial mass of that same material, see Orlando et al.
(2005). For $f=0$ no mixing at all has taken place inside of the remnant while
for $f=1$ most ejecta material is dispersed into volume elements to which it
contributes $<50\%$ in number density.

black At the beginning of the simulation, the main-sequence material is
already fairly mixed to a fraction $\geq 95\%$, regardless of the motion of
the progenitor star (Fig. 7a). blackWe consider only the circumstellar medium
at distance $<50\,\rm pc$ from the location of the explosion and use it to
initialise the simulations, which ignores the long trail of main-sequence
material behind the star (Fig.3). Hence, the main-sequence material in the
computational domain only represents that expelled by the star before the
onset of the red supergiant phase. The mixing factor increases $15-20$ kyr
after the explosion, when the ejecta, channelled into the low-density cavity,
pushes the main-sequence wind out of the box, and $f_{\rm MS}=1$ for
$v_{\star}=20\,\rm km\,\rm s^{-1}$, see Fig. 7a. Conversely, the red
supergiant material is distributed between the Wolf-Rayet gas and the main-
sequence bow shock, and it is compressed and mixed as the blastwave
propagates. Since the ring produced by the Wolf-Rayet material is out of its
stellar wind bubble, the space velocity of the progenitor star matters very
little, and the evolution of $f_{\rm RSG}$ is similar in both cases, see Fig.
7b.

black The same Wolf-Rayet wind develops at early times in our supernova
remnants and the ovoidal ring it carves is similar in both cases, and so is
the circumstellar distribution of Wolf-Rayet material. The ring is hit by a
symmetric supernova blastwave resulting from a spherical explosion, and the
mixing factor $f_{\rm WR}$ blackevolves similarly in both models. The material
is rapidly mixed, and blackeventually the Wolf-Rayet gas is distributed in
small islands and clumps (Fig. 7c). The ejecta mixing factor in the supernova
remnant with velocity $v_{\star}=20\,\rm km\,\rm s^{-1}$ stays $f_{\rm
ej}<0.15$ during the early $50\,\rm kyr$ of its evolution and evolves little
over the first $20\,\rm kyr$, because once the shock wave hits the
circumstellar medium and is slowly reflected towards the center of the
location (Fig 3), no further mixing happens. In the model with velocity
$v_{\star}=40\,\rm km\,\rm s^{-1}$, these reflections of the blastwave are
facilitated by the bow shock, which permits the shock wave to rapidly undergo
further propagation to the center of the remnant, inducing more mixing of the
supernova ejecta with the stellar wind (Fig 7d). At time $30\,\rm kyr$ the
mixing factor has reached $f_{\rm ej}=0.4$ to further increase as the waves
propagates into the remnant (Fig 3f). All in all, within a few tens of parsecs
around the center of the explosion, stellar motion principally affects the
mixing of supernova ejecta rather than that of the stellar winds.

black A comparison blackof the mixing of material in supernova remnants
between static and runaway progenitor stars is plotted as a yellow solid line
in Fig. 7a-d. When the progenitor star does not move but evolves under the
influence of the magnetic field of the ISM, the region of unshocked wind
bordered by the wind/ISM contact discontinuity is elongated along the
direction of the local ISM magnetic field, see (Meyer et al., 2022). The main
difference blackbetween static and runaway models lies in the blackpaucity of
($<70\%$) mixing of main-sequence material blackfor $v_{\star}=0\,\rm km\,\rm
s^{-1}$ (Fig. 7a). This also applies to the blackred supergiant and Wolf-Rayet
wind material, see Fig. 7b,c, blackwhose mixing efficiencies are much below
that of the remnant models with runaway progenitors. The situation is
different regarding the supernova ejecta, as its mixing efficiency is larger
in the model with $v_{\star}=0\,\rm km\,\rm s^{-1}$ compare to that of with
$v_{\star}=20\,\rm km\,\rm s^{-1}$ for times $>12\,\rm kyr$ (Fig. 7d). The
reason for this is the proximity of the center of the supernova explosion with
the shocked wind regions as an effect of the stretched contact discontinuity
of stellar wind bubbles around static massive stars, see also van Marle et al.
(2015).

Figure 9:  black Comparison of the mixing of different materials in the
supernova remnants for progenitor velocity $v_{\star}=20\,\rm km\,\rm s^{-1}$
(top) and $v_{\star}=40\,\rm km\,\rm s^{-1}$ (bottom) at the time of the
supernova explosion and $16\,\rm kyr$ later. The figures distinguish the
stellar wind material (red) and the ejecta (green). The circumstellar material
of the progenitor star is further divided into main-sequence (solid line), red
supergiant (dotted line) and Wolf-Rayet materials (dashed line), respectively.

black Fig. 8 shows the blackvolumetric distribution of the gas number density
of the supernova remnant, plotted at specific times from the supernova
explosion to the age of $44\,\rm kyr$. The figures distinguishe the stellar
wind material (red) and the ejecta (green). The circumstellar material of the
progenitor star is further divided in main-sequence (solid line), red
supergiant (dotted line) and Wolf-Rayet (dashed line) materials, respectively.
The blackejecta have a broad density distribution. It is already interacting
with the stellar wind and blackdisplays a two-components distribution, while
the stellar wind distribution blackhas a single peak that is dominated by
main-sequence gas. The red supergiant gas already exhibits a rather complex
distribution, blackreflecting the complicated and knotty morphology of the
circumstellar medium in which the star died (Fig. 8a). At time $2\,\rm kyr$,
the ejecta density distribution has evolved as a result of the interaction
with the stellar wind and its expansion (Fig. 8b). At time $4\,\rm kyr$ both
the stellar wind and the supernova ejecta distributions adopt the shape of a
two-component structure. This is due to the interaction of the supernova shock
wave with the Wolf-Rayet material, whose number density increases. Similarly,
the red supergiant material distribution splits into two components, and
blackthe two peaks become clearly separated when the supernova ejecta are
channelled into the low-density cavity of the circumstellar medium. As the
blastwave further interacts with the former bow shock, the density of shocked
Wolf-Rayet gas further increases and the distribution of stellar wind blackgas
shifts to the higher densities (Fig. 8d). At time $25\,\rm kyr$ the red
supergiant material is efficiently mixed a pattern inside of the supernova
remnant, leading to multiple components in the red supergiant wind
distribution (Fig. 8e), while the Wolf-Rayet material globally increases in
density. At later time ($44\,\rm kyr$) only some Wolf-Rayet material remains
unmixed, the other kinds of stellar wind being largely melted with each other.
At this time, the supernova ejecta are almost entirely dispersed into the
remnant.

black Fig. 9 compares the distribution of the gas density at selected times in
the evolution of supernova remnants generated by a progenitor moving with
velocity $v_{\star}=20\,\mathrm{km\,s}^{-1}$ (top) and
$v_{\star}=40\,\mathrm{km\,s}^{-1}$ (bottom). The figures distinguish the
stellar wind material (red) and the ejecta (green). The circumstellar material
of the progenitor star is further divided in main-sequence, red supergiant and
Wolf-Rayet materials, respectively. At the moment of the explosion, the
distribution of ejecta does not depend on the progenitor speed since we used
the same explosion properties for both models (Fig. 9 a,c), differences being
exclusively in the distribution of the circumstellar material and shocked ISM.
The stellar wind material (especially red supergiant) already exhibits a
different level of mixing because a faster motion of the star induces more
mixing. blackIn the model with $v_{\star}=20\,\rm km\,\rm s^{-1}$, this
dichotomy increases at later times, when the red supergiant material enters
the low-density cavity and the density distribution of the stellar wind and
the gas distribution splits into two components (Fig. 9b), whereas the
distribution of Wolf-Rayet gas blackappears similar in both models (green
histograms of Fig. 9b,d). The main difference in the density distribution is
the manner the main-sequence and red supergiant gas melt into the supernova
remnant with the other kind of materials, see solid and dotted lined
histograms of Fig. 9b,d. The slightly more peaked distribution of supernova
ejecta in the fast-moving model reflects the different reverberation
properties of the supernova shock wave against the circumstellar medium (Fig.
9b,d).

Figure 10:  black Effect of spatial resolution on the mixing of the different
materials in the supernova remnants, plotted as a function of time (in $\rm
kyr$) for the models with stellar velocity $v_{\star}=20\,\rm km\,\rm s^{-1}$.
The panels show the mixing efficiency in simulations performed with spatial
resolutions $1000\times 2000$ (dashed red line), $2000\times 4000$ (dotted
dashed blue line), $4000\times 8000$ (solid orange line) and $8000\times
16000$ (solid black line). The figures plot the evolution of the mixing
efficiency during the early $50\,\rm kyr$ after the supernova explosion for
the main-sequence (a), red supergiant (b), Wolf-Rayet (c) stellar winds and
for the supernova ejecta (d), respectively.
Figure 11:  black Mixing of material in the supernova remnants generated by a
progenitor moving with velocity $20\,\rm km\,\rm s^{-1}$ (left) and $40\,\rm
km\,\rm s^{-1}$ (right), shown at times $8\,\rm kyr$ (top) and $80\,\rm kyr$
(bottom). The figures plot the value of the quantity $Q_{\rm RSG}$
representing the proportion of red supergiant (left part of the panels) and
the quantity $Q_{\rm WR}$ representing the proportion of Wolf-Rayet (right
part of the panels) materials, respectively. The blue contours of the passive
scalar field $Q_{\rm EJ}$ indicate a 10 per cent contribution of ejecta in
number density.

##  4 Discussion

In this last section, we recall the limitation of the models and discuss our
models in the context of available observations.

###  4.1 Caveats of the method

Although we introduce stellar rotation and magnetised stellar wind flows in
our simulations of core-collapse supernova remnants, our models are
nevertheless affected by a number of limitations. The most important of them
is the two-dimensional nature of the calculations, potentially impacting on
the stability of the wind-ISM interaction. Three-dimensional simulations can
solve this issue, however at much lower spatial resolution (Katushkina et al.,
2017, 2018; Gvaramadze et al., 2018; Meyer et al., 2021c). The treatment of
the ISM is also an element which can be subject to improvements, especially in
terms of intrinsic turbulence and granulosity of the ambient medium in which
the progenitor stars evolve (Martizzi et al., 2015; Haid et al., 2016).
Similarly, the interaction of the supernova shock wave with dense molecular
regions (McKee & Cowie, 1975; van Dishoeck et al., 1993) or low-density parts
(Arias et al., 2019) of the ISM which might modify the propagation of the
forward shock of the expanding supernova should also be investigated in future
studies as this may be a key factor in the shaping of middle-age supernova
remnants (Ferreira & de Jager, 2008; Rathjen et al., 2021). Last, several
physical mechanisms must be incorporated into the simulations of our study,
such as anisotropic heat conduction (Bedogni & D’Ercole, 1988; Velázquez et
al., 2004; Orlando et al., 2005; Vieser & Hensler, 2007; Balsara et al.,
2008a; Balsara et al., 2008b) or an explicit treatment of the photoionization
of the progenitor star in order to simulate the degree of ionization of the
gas in the remnant (see e.g. García-Segura et al., 1999; van Marle et al.,
2004; Arthur & Hoare, 2006; Toalá & Arthur, 2011).

###  4.2 Effects of spatial resolution

Fig. 10 repeats the analysis presented in Fig. 7 by investigating the effect
of spatial resolution on the time evolution of the mixing of the different
materials in the supernova remnant models with stellar velocity
$v_{\star}=20\,\rm km\,\rm s^{-1}$. The panels show the mixing efficiency in
simulations performed with spatial resolutions $1000\times 2000$ (dashed red
line), $2000\times 4000$ (dotted dashed blue line), $4000\times 8000$ (dashed
dotted orange line) and $8000\times 16000$ (solid black line). The tracers in
the simulation are used to separate the mixing efficiency during the early
$50\,\rm kyr$ after the supernova explosion for the main-sequence (a), red
supergiant (b), Wolf-Rayet (c) stellar winds and for the supernova ejecta (d),
respectively. No notable differences in terms of mixing efficiency of the
main-sequence material happens, especially because most of this kind of gas
has left the computational domain at times $\geq 25\,\rm kyr$. black Note that
because of the definition that we adopted, the mixing factor equals to unity
once all the material has left the computational domain (Fig. 10a). The mixing
of the red supergiant material is the most sensitive to spatial resolution.
The variability of the mixing factor diminishes with the increasing spatial
resolution, and its value is slightly overestimated when the resolution is too
coarse (Fig. 10b). The models show that it has not yet converged, except at
longer times $\geq 45\,\rm kyr$, and that even higher-resolution simulations
are needed to carefully track the mixing of the dusty supergiant material into
supernova remnants. Inversely, the Wolf-Rayet gas mixes quickly into the
remnants, and has already converged in our low-resolution models (Fig. 10c).
The ejecta mixing behaves differently: the mixing efficiency increases by a
factor of $\approx 2$ when the resolution is multiplied by a factor of 4 (Fig.
10d). High-resolution is therefore a key ingredient of core-collapse supernova
mixing modelling (Orlando et al., 2019).

###  4.3 Spatial distribution of the mixed material

black Fig. 11 displays the distribution of mixed materials in the supernova
remnant generated by a progenitor moving with velocities $20\,\rm km\,\rm
s^{-1}$ (left panels) and $40\,\rm km\,\rm s^{-1}$ (right panels), shown at
times $8\,\rm kyr$ (top panels) and $80\,\rm kyr$ (bottom panels). The figures
plot the value of the quantity $Q_{\rm RSG}$ representing the proportion of
red supergiant (left part of the panels) and the quantity $Q_{\rm WR}$
representing the proportion of Wolf-Rayet (right part of the panels)
materials, respectively. The blue contours of the passive scalar field $Q_{\rm
EJ}$ indicate a 10 per cent contribution of ejecta in number density. At times
$8\,\rm kyr$ blackthere are few differences between the models. The enriched
gases are distributed as a ring surrounding the expanding supernova blastwave
that is interacting with the pre-supernova circumstellar medium (Fig. 11a,b).
The low-density cavity of the model with $20\,\rm km\,\rm s^{-1}$ mostly
contains red supergiant material channelled into it (Fig. 11a). When the
supernova remnant is much older ($80\,\rm kyr$), the manner the blastwave
reflected towards the center of the explosion also governs blackthe
distribution of the stellar winds. Massive bow shocks have stopped the
blastwave after its first reverberation (Meyer et al., 2015, 2021b), and both
the red supergiant and Wolf-Rayet material are located around the location
where the supernova shock wave interacted with the circumstellar medium,
producing a Cygnus-loop-like region filled with ejecta and post-main-sequence
materials (Fig. 11c). In the other model with $40\,\rm km\,\rm s^{-1}$, the
red supergiant gas fills the lower ($z<0$) region of the supernova remnant, in
which clumps of Wolf-Rayet-rich gas are scattered. This is a major difference
between remnants from slow- and fast-moving progenitors, which has observable
consequences.

###  4.4 Motivation of the choice of an initial $35\,\rm M_{\odot}$ progenitor

This work constitutes the method paper of a series of studies devoted to the
comprehensive modelling of supernova remnants of static and runaway massive
stars. A single, non-rotating model has already been published, presenting
thermal and non-thermal emission maps of the transient rectangular-like
supernova remnant generated when a supernova blastwave interacts with the
elongated termination shock of a stellar wind bubble forming around a massive
star located in a uniformly magnetised environment. It underlined how this
mechanism may participate in the shaping of the rectangular supernova remnant
Puppis A (Meyer et al., 2022). The choice of a $35\,\rm M_{\odot}$ progenitor
had been motivated by conclusions from observations of Puppis A (Paron et al.,
2008), although alternative scenarios with lower-mass scenarios are currently
considered (Velázquez et al., 2023). Our work continues this numerical effort
by exploring the bulk motion of the same Wolf-Rayet-evolving progenitor. The
zero-age main-sequence rotation rate of the massive star used to simulate our
supernova remnants is taken to $10\%$ of its break-up velocity. It is
necessary to include rotation in the simulations otherwise the toroidal
component of the stellar magnetic field vanishes and only the insignificantly
weak radial component will remain. This will let the unshocked supernova
ejecta unmagnetised and make the exercise of Fig. 4 impossible with respect to
the magnetic energy. Further models will be performed in the future in order
to continue investigating how the surroundings of defunct high-mass stars look
like, e.g. by considering different initial progenitor masses and/or
environments. A global picture regarding the mixing of enriched materials via
core-collapse supernova remnants will hence emerge.

This collection of remnants will serve as a basis for the detailed study of
several potential individual targets. Amongst others, one can mention a few
supernova remnants of constrained initial mass consistent with a formation
scenario involving a Wolf-Rayet progenitor: Kes 73, G292.2-0.5, G349.7+0.2,
G292.0+1.8, Kes 79, Kes 17,in the Milky Way, N132D in the Large Magellanic
Cloud or 1E0102.2-7219 in the Small Magellanic Cloud. Note that our precedent
work (Meyer et al., 2021b) has established that supernova progenitors evolving
through a Wolf-Rayet phase can produce bilateral and the shell-type remnants
such as G296.5+10.0 and CTB109, respectively, at least as seen in the radio
waveband by non-thermal emission. Such Galactic and extragalactic objects
which are the perfect targets for tailoring our results to specific core-
collapse supernova remnants. Last, let also note that the Cherenkov Telescope
Array (CTA) will help in discovering new core-collapse supernova remnants to
which understanding will require numerical models.

##  5 Conclusion

This paper investigates the evolution of young to middle-age core-collapse
supernova remnants generated by a $35\,\rm M_{\odot}$ massive runaway star
with particular emphasis on the mixing of materials. Our 2.5D magneto-
hydrodynamical numerical simulations account for the rotation (Ekström et al.,
2012) and magnetisation (Meyer et al., 2021c; Baalmann et al., 2020, 2021) of
the progenitor star that is moving through the warm, uniform, magnetised phase
of the interstellar medium of the Milky Way. A system of passive scalar
tracers permits us to follow the advection of the different kind of wind
materials throughout the pre-supernova evolution of the progenitor. The
simulations are conducted in the fashion of the precedent papers of this
series, by first modelling the detailed circumstellar medium of the progenitor
star at the pre-supernova time, and by releasing into it a separately
calculated solution for the ejecta-wind that is simulated up to $100\,\rm
kyr$, see (Meyer et al., 2015; Meyer et al., 2020a; Meyer et al., 2021b). The
high-resolution calculations are performed using the pluto code (Mignone et
al., 2007, 2012; Vaidya et al., 2018) that has been validated for
circumstellar medium and supernova remnants studies (Meyer et al., 2014,
2017).

It is found that the manner the stellar winds associated to the various
evolutionary phases of the progenitor, as well as the supernova ejecta,
distribute into the supernova remnant differs greatly depending on the bulk
motion of the star, the level of wind mixing increasing with it, i.e. faster-
moving progenitors generate older remnants hosting more complex, dense
filamentary internal patterns of mixed species, where wind mixes less but
ejecta does it better. Very little main-sequence material ($\leq 5\%$) is
present in the remnants, as it is channelled into the tubular low-density
cavity of wind produced by the progenitor’s stellar motion. The red supergiant
material is the least mixed species, while the Wolf-Rayet gas efficiently
mixes in the early $20\,\rm kyr$ after the explosion.

black Our results show that the spatial distribution of the post-main-sequence
stellar winds blackdepends on the progenitor’s spatial motion. At later times
after the supernova explosion ($>40\,\rm kyr$), the reflection of the
supernova blastwave spreads blackwind material to the interior of the remnant
nebula, distributing Wolf-Rayet clumps inside of it, if the progenitor’s
motion is important ($v_{\star}=40\,\rm km\,\rm s^{-1}$). The enriched stellar
wind material remains in the region of the pre-supernova circumstellar medium
for slower progenitors ($v_{\star}=20\,\rm km\,\rm s^{-1}$). Interestingly,
the mixing efficiency of ejecta is less important in the context of runaway
progenitors than in the case of a massive star (Meyer et al., 2022). The
spatial resolution is very important in the numerical study of mixing
processes into core-collapse supernova remnants, as a low resolution tends to
overestimate that of the red supergiant and Wolf-Rayet stellar winds and
underestimates that of the supernova ejecta. Indeed, the regions filled by the
winds turned to be isolate islands located around the shell made of the
circumstellar medium shocked by the blastwave, where the material is both
subsonic and superalfvenic. Such work should be extended to the stellar
clusters made of multiple massive stars (Badmaev et al., 2022).

## Acknowledgements

The authors acknowledge the North-German Supercomputing Alliance (HLRN) for
providing HPC resources that have contributed to the research results reported
in this paper. M. Petrov acknowledges the Max Planck Computing and Data
Facility (MPCDF) for providing data storage resources and HPC resources which
contributed to test and optimise the pluto code.

## Data availability

This research made use of the pluto code developed at the University of Torino
by A. Mignone (http://plutocode.ph.unito.it/). 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

  * Ahmadi et al. (2018) Ahmadi A., et al., 2018, A&A, 618, A46
  * Aloy & Obergaulinger (2021) Aloy M. Á., Obergaulinger M., 2021, MNRAS, 500, 4365
  * Arias et al. (2019) Arias M., et al., 2019, A&A, 622, A6
  * Arnaud (1996) Arnaud K. A., 1996, in Jacoby G. H., Barnes J., eds, Astronomical Society of the Pacific Conference Series Vol. 101, Astronomical Data Analysis Software and Systems V. p. 17
  * Arthur & Hoare (2006) Arthur S. J., Hoare M. G., 2006, ApJS, 165, 283
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Baalmann et al. (2020) Baalmann L. R., Scherer K., Fichtner H., Kleimann J., Bomans D. J., Weis K., 2020, A&A, 634, A67
  * Baalmann et al. (2021) Baalmann L. R., Scherer K., Kleimann J., Fichtner H., Bomans D. J., Weis K., 2021, A&A, 650, A36
  * Badmaev et al. (2022) Badmaev D. V., Bykov A. M., Kalyashova M. E., 2022, MNRAS, 517, 2818
  * Balsara et al. (2008a) Balsara D. S., Tilley D. A., Howk J. C., 2008a, MNRAS, 386, 627
  * Balsara et al. (2008b) Balsara D. S., Bendinelli A. J., Tilley D. A., Massari A. R., Howk J. C., 2008b, MNRAS, 386, 642
  * Bandiera et al. (2023) Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2023, MNRAS,
  * Baranov et al. (1971) Baranov V. B., Krasnobaev K. V., Kulikovskii A. G., 1971, Soviet Physics Doklady, 15, 791
  * Bedogni & D’Ercole (1988) Bedogni R., D’Ercole A., 1988, A&A, 190, 320
  * Benaglia et al. (2010) Benaglia P., Romero G. E., Martí J., Peri C. S., Araudo A. T., 2010, A&A, 517, L10
  * Blaauw (1961) Blaauw A., 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Blondin et al. (1996) Blondin J. M., Lundqvist P., Chevalier R. A., 1996, ApJ, 472, 257
  * Brighenti & D’Ercole (1995a) Brighenti F., D’Ercole A., 1995a, MNRAS, 273, 443
  * Brighenti & D’Ercole (1995b) Brighenti F., D’Ercole A., 1995b, MNRAS, 277, 53
  * Broersen et al. (2014) Broersen S., Chiotellis A., Vink J., Bamba A., 2014, MNRAS, 441, 3040
  * Brott et al. (2011a) Brott I., et al., 2011a, A&A, 530, A115
  * Brott et al. (2011b) Brott I., et al., 2011b, A&A, 530, A116
  * Castro et al. (2015) Castro N., et al., 2015, A&A, 581, A81
  * Castro et al. (2017) Castro N., et al., 2017, A&A, 597, L6
  * Chevalier (1982) Chevalier R. A., 1982, ApJ, 258, 790
  * 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
  * Chiotellis et al. (2012) Chiotellis A., Schure K. M., Vink J., 2012, A&A, 537, A139
  * Chiotellis et al. (2013) Chiotellis A., Kosenko D., Schure K. M., Vink J., Kaastra J. S., 2013, MNRAS, 435, 1659
  * Chiotellis et al. (2019) Chiotellis A., Boumis P., Derlopa S., Steffen W., 2019, arXiv e-prints, p. arXiv:1909.08947
  * Chiotellis et al. (2020) Chiotellis A., Boumis P., Spetsieri Z. T., 2020, Galaxies, 8, 38
  * Chiotellis et al. (2021) Chiotellis A., Boumis P., Spetsieri Z. T., 2021, MNRAS, 502, 176
  * Chita et al. (2008) Chita S. M., Langer N., van Marle A. J., García-Segura G., Heger A., 2008, A&A, 488, L37
  * Comerón & Kaper (1998) Comerón F., Kaper L., 1998, A&A, 338, 273
  * Crowther (2007) Crowther P. A., 2007, ARA&A, 45, 177
  * Das et al. (2022) Das S., Brose R., Meyer D. M. A., Pohl M., Sushch I., Plotko P., 2022, A&A, 661, A128
  * Davies (2017) Davies B., 2017, Philosophical Transactions of the Royal Society of London Series A, 375, 20160270
  * De Becker et al. (2017) De Becker M., del Valle M. V., Romero G. E., Peri C. S., Benaglia P., 2017, MNRAS, 471, 4452
  * Derlopa et al. (2019) Derlopa S., Boumis P., Chiotellis A., Steffen W., Akras S., 2019, arXiv e-prints, p. arXiv:1909.06131
  * 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., et al., 2012, A&A, 537, A146
  * Eldridge et al. (2006) Eldridge J. J., Genet F., Daigne F., Mochkovitch R., 2006, MNRAS, 367, 186
  * Fang et al. (2017) Fang J., Yu H., Zhang L., 2017, MNRAS, 464, 940
  * Ferreira & de Jager (2008) Ferreira S. E. S., de Jager O. C., 2008, A&A, 478, 17
  * Finn et al. (2016) Finn K., Bianco F. B., Modjaz M., Liu Y.-Q., Rest A., 2016, ApJ, 830, 73
  * Fossati et al. (2015) Fossati L., et al., 2015, A&A, 574, A20
  * 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. (1996a) Garcia-Segura G., Mac Low M.-M., Langer N., 1996a, A&A, 305, 229
  * Garcia-Segura et al. (1996b) Garcia-Segura G., Langer N., Mac Low M.-M., 1996b, 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
  * 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
  * Gies (1987) Gies D. R., 1987, ApJS, 64, 545
  * Green et al. (2019) Green S., Mackey J., Haworth T. J., Gvaramadze V. V., Duffy P., 2019, A&A, 625, A4
  * Groh et al. (2014) Groh J. H., Meynet G., Ekström S., Georgy C., 2014, A&A, 564, A30
  * Gull & Sofia (1979) Gull T. R., Sofia S., 1979, ApJ, 230, 782
  * Gvaramadze (2006) Gvaramadze V. V., 2006, A&A, 454, 239
  * Gvaramadze et al. (2012) Gvaramadze V. V., Weidner C., Kroupa P., Pflamm-Altenburg J., 2012, MNRAS, 424, 3037
  * Gvaramadze et al. (2018) Gvaramadze V. V., Alexashov D. B., Katushkina O. A., Kniazev A. Y., 2018, MNRAS, 474, 4421
  * Haid et al. (2016) Haid S., Walch S., Naab T., Seifried D., Mackey J., Gatto A., 2016, MNRAS, 460, 2962
  * Hamann (1992) Hamann W. R., 1992, in Heber U., Jeffery C. S., eds, , Vol. 401, The Atmospheres of Early-Type Stars. p. 87
  * Harten et al. (1983) Harten A., Lax P. D., van Leer B., 1983, SIAM Review, 25, 35
  * Herbst et al. (2020) Herbst K., et al., 2020, ApJ, 897, L27
  * Herbst et al. (2022) Herbst K., et al., 2022, Space Sci. Rev., 218, 29
  * Holland-Ashford et al. (2020) Holland-Ashford T., Lopez L. A., Auchettl K., 2020, ApJ, 889, 144
  * Hummer (1994) Hummer D. G., 1994, MNRAS, 268, 109
  * Ilee et al. (2016) Ilee J. D., Cyganowski C. J., Nazari P., Hunter T. R., Brogan C. L., Forgan D. H., Zhang Q., 2016, MNRAS, 462, 4386
  * Jardine et al. (2022) Jardine R., Powell J., Müller B., 2022, MNRAS, 510, 5535
  * Johnston et al. (2015) Johnston K. G., et al., 2015, ApJ, 813, L19
  * Kaper et al. (1997) Kaper L., van Loon J. T., Augusteijn T., Goudfrooij P., Patat F., Waters L. B. F. M., Zijlstra A. A., 1997, ApJ, 475, L37
  * Katsuda et al. (2018) Katsuda S., Takiwaki T., Tominaga N., Moriya T. J., Nakamura K., 2018, ApJ, 863, 127
  * Katushkina et al. (2017) Katushkina O. A., Alexashov D. B., Izmodenov V. V., Gvaramadze V. V., 2017, MNRAS, 465, 1573
  * Katushkina et al. (2018) Katushkina O. A., Alexashov D. B., Gvaramadze V. V., Izmodenov V. V., 2018, MNRAS, 473, 1576
  * Kervella et al. (2018) Kervella P., et al., 2018, A&A, 609, A67
  * Kifonidis et al. (2006) Kifonidis K., Plewa T., Scheck L., Janka H. T., Müller E., 2006, A&A, 453, 661
  * Kobulnicky et al. (2016) Kobulnicky H. A., et al., 2016, ApJS, 227, 18
  * Kobulnicky et al. (2017) Kobulnicky H. A., Schurhammer D. P., Baldwin D. J., Chick W. T., Dixon D. M., Lee D., Povich M. S., 2017, AJ, 154, 201
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Langer (2012) Langer N., 2012, ARA&A, 50, 107
  * Levesque (2010) Levesque E. M., 2010, New Astron. Rev., 54, 1
  * Mackey et al. (2012) Mackey J., Mohamed S., Neilson H. R., Langer N., Meyer D. M.-A., 2012, ApJ, 751, L10
  * Maeder & Meynet (2000) Maeder A., Meynet G., 2000, ARA&A, 38, 143
  * Martizzi et al. (2015) Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504
  * McKee & Cowie (1975) McKee C. F., Cowie L. L., 1975, ApJ, 195, 715
  * Meyer (2021) Meyer D. M. A., 2021, MNRAS, 507, 4697
  * Meyer & Meliani (2022) Meyer D. M. A., Meliani Z., 2022, MNRAS, 515, L29
  * Meyer et al. (2014) Meyer D. M.-A., Mackey J., Langer N., Gvaramadze V. V., Mignone A., Izzard R. G., Kaper L., 2014, MNRAS, 444, 2754
  * Meyer et al. (2015) Meyer D. M.-A., Langer N., Mackey J., Velázquez P. F., Gusdorf A., 2015, MNRAS, 450, 3080
  * Meyer et al. (2016) Meyer D. M.-A., van Marle A.-J., Kuiper R., Kley W., 2016, MNRAS, 459, 1146
  * Meyer et al. (2017) Meyer D. M.-A., Vorobyov E. I., Kuiper R., Kley W., 2017, MNRAS, 464, L90
  * Meyer et al. (2019) Meyer D. M. A., Vorobyov E. I., Elbakyan V. G., Stecklum B., Eislöffel J., Sobolev A. M., 2019, MNRAS
  * Meyer et al. (2020a) Meyer D. M. A., Petrov M., Pohl M., 2020a, MNRAS, 493, 3548
  * Meyer et al. (2020b) Meyer D. M. A., Oskinova L. M., Pohl M., Petrov M., 2020b, MNRAS, 496, 3906
  * Meyer et al. (2021a) Meyer D. M. A., Vorobyov E. I., Elbakyan V. G., Eislöffel J., Sobolev A. M., Stöhr M., 2021a, MNRAS, 500, 4448
  * Meyer et al. (2021b) Meyer D. M. A., Pohl M., Petrov M., Oskinova L., 2021b, MNRAS, 502, 5340
  * Meyer et al. (2021c) Meyer D. M. A., Mignone A., Petrov M., Scherer K., Velázquez P. F., Boumis P., 2021c, MNRAS, 506, 5170
  * Meyer et al. (2022) Meyer D. M. A., et al., 2022, MNRAS, 515, 594
  * 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
  * Moffat et al. (1998) Moffat A. F. J., Marchenko S. V., Seggewiss W., van der Hucht et al. K. A. e. a., 1998, A&A, 331, 949
  * Müller & Janka (2015) Müller B., Janka H. T., 2015, MNRAS, 448, 2141
  * Müller et al. (2012) Müller E., Janka H. T., Wongwathanarat A., 2012, A&A, 537, A63
  * Olmi & Bucciantini (2019) Olmi B., Bucciantini N., 2019, MNRAS, 490, 3608
  * Olmi & Bucciantini (2023) Olmi B., Bucciantini N., 2023, arXiv e-prints, p. arXiv:2301.12903
  * Orlando et al. (2005) Orlando S., Peres G., Reale F., Bocchino F., Rosner R., Plewa T., Siegel A., 2005, A&A, 444, 505
  * Orlando et al. (2019) Orlando S., et al., 2019, A&A, 622, A73
  * Orlando et al. (2020) Orlando S., et al., 2020, A&A, 636, A22
  * Osterbrock & Bochkarev (1989) Osterbrock D. E., Bochkarev N. G., 1989, Soviet Ast., 33, 694
  * Parker (1958) Parker E. N., 1958, ApJ, 128, 664
  * Paron et al. (2008) Paron S., Dubner G., Reynoso E., Rubio M., 2008, A&A, 480, 439
  * Peri et al. (2012) Peri C. S., Benaglia P., Brookes D. P., Stevens I. R., Isequilla N. L., 2012, A&A, 538, A108
  * Peri et al. (2015) Peri C. S., Benaglia P., Isequilla N. L., 2015, A&A, 578, A45
  * 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
  * Powell (1997) Powell K. G., 1997, An Approximate Riemann Solver for Magnetohydrodynamics. Springer Berlin Heidelberg, Berlin, Heidelberg, pp 570–583, doi:10.1007/978-3-642-60543-7_23, https://doi.org/10.1007/978-3-642-60543-7_23
  * Przybilla et al. (2016) Przybilla N., et al., 2016, A&A, 587, A7
  * Raga et al. (1997) Raga A. C., Mellema G., Lundqvist P., 1997, ApJS, 109, 517
  * Rathjen et al. (2021) Rathjen T.-E., et al., 2021, MNRAS, 504, 1039
  * Rozyczka & Franco (1996) Rozyczka M., Franco J., 1996, ApJ, 469, L127
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * 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
  * Shu (1992) Shu F. H., 1992, The physics of astrophysics. Volume II: Gas dynamics.
  * Smartt (2009) Smartt S. J., 2009, ARA&A, 47, 63
  * Smith (2014) Smith N., 2014, ARA&A, 52, 487
  * Szécsi et al. (2022) Szécsi D., Agrawal P., Wünsch R., Langer N., 2022, A&A, 658, A125
  * Telezhinsky et al. (2012a) Telezhinsky I., Dwarkadas V. V., Pohl M., 2012a, Astroparticle Physics, 35, 300
  * Telezhinsky et al. (2012b) Telezhinsky I., Dwarkadas V. V., Pohl M., 2012b, A&A, 541, A153
  * Telezhinsky et al. (2013) Telezhinsky I., Dwarkadas V. V., Pohl M., 2013, A&A, 552, A102
  * Tenorio-Tagle et al. (1990) Tenorio-Tagle G., Bodenheimer P., Franco J., Rozyczka M., 1990, MNRAS, 244, 563
  * Tenorio-Tagle et al. (1991) Tenorio-Tagle G., Rozyczka M., Franco J., Bodenheimer P., 1991, MNRAS, 251, 318
  * Toalá & Arthur (2011) Toalá J. A., Arthur S. J., 2011, ApJ, 737, 100
  * Toledo-Roy et al. (2014) Toledo-Roy J. C., Esquivel A., Velázquez P. F., Reynoso E. M., 2014, MNRAS, 442, 229
  * Truelove & McKee (1999) Truelove J. K., McKee C. F., 1999, ApJS, 120, 299
  * Tutone et al. (2020) Tutone A., et al., 2020, A&A, 642, A67
  * Uchida et al. (2009) Uchida H., Tsunemi H., Katsuda S., Kimura M., Kosugi H., 2009, PASJ, 61, 301
  * Underhill (1968) Underhill A. B., 1968, ARA&A, 6, 39
  * Vaidya et al. (2018) Vaidya B., Mignone A., Bodo G., Rossi P., Massaglia S., 2018, ApJ, 865, 144
  * Velázquez et al. (2004) Velázquez P. F., Martinell J. J., Raga A. C., Giacani E. B., 2004, ApJ, 601, 885
  * Velázquez et al. (2006) Velázquez P. F., Vigh C. D., Reynoso E. M., Gómez D. O., Schneiter E. M., 2006, ApJ, 649, 779
  * Velázquez et al. (2023) Velázquez P. F., Meyer D. M. A., Chiotellis A., Cruz-Álvarez A. E., Schneiter E. M., Toledo-Roy J. C., Reynoso E. M., Esquivel A., 2023, MNRAS, 519, 5358
  * Vieser & Hensler (2007) Vieser W., Hensler G., 2007, A&A, 472, 141
  * Vigh et al. (2011) Vigh C. D., Velázquez P. F., Gómez D. O., Reynoso E. M., Esquivel A., Matias Schneiter E., 2011, ApJ, 727, 32
  * Vink et al. (2000) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2000, A&A, 362, 295
  * 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
  * Vorobyov (2009) Vorobyov E. I., 2009, ApJ, 704, 715
  * Vorobyov (2010) Vorobyov E. I., 2010, ApJ, 723, 1294
  * Vorobyov & Basu (2006) Vorobyov E. I., Basu S., 2006, ApJ, 650, 956
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Weber & Davis (1967) Weber E. J., Davis Leverett J., 1967, ApJ, 148, 217
  * Weiler & Sramek (1988) Weiler K. W., Sramek R. A., 1988, ARA&A, 26, 295
  * Whalen et al. (2008) Whalen D., van Veelen B., O’Shea B. W., Norman M. L., 2008, ApJ, 682, 49
  * Wiersma et al. (2009) Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99
  * Wilkin (1996) Wilkin F. P., 1996, ApJ, 459, L31
  * Woosley & Weaver (1986) Woosley S. E., Weaver T. A., 1986, ARA&A, 24, 205
  * Yusef-Zadeh et al. (2003) Yusef-Zadeh F., Wardle M., Rho J., Sakano M., 2003, ApJ, 585, 319
  * Zavlin et al. (1998) Zavlin V. E., Pavlov G. G., Trumper J., 1998, A&A, 331, 821
  * del Valle et al. (2013) del Valle M. V., Romero G. E., De Becker M., 2013, A&A, 550, A112
  * van Dishoeck et al. (1993) van Dishoeck E. F., Jansen D. J., Phillips T. G., 1993, A&A, 279, 541
  * van Marle et al. (2004) van Marle A. J., Langer N., García-Segura G., 2004, in G. Garcia-Segura, G. Tenorio-Tagle, J. Franco, & H. W. Yorke ed., Revista Mexicana de Astronomia y Astrofisica Conference Series Vol. 22, Revista Mexicana de Astronomia y Astrofisica Conference Series. pp 136–139
  * 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. (2010) van Marle A. J., Smith N., Owocki S. P., van Veelen B., 2010, MNRAS, 407, 2305
  * van Marle et al. (2011) van Marle A. J., Meliani Z., Keppens R., Decin L., 2011, ApJ, 734, L26
  * van Marle et al. (2014) van Marle A. J., Decin L., Meliani Z., 2014, A&A, 561, A152
  * van Marle et al. (2015) van Marle A. J., Meliani Z., Marcowith A., 2015, A&A, 584, A49
  * van Veelen et al. (2009) van Veelen B., Langer N., Vink J., García-Segura G., van Marle A. J., 2009, A&A, 503, 495