# On the plerionic rectangular supernova remnants of static progenitors

Meyer D. M.-A.1, Meliani Z.2, Velázquez P. F.3, Pohl M.1,4 and Torres D.
F.5,6,7
1 Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-
Strasse 24/25, D-14476 Potsdam, Germany,
E-mail<EMAIL_ADDRESS>
2 Laboratoire Univers et Théories, Observatoire de Paris, Université PSL,
Université de Paris, CNRS, F-92190 Meudon, France
3 Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, CP
04510, Mexico City, Mexico
4 Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen,
Germany
5 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans
s/n, 08193 Barcelona, Spain
6 Institut d’Estudis Espacials de Catalunya (IEEC), Gran Capità 2-4, 08034
Barcelona, Spain
7 Institució Catalana de Recerca i Estudis Avançats (ICREA), 08010 Barcelona,
Spain

(Last updated 2020 June 10; in original form 2013 September 5)

###### Abstract

Pulsar wind nebulae are a possible final stage of the circumstellar evolution
of massive stars, where a fast-rotating, magnetised neutron star produces a
powerful wind that interacts with the supernova ejecta. The shape of these so-
called plerionic supernova remnants is influenced by the distribution of
circumstellar matter at the time of the explosion, itself impacted by the
magnetic field of the ambient medium responsible for the expansion of the
circumstellar bubble of the progenitor star. To understand the effects of
magnetization on the circumstellar medium and resulting pulsar nebulae, we
conduct 2D magneto-hydrodynamical simulations. Our models explore the impact
of the interstellar medium magnetic field on the morphology of a supernova
remnant and pulsar wind nebula that develop in the circumstellar medium of
massive star progenitor in the warm phase of the Milky Way’s interstellar
medium. Our simulations reveal that the jet-like structures formed on both
sides perpendicularly to the equatorial plane of the pulsar, creating complex
radio synthetic synchrotron emissions. This morphology is characterized by a
rectangular-like remnant, which is typical of the circumstellar medium of
massive stars in a magnetized medium, along with the appearance of a spinning
top structure within the projected rectangle. We suggest that this mechanism
may be partially responsible for the complex morphologies observed in pulsar
wind nebulae that do not conform to the typical torus/jet or bow shock/tail
shapes observed in most cases.

###### keywords:

methods: MHD – stars: evolution – stars: massive – pulsars: general – ISM:
supernova remnants.

††pubyear: 2022††pagerange: On the plerionic rectangular supernova remnants of
static progenitors–On the plerionic rectangular supernova remnants of static
progenitors

##  1 Introduction

Pulsar wind nebulae (PWN) are intriguing astrophysical phenomena powered by
pulsars, rotating neutron stars with high magnetic field. PWN emerge from the
intricate interplay between the pulsar’s powerful wind and its surrounding
environment, a phenomenon comprehensively elucidated in the review articles by
Gaensler & Slane (2006); Kargaltsev et al. (2017); Olmi & Bucciantini (2023).
These enigmatic nebulae manifest across the electromagnetic spectrum,
predominantly characterized by synchrotron radiation and inverse Compton
scattering, complemented by distinctive emission lines (Reynolds et al.,
2017). Their presence is evident from radio observations, exemplified by PWN
G54.1+0.3 (Driessen et al., 2018), to X-ray emissions typified by G11.2-0.3
(Borkowski et al., 2016), and extends all the way to gamma-rays (Aharonian et
al., 2006; Abdalla et al., 2021), reaching PeV $\gamma$-ray energies (up to
$10^{15}$ eV), as prominently demonstrated by the Crab Nebula (Lhaaso
Collaboration et al., 2021), see also Acero et al. (2023).

Pulsars represent the final evolutionary phase of massive stars that do not
directly collapse into black holes. Understanding the physics of a pulsar and
its interaction with the surrounding medium requires knowledge of various
physical processes, including high-energy phenomena, fluid dynamics, general
relativity, and nuclear physics, see, e.g., Weber (1999); Bucciantini (2011);
Steiner et al. (2005); Lasky (2015); Pnigouras & Kokkotas (2015, 2016);
Pnigouras (2019). Pulsars have very powerful magnetospheres with strong
magnetic fields on the order of kilogauss (kG), which play a crucial role in
the evolution see e.g., Mestel et al. (1985). The rotating magnetospheres
extract energy from the pulsar and generate a powerful wind, see, e.g., Pétri
(2022). The interaction of the pulsar wind with the ambient medium produces
the so-called pulsar nebulae, which can be located inside or outside the
supernova remnant of the progenitor star, depending on whether the supernova
explosion had kicked off the pulsar. Various observations have triggered
investigations into such phenomena, see, amongst others, the studies of Pavan
et al. (2016); Kargaltsev et al. (2017); de Vries & Romani (2020); Igoshev
(2020); de Vries et al. (2021).

A particular class of supernova remnants containing a pulsar exhibit a
succession of structured shocks powered by the pulsar’s magnetic wind,
producing multi-wavelength polarized non-thermal emission. Examples of such
plerions include the Crab Nebula (Hester, 2008), as well as the Geminga pulsar
residing within the Vela supernova remnant (Bock et al., 1998; Popov et al.,
2019). Additionally, one can observe youthful supernova remnants hosting both
a pulsar and a pulsar-wind nebula, such as B0540-693 (Williams et al., 2008)
and G11.2-0.3 (Borkowski et al., 2016).

The modelling of PWN has been a long-standing challenge for several reasons.
First, the physics involved in PWN is inherently complex, involving the
interaction between the pulsar’s relativistic wind and the surrounding medium.
This requires a multi-disciplinary approach. Second, the environment in which
the pulsar wind is launched is often structured, as it depends on the
supernova remnant’s properties, the progenitor star’s circumstellar medium, or
the interstellar medium in which the pulsar resides. The properties of the
surrounding medium can significantly affect the dynamics and emission of the
PWN. These factors together make the modelling of PWN a complex and multi-
faceted problem, requiring sophisticated theoretical models and numerical
simulations to understand the physics at play fully. The Crab Nebula stands
out as a prominent example of a plerion. Extensive research, both theoretical
(Kennel & Coroniti, 1984; Coroniti, 1990; Begelman & Li, 1992; Begelman, 1998)
and numerical, has been dedicated to studying PWN like the Crab Nebula. These
investigations encompass relativistic axisymmetric 2D simulations (e.g.,
Komissarov & Lyubarsky, 2003, 2004; Komissarov, 2006; Del Zanna et al., 2006;
Camus et al., 2009; Komissarov & Lyutikov, 2011; Olmi et al., 2014) as well as
relativistic 3D simulations (e.g., Mignone et al., 2013; Porth et al., 2014;
Olmi et al., 2016).

Table 1:  List of models in this study. All simulations assume a rotating static massive star of mass $M_{\star}$ at solar metallicity, in a medium of number density $n_{\rm ISM}$ and organised magnetic field strength $B_{\rm ISM}$. The initial rotation rate of the central massive star is $\Omega_{\star}/\Omega_{\rm K}=0.1$.  ${\rm{Model}}$ |  $M_{\star}$ ($\rm M_{\odot}$) |  $n_{\rm ISM}$ ($\rm cm^{-1}$) |  $B_{\rm ISM}$ ($\mu\,\rm G$) | Description
---|---|---|---|---
Run-35-HD-0.79-PWN | $35$ | $0.79$ | 0 |  $35$$-$$\rm M_{\odot}$$-$induced PWN in unmagnetised warm ISM
Run-35-MHD-0.79-PWN | $35$ | $0.79$ | 7 |  $35$$-$$\rm M_{\odot}$$-$induced PWN in magnetised warm ISM

Moreover, the Crab Nebula is pivotal in advancing our understanding of pulsar
physics and their interactions with supernova remnants. Notably, the dynamics
and morphology of pulsar wind nebulae experience significant transformations
as they expand within a supernova remnant. This influence can be even more
pronounced when the pulsar receives a kick during a supernova explosion, as
observed in the case of PWN CTB 87 (Matheson et al., 2013). Extensive studies
have been conducted on the interaction between PWN and supernova remnants,
focusing on no-moving pulsars in 1D and 2D scenarios (van der Swaluw et al.,
2001; Blondin & Chevalier, 2017), including a mock complex surrounding the
neutron star (van der Swaluw, 2003; Blondin et al., 2001).

These studies have also been extended to moving pulsars inside the supernova
ejecta, revealing the development of strong bow-shocks between the pulsar wind
and supernova remnant (van der Swaluw et al., 2003; van der Swaluw et al.,
2004; Temim et al., 2017; Kolb et al., 2017; Temim et al., 2022) and the study
was extended using relativistic MHD (Bucciantini et al., 2004). In some cases,
the strong interaction between the PWN and the reverse shock of the supernova
remnant can result in compression. This interaction phase is referred to as
reverberation, during and after which the morphology, spectrum, and dynamics
of PWN could undergo significant changes, see the recent studies by Torres &
Lin (2018); Bandiera et al. (2020, 2023a, 2023b).

In later phases, as the moving pul sar leaves the supernova remnant and begins
to interact with the interstellar medium, a bow shock nebula forms around the
runaway pulsar. This intriguing phenomenon has been extensively studied
analytically (Bucciantini & Bandiera, 2001) and through numerical simulations
in two dimensions with relativistic considerations (Bucciantini, 2018).
Furthermore, research into this phenomenon has delved into three-dimensional
simulations, encompassing non-relativistic pulsar winds (Toropina et al.,
2019) and relativistic pulsar winds (Barkov et al., 2019a, b, 2020; Olmi &
Bucciantini, 2019). Despite the high complexity of these simulations and the
numerous questions they leave unanswered (Olmi & Bucciantini, 2023), all of
these studies still neglect the effects of the circumstellar medium of the
defunct star in which the supernova remnant and the PWN expand during the
initial phases.

The circumstellar medium of a defunct star is formed through the interaction
between the star’s wind and luminosity with the surrounding interstellar
medium (ISM). The shape and properties of the circumstellar medium depend on
various factors, including the evolution of the star, such as its mass, age,
and stage of evolution, as well as the characteristics of the surrounding ISM,
such as density, temperature, and magnetic field van Marle et al. (2015a);
Meyer et al. (2022a). In the context of massive stars, the circumstellar
medium undergoes successive structural changes. During its early life, it
forms an accretion disc (Liu et al., 2020; Meyer et al., 2022b; Elbakyan et
al., 2023; Burns et al., 2023). In its main sequence phase, it expands into a
wind bubble (Weaver et al., 1977; Gull & Sofia, 1979; Wilkin, 1996a). Later
on, it evolves into expanding shells (Stock & Barlow, 2010; Cox et al., 2012;
Decin, 2012; Decin et al., 2012). If a supernova explosion occurs, it leaves
behind an expanding remnant shell (Aschenbach & Leahy, 1999; Yusef-Zadeh et
al., 2003; Katsuda et al., 2018; Arias et al., 2019; Chiotellis et al., 2019;
Derlopa et al., 2019). Once the pulsar emits a relativistic and powerful wind,
it initially interacts with the surrounding supernova ejecta (Cox et al.,
1999; Sun et al., 1999; Crawford et al., 2001; Olmi & Bucciantini, 2023). As
the PWN passes through the supernova ejecta, it subsequently interacts with
the circumstellar medium of the defunct star. The distribution of ejecta,
stellar wind, and ISM gas acts as a matrix that channels the expansion of the
pulsar wind (Kolb et al., 2017; Temim et al., 2022). This is particularly
important when the supernova progenitor is a runaway star, as the bow shock
created by its surrounding stellar wind can influence the subsequent evolution
and emission of the supernova ejecta and the PWN (Meyer & Meliani, 2022).

This study aims to investigate how a magnetised ambient medium influences the
dynamics, morphologies, and emission properties of PWN with static massive
star progenitors. The multi-dimensional magnetohydrodynamical (MHD)
simulations conducted by van Marle et al. (2015b) have revealed that the
circumstellar medium of high-mass stars is significantly influenced by the
organized magnetic field of its ambient medium. This finding has profound
implications on the understanding of stellar wind bubbles around massive
stars, as previously studied by Freyer et al. (2003); Dwarkadas (2005); Freyer
et al. (2006); Dwarkadas (2007). The presence of a magnetic field can cause
expanding stellar wind bubbles to become elongated and adopt an oblong
morphology along the direction of the magnetic field lines. Our previous work
(Meyer et al., 2022a) has shown that such asymmetric pre-supernova
environments can result in a peculiar reflection of the supernova shock wave,
forming rectangular-shaped remnants like Puppis A. In this study, we further
investigate the effects of the reflection of the supernova blastwave in
asymmetric, magnetized wind bubbles that are generated by a static, rotating
star in the warm phase of the Galactic plane and how it may impact the
evolution of plerionic pulsar wind nebulae.

The paper is structured as follows. In Section 2, we present the modelling
methods used in this study. This includes the description of the numerical
simulations of pulsar wind nebulae, which are detailed in Section 3. We then
discuss the outcomes of our study in Section 4, and present our conclusions in
Section 5.

##  2 Method

In this section, we will provide a comprehensive review of the numerical setup
used in this study to generate models of PWN from static massive stars. We
will summarize the initial conditions, including both the initial and boundary
conditions, in the following paragraphs. Additionally, we will describe the
numerical methods employed in the simulations.

###  2.1 Initial conditions and boundaries conditions

This paper presents models that simulate the interaction between a star’s wind
and ejecta at all phases of its evolution with the warm ISM in the Milky Way
galaxy. The total number density of the ISM is taken to be $n_{\rm
ISM}=0.79\,\rm cm^{-3}$, while the magnetic field of the ISM is uniform and
structured, with a strength of $B_{\rm ISM}=7\,{\rm\mu G}$. In these models,
we assume that the ionized gas has a temperature of $8000\,\rm K$ (table 1).
The ambient medium is in equilibrium between the photoheating provided by the
reionizing gas around the star, as described in Osterbrock & Bochkarev (1989)
and Hummer (1994), and the radiative losses from optically-thin cooling
processes, as outlined in Wolfire et al. (2003). The cooling law used in this
study is based on the work of Wiersma et al. (2009), which is suitable for a
solar metallicity environment (Asplund et al., 2009). The cooling law accounts
for hydrogen and helium as the primary coolants at temperatures $T<10^{6}\,\rm
K$, and various metals’ emission lines at temperatures $T\geq 10^{6}\,\rm K$.
The cooling curve is further enhanced with [Oiii] $\lambda,5007$ collisionally
excited forbidden lines, as described in Asplund et al. (2009) and Henney et
al. (2009).

This paper presents a model that captures the evolution of the circumstellar
medium surrounding a static massive star with an initial mass of $35\,\rm
M_{\odot}$ at the zero-age main sequence. The star is considered to be
rotating with an angular velocity ratio of $\Omega_{\star}/\Omega_{\rm
K}=0.1$, where $\Omega_{\star}$ represents the star’s initial angular
frequency and $\Omega_{\rm K}$ is its equatorial Keplerian angular velocity.
Consequently, the equatorial velocity of the star can be expressed as,

$v_{\rm rot}(t)=\Omega_{\star}(t)R_{\star}(t).$ (1)

Here, $R_{\star}(t)$ denotes the stellar radius, and the time-dependence
signifies the variation in surface properties throughout the star’s entire
lifespan. The model tracks the complete evolution of the circumstellar medium
surrounding the static star, ranging from the onset of the zero-age main
sequence to the pre-supernova phase. This comprehensive approach encompasses
various stages, including the main sequence, red supergiant, and final Wolf-
Rayet phase.

Regarding the stellar wind throughout the evolution phase, we assume that the
stellar wind maintains spherical symmetry throughout the entire lifespan of
the supernova progenitor, with the axis of rotation of the rotating star
aligned with the axis of symmetry of the domain. To determine the wind’s
characteristics, we use the one-dimensional stellar evolution model provided
by the geneva library, as described in Ekström et al.
(2012)111https://www.unige.ch/sciences/astro/evolution/en/database/syclist/.
Specifically, we extract the mass-loss rate $\dot{M}(t)$ and the effective
temperature $T_{\rm eff}(t)$ of the star at each stage of evolution from this
database, and derive from it the wind density,

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

In this equation, $r$ represents the radial distance from the star, and
$\dot{M}(t)$ corresponds to the mass loss rate of the star at time $t$.

The terminal velocity of the stellar wind, denoted as $v_{\rm w}(t)$, is
calculated based on the escape velocity $v_{\rm esc}(t)$. The escape velocity
depends on the star’s effective temperature $T_{\rm eff}$ and is determined
using the conversion law,

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

where $G$ represents the gravitational constant, and $\beta(T)$ is a the
normalisation factor introduced by Eldridge et al. (2006).

We adopt the time-dependent evolution of the surface magnetic field
$B_{\star}$ of the supernova progenitor as derived in Meyer et al. (2023)
where the magnetic field strength at the surface of the star are scaled to
that of the Sun, as described in Scherer et al. (2020); Herbst et al. (2020);
Baalmann et al. (2020, 2021); Meyer et al. (2021b). Specifically, we assume a
magnetic field strength at star surface of $B_{\star}=500\,\rm G$ during the
main-sequence phase (Fossati et al., 2015; Castro et al., 2015; Przybilla et
al., 2016; Castro et al., 2017), to a Betelgeuse-like field of
$B_{\star}=0.2\,\rm G$ for the red supergiant phase (Vlemmings et al., 2002,
2005; Kervella et al., 2018) and $B_{\star}=100\,\rm G$ during the Wolf-Rayet
phase (Meyer, 2021). Concerning the stellar magnetic field structure, we
utilize a Parker spiral made of a radial component,

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

and a toroidal component,

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

respectively, with,

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

being the latitude-dependent surface velocity of the rotating massive star
(Parker, 1958; Weber & Davis, 1967; Pogorelov & Semenov, 1997; Pogorelov &
Matsuda, 2000; Chevalier & Luo, 1994; Rozyczka & Franco, 1996).

At the end of a star’s evolution, it enters the supernova phase, during which
we model the expanding supernova ejecta as a spherically symmetric
distribution within a radius of $r_{\rm max}$. The ejecta has a total energy
of $E_{\rm SN}=10^{51}$ erg and a mass of $M_{\rm SN}$, which takes into
account the star’s mass loss throughout its entire evolution until the
immediate pre-supernova time $t_{\rm psn}$, as well as the mass $M_{\rm NS}$
of the neutron star that forms at the centre. Specifically, we set,

$M_{\rm sn}=M_{\star}-\int_{0}^{t_{\rm psn}}\dot{M}(t)dt-M_{\rm
NS}=10.12M_{\odot},$ (7)

with $t_{\rm psn}$ and $M_{\rm NS}=1.4\,\rm M_{\odot}$ (Das et al., 2022).

In our study, we adopt a density and velocity profile for the freely expanding
supernova ejecta based on the work by Truelove & McKee (1999). This profile
consists of two distinct regions (Bandiera et al., 2021). The first region is
a uniform density core extending from 0 to $r_{\rm core}$, where $r_{\rm
core}$ represents the core radius. In this region, the density decreases with
time following a power-law relationship of $t^{-3}$, where $t$ denotes the
time after the explosion. The second region is the outer edge, extending from
$r_{\rm core}$ to $r_{\rm max}$, where $r_{\rm max}$ corresponds to the
maximum radius. In this region, the density decreases steeply with radius,
following a power-law relationship of $\rho\propto r^{-n}$, with the exponent
$n$ set to 11. Additionally, the density in the outer edge region decreases
with time as $t^{-(3+n)}$. These density profiles can be expressed as follows:

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

and,

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

respectively. These density profiles are commonly used for core-collapse
supernovae (Chevalier, 1982).

For the velocity, we utilize a homologous radial profile for the supernova
ejecta, given by $v=t/r$, across all regions from 0 to $r_{\rm max}$ at time
$t_{\rm max}$. The characteristics of the supernova ejecta profile are
computed following the methodology outlined in Truelove & McKee (1999) and
Whalen et al. (2008).

The velocity at the core radius, denoted as $v_{\rm core}(r_{\rm core})$, is
determined as,

$v_{\rm core}(r_{\rm core})=\left(\frac{10(n-5)E_{\rm sn}}{3(n-3)M_{\rm
sn}}\right)^{1/2},$ (10)

where $E_{\rm sn}$ represents the total energy of the supernova ejecta and
$M_{\rm sn}$ is the total mass of the ejecta. This equation ensures
conservation of both mass and energy in the supernova ejecta. The maximum
speed, denoted as $v_{\rm max}$, is set to:

$v_{\rm max}=\frac{r_{\rm max}}{t_{\rm max}}=3\times 10^{4}\,\rm km\,\rm
s^{-1},$ (11)

This choice of $v_{\rm max}$ maintains the conservation of total mass and
energy in the supernova ejecta (van Veelen et al., 2009).

As the supernova ejecta are expelled, we set a radial pulsar’s wind that
emanates from the centre, as described by Meyer et al. (2022a). This wind has
a total power that is assumed to evolve over time, $t$, according to the
following equation:

$\dot{E}_{\rm psw}=\dot{E}_{\rm
psw,0}\left(1+\frac{t}{\tau_{0}}\right)^{-\frac{n+1}{n-1}}.$ (12)

with $\tau_{\rm o}$ is the initial spin-down of the pulsar, defined as,

$\tau_{\rm o}=\frac{P_{\rm o}}{(n-1)\dot{P}{\rm o}},$ (13)

where the initial spin period of the pulsar is set to $P_{\rm o}=0.3\,\rm s$,
and its time derivative is set to $\dot{P}_{\rm o}=10^{-17}\,\rm s\,\rm
s^{-1}$. The braking index is assumed to be $n=3$, which corresponds to
magnetic dipole spin-down, as outlined in Pacini (1967).

Furthermore, we assume that the pulsar’s wind maintains a constant speed of
$v_{\rm psw}=10^{-2}c$, where $c$ denotes the speed of light in vacuum. It’s
important to acknowledge that this speed is significantly lower than the
realistic pulsar wind speeds, which can reach $c$, corresponding to a Lorentz
factor of $10^{6}$ (as demonstrated in Kennel & Coroniti, 1984). This decision
to employ a reduced pulsar wind speed can lead to noticeable alterations in
the properties of its termination shock. These changes encompass compression
rates, speeds, and, subsequently, shock positions and influence the
development of associated instabilities. It is crucial to emphasize that our
paper’s primary objective is to replicate the overall evolution of the PWN
accurately. This evolution is predominantly governed by the wind’s momentum
flux (Wilkin, 1996b). In terms of magnetization, we have opted for a low value
of $\sigma=10^{-3}$ in this study, a choice in line with descriptions found in
Rees & Gunn (1974), Kennel & Coroniti (1984), Slane (2017), Begelman & Li
(1992) and Torres et al. (2014). This magnetization value implies that a
significant portion of the magnetic field is converted into kinetic energy. It
is worth noting that recent multi-dimensional simulations have demonstrated
that larger magnetization values, such as $\sigma=0.01$ in 2D (e.g.,
Komissarov & Lyubarsky, 2003, 2004; Del Zanna et al., 2004; Del Zanna et al.,
2006) and even $\sigma>1$ in 3D (Porth et al., 2014; Barkov et al., 2019a),
can accurately reproduce the features of termination shocks of PWN.

However, it is essential to recognize that the value of pulsar wind
magnetization remains a topic of debate, as it significantly influences PWN
termination shock strength and, consequently, particle acceleration. Moreover,
the magnetization of the equatorial wind zone may decrease, leading to lower
magnetization values due to the annihilation of equatorial wind magnetic
stripes (Coroniti, 2017). By selecting such a low magnetization, as in
Bucciantini et al. (2004), the Pulsar Wind Nebula tends to expand more in the
equatorial plane, resulting in a stronger termination shock.

Furthermore, Komissarov and Lyubarsky in 2003-2004 and Del Zanna et al. 2004,
showed that the properties of the inner nebula can only be recovered with 2D
simulations if the injected magnetization is larger than 0.01.

The magnetic field is assumed to have only a toroidal component. The total
kinetic energy, magnetic field strength, and kinetic energy are functions of
the radial distance $r$ and polar angle $\theta$, as described in Komissarov &
Lyubarsky (2004).

$\displaystyle\dot{E}_{\rm total}=\frac{\dot{E}_{\rm
psw}}{r^{2}}\left(\sin^{2}{\theta}+1/\sigma\right),$ (14) $\displaystyle
B=\sqrt{\frac{4\pi}{c}}\frac{1}{r}\sin{\theta}\left(1-2\theta/\pi\right),$
(15) $\displaystyle\dot{E}_{\rm kinetic}=\dot{E}_{\rm
total}-\frac{B^{2}}{4\pi}\,c\,.$ (16)

Our choice of a spherically symmetric supernova explosion allows us to assume
that the neutron star is at rest at the location of the explosion and neglects
any potential kick velocity resulting from asymmetries in the explosion.

###  2.2 Numerical methods

To investigate the evolution of the PWN within the circumstellar medium of its
static progenitor star that is surrounded by a magnetized external medium, we
follow the strategy we used in Meyer et al. (2015); Meyer et al. (2020) and
that we extended after to PWN in Meyer & Meliani (2022). The magneto-
hydrodynamical simulations are conducted in a 2.5-dimensional, axisymmetric
cylindrical coordinate system. The simulation box extends over the range
$[O;R_{\rm max}]\times[z_{\rm min};z_{\rm max}]$ and is discretized using a
uniform grid of $N_{\rm R}\times N_{\rm z}$ cells. Consequently, the spatial
resolution is consistent along both directions, with each grid cell having a
size of $\Delta=R_{\rm max}/N_{\rm R}$. We employ two different spatial
resolutions throughout the evolutionary process. During the progenitor star
wind phases, the circumstellar medium is resolved using a grid resolution of
$N_{\rm R}=2000$ and $N_{\rm z}=4000$ cells. The stellar wind is implemented
as an internal boundary condition within a sphere centred at the origin of the
computational domain, with a radius of $20\Delta$, following the standard
procedure outlined in Comerón & Kaper (1998).

At the immediate pre-supernova stage, we remap the solution for the
circumstellar medium onto a finer grid with $N_{\rm R}=3000$ and $N_{\rm
z}=6000$ cells. The supernova ejecta is confined within a central sphere of
radius $r_{\rm max}$, as described in section 2.1. Simultaneously, the pulsar
wind is imposed within a sphere of radius $r_{\rm ns\,wind}=20\Delta$, also
detailed in section 2.1. Due to our choice of an asymmetric coordinate system,
we are compelled to align the directions of the pulsar spin axis and the
symmetry axis of the computational domain to be the same.

In this paper, we study the evolution of the circumstellar medium influenced
by the magnetized wind emitted by a massive star with a mass of $35\,\rm
M_{\odot}$ in two distinct types of external medium: the magnetized and
unmagnetized warm phases of the Galactic plane in the Milky Way. We refer to
these models as Run-35-HD-0.79-PWN and Run-35-MHD-0.79-PWN. In the magnetised
external medium case, the adopted strength of the background magnetic field is
set to that measured in the spiral arms of the Galaxy, with an average
strength of $B_{\rm ISM}=7\,\mu\rm G$ (see Draine, 2011). The main parameters
utilized in the two cases investigated in this paper are provided in Table 1.
For a more comprehensive description of the model and the implemented
strategy, please refer to Meyer et al. (2023) and Meyer & Meliani (2022),
where detailed explanations can be found.

Figure 1:  Number density fields in our magneto-hydrodynamical simulation of
the pulsar wind nebula forming in the supernova remnant of a static $35\,\rm
M_{\odot}$ star rotating with $\Omega_{\star}/\Omega_{\rm K}=0.1$ in an
unmagnetised (left) and magnetised (right) ISM. The red contours highlight the
region with a $50\%$ contribution of pulsar wind material, i.e. the contact
discontinuity. The streamlines in the right-hand side of panels b,d,f mark the
ISM magnetic field lines.

The numerical simulations are conducted using the pluto code (Mignone et al.,
2007, 2012; Vaidya et al., 2018)222http://plutocode.ph.unito.it/ and we solve
the following set of equations,

$\frac{\partial\rho}{\partial
t}+\mn@boldsymbol{\nabla}\cdot\big{(}\rho\mn@boldsymbol{v})=0,$ (17)
$\frac{\partial\mn@boldsymbol{m}}{\partial
t}+\mn@boldsymbol{\nabla}\cdot\Big{(}\mn@boldsymbol{m}\otimes\mn@boldsymbol{v}-\mn@boldsymbol{B}\otimes\mn@boldsymbol{B}+\mn@boldsymbol{\hat{I}}p_{\rm
t}\Big{)}\mn@boldsymbol{0},$ (18) $\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),$
(19)

and

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

with the gas density $\rho$, velocity $v$, momentum
$\mn@boldsymbol{m}=\rho\mn@boldsymbol{v}$ and magnetic field
$\mn@boldsymbol{B}$, as well as the the total pressure $p_{t}$ and the energy
of the gas

$E=\frac{p}{(\gamma-1)}+\frac{\mn@boldsymbol{m}\cdot\mn@boldsymbol{m}}{2\rho}+\frac{\mn<EMAIL_ADDRESS>(21)

The sound speed of the medium reads

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

where the adiabatic index is $\gamma=5/3$. Last, radiative cooling by
optically-thin processes and photo-heating are included into the equations via
the term $\Phi(T,\rho)$, with the gas temperature $T$, accounting for the
prescriptions of Meyer et al. (2014). Regarding the cooling/heating processes
of the gas, we assume the gas to be optically thin throughout the entire
progenitor’s life. After this point, with the launch of the pulsar wind, the
cooling and heating mechanisms are disabled. We employ a Godunov-type
numerical scheme with the Harten-Lax-van Leer approximate Riemann solver (HLL)
and utilize the 8-waves magnetic field formulation (Powell, 1997). For time
integration, a third-order Runge-Kutta scheme is employed, with the time step
controlled by the Courant-Friedrichs-Lewy (CFL) number.

The numerical simulations are performed at the North-German Supercomputing
Alliance (HLRN333https://www.hlrn.de/) using the LISE cluster in Berlin, which
is equipped with Cray XC40/30 processors.

##  3 Results

In this section, we will analyze the results of the evolution of the PWN
within the supernova remnant and circumstellar medium of the progenitor star
in both the unmagnetized and magnetized cases. Our focus will be on
investigating the influence of the magnetic field of the progenitor star and
the external medium on the shape and dynamics of the PWN.

###  3.1 Model with unmagnetised ISM

In Figure 1, the density contour is shown for the unmagnetized case
Run-35-HD-0.79-PWN (left panels) and the magnetized case Run-35-MHD-0.79-PWN
(right panels) at different evolution times, from to top to the bottom. The
density contour is represented in logarithmic scale in $\rm cm^{-3}$ units. In
both cases, the red contour marks indicate the region of the plerionic
supernova remnants where the contribution of the pulsar wind reaches $50$
times of the number density.

In Figure 1a (top-left), we present the pre-supernova circumstellar medium. At
this stage, it forms a large-scale quasi-spherical stellar bubble (Weaver et
al., 1977), and its spherical forward shock extends to distances of
approximately $90\,\rm pc$. Throughout the star’s evolution, the stellar wind
interacts strongly with the ambient medium. Each phase of evolution
contributes to the formation of successive shock structures, which appear in
order from the farthest to the nearest region to the star. The thick and dense
shell located farthest from the star, with a radial extent of $\geq 50\,\rm
pc$, is the result of the interaction between the stellar wind and the ISM,
and it occurs mainly during the main-sequence phase (Freyer et al., 2003;
Dwarkadas, 2005; Freyer et al., 2006; Dwarkadas, 2007). In the central region,
within a radius of less than $20\,\rm pc$, a low-density cavity is formed due
to the continuous outflow of the free stellar wind during the Wolf-Rayet
phase. This cavity is surrounded by successive dense shells resulting from the
interactions between the Wolf-Rayet wind and the slower wind from the
preceding red-giant phase. The first shell, extending to approximately
$35\,\rm pc$, is dense and exhibits unstable behaviour. Subsequently, a
second, less dense shell is formed due to the interaction between the red-
giant wind and the main-sequence wind. Additionally, the main-sequence wind
interacts with the surrounding ambient medium, forming an external dense shell
that is limited by the contact discontinuity surface.

It is worth noting that the contact discontinuity, which marks the interface
between the wind and the ISM, exhibits a slightly aspherical morphology,
particularly in the region close to the symmetry axis. This aspherical shape
is influenced by the presence of the magnetic field and the star’s rotation.
The variations between the bubbles depicted in Fig. 1a and Fig. 1 of Meyer et
al. (2022a) highlight this effect. Furthermore, the grid’s proximity to the
near-polar axis amplifies this asymmetry. Moving on to Fig. 1c (middle-left),
we can observe the supernova remnant at $25\,\rm kyr$ after the explosion. The
expanding shock wave from the supernova remnant propagates outward, sweeping
up and pushing away all the previously formed dense shells associated with the
successive stellar winds. As the shock wave reaches the contact discontinuity
surface between the main-sequence stellar wind and the ISM, it interacts with
this surface, causing reflection, as described in Meyer et al. (2015, 2021a);
Meyer & Meliani (2022). This interaction and reverberation contribute to the
observed structure and morphology of the supernova remnant.

After the supernova explosion, a pulsar wind with high initial mechanical
luminosity, $\dot{E}_{\rm psr,0}=10^{38}\rm ergs^{-1}$, is launched. However,
this luminosity decreases over time according to Eq. 21 of Pacini (1967). This
pulsar wind interacts with the dense supernova ejecta (van der Swaluw et al.,
2004), resulting in the formation of a complex structure as described in Meyer
& Meliani (2022) for a runaway progenitor star with a zero-age main-mass of
$20\,\rm M_{\odot}$. Within this structure, the central region of the plerion
is occupied by the freely-expanding pulsar wind. Surrounding the central
region, a shell of shocked pulsar wind is formed, resulting from the
interaction of the pulsar wind with the expanding supernova remnant. A pulsar
wind termination shock is formed at the interface between the unperturbed
pulsar wind and the shocked pulsar wind. The outermost region of the pulsar
wind nebula behind the termination shock contains the contact discontinuity.
This contact discontinuity marks the interface between the supernova ejecta
and the shocked pulsar wind (depicted by the red contour in Fig. 1). Beyond
the contact discontinuity, a transmitted pulsar wind forward shock propagates
through the still unshocked supernova ejecta and further travels into the
surrounding medium.

The pulsar wind contact discontinuity undergoes expansion to larger radii due
to the fast rotation of the magnetized neutron star. This expansion leads to
the characteristic shape with an equatorial torus and an elongated polar jet,
as found by Komissarov & Lyubarsky (2004); Del Zanna et al. (2006); Porth et
al. (2014); Olmi et al. (2016). However, it’s important to note that due to
limitations in the numerical scheme applied to the 2D symmetry axis, the jet
along the polar axis may appear more elongated than it would in a full 3D
simulation. Nevertheless, despite these limitations, the general behaviors of
the PWN remains accurate. This shape can be observed at a later time,
specifically $45\,\rm kyr$ after the explosion, as shown in Fig. 1e. As the
contact discontinuity surface expands, it encounters Rayleigh-Taylor
instabilities due to the significant differences in density and velocity
between the pulsar wind and the supernova ejecta. These instabilities are
further amplified by the reverberation of the reverse shock from the supernova
ejecta, as illustrated in Fig. 1e.

###  3.2 Model with magnetized ISM

During the main-sequence phase of a massive star, the influence of the ISM
magnetic field becomes particularly significant. During this phase, the
interaction between the stellar wind and the magnetized ISM carves out a
large-scale circumstellar wind bubble. This wind bubble plays a crucial role
in shaping the propagation of the supernova forward shock. Additionally, the
wind bubble’s presence influences the pulsar wind’s dynamics, further
highlighting the interplay between the stellar wind, the ISM magnetic field,
and the subsequent evolution of the system. We will describe it in detail in
the following. In Fig. 1b (top-right), we can observe the circumstellar medium
surrounding the massive star in the presence of a magnetized ISM, as
represented in the model Run-35-MHD-0.79-PWN. The black arrows indicate the
magnetic field lines of the ISM, which are initially aligned with the polar
axis. The overall structure of the circumstellar medium in the presence of the
magnetized ISM remains similar to the unmagnetized model (Run-35-HD-0.79-PWN,
Fig. 1a). However, the morphology of the shocked shells within the low-density
cavity, up to the contact discontinuity between the shocked stellar wind and
the shocked ISM, appears to be more elongated along the polar axis due to the
influence of the ISM magnetic field.

Indeed, as the expanding stellar bubble interacts with the magnetized ISM, it
compresses the magnetic field lines, increasing magnetic pressure and tension
along the polar axis. This phenomenon has been extensively studied and
described in detail in van Marle et al. (2015b). During the last evolution
phase, when the Wolf-Rayet wind material reaches the main-sequence termination
shock, it undergoes reflection near the equator. This anisotropic reflection
causes a change in the direction of propagation of the shocked material,
resulting in the loss of the initially spherical shape of the shocked shell
from the Wolf-Rayet wind. The interaction with the magnetized ISM further
influences the shape and dynamics of the shocked shell, leading to the
observed rectangular morphology of the resulting supernova ejecta.
Furthermore, as the expanding supernova blast wave propagates within the
elongated cavity (as shown in the left panel of Fig. 1), it interacts with the
reflected dense shells resulting from the Wolf-Rayet wind and the elongated
contact discontinuity. These interactions lead to anisotropic reverberation at
the contact discontinuity of the supernova ejecta. As a result, the shape of
the supernova ejecta becomes rectangular, reflecting the influence of the
asymmetric interactions with the elongated structure induced by the magnetized
circumstellar medium. This mechanism is specifically described within the
context of the remnant Puppis A in Meyer et al. (2022a).

In Fig. 1d and f, the influence of the magnetized ISM on the shaping of the
PWN can be observed. The ISM magnetic field, which plays a significant role in
determining the morphology of the circumstellar medium and supernova
blastwave, also affects the confinement and shape of the pulsar wind. Under
the influence of the ISM magnetic field, the reflected and the supernova
blastwave adopts a rectangular morphology along the direction perpendicular to
the magnetic field. This happens because the ram pressure of the supernova
ejecta is directed towards the polar axis, causing compression and confinement
of the pulsar wind in that direction. In contrast, in the direction parallel
to the magnetic field, the pressure exerted by the supernova ejecta is lower,
resulting in a more extended shape of the PWN. This interplay between the
magnetic field of the ISM, the reverse shock, and the pulsar wind contributes
to the complex and asymmetric morphology observed in the PWN, as depicted in
Fig. 1d and f.

Indeed, the presence of a magnetized ISM influences the propagation of the
PWN, resulting in distinct behavior compared to an unmagnetized ISM. In the
magnetized ISM model (Run-35-MHD-0.79-PWN), the expansion of the PWN is less
pronounced in the equatorial plane compared to the hydrodynamical simulation
(Run-35-HD-0.79-PWN), as illustrated in Fig.1c and d. As time progresses, at a
later evolution time of $45\,\rm kyr$ as depicted in Fig.1f, the pulsar wind
continues to be channeled along the direction of the ISM’s magnetic field,
leading to the formation of a stretched PWN. The presence of the ISM magnetic
field affects the dynamics of the PWN and leads to enhanced instabilities at
the termination shock of the pulsar wind. These instabilities, which arise
from the interaction between the pulsar wind and the magnetized ISM, are more
pronounced in the magnetized ISM model (Run-35-MHD-0.79-PWN) compared to the
hydrodynamical simulation (Run-35-HD-0.79-PWN).

Our models provide compelling evidence that the morphology of the PWN inside a
subsequent supernova, when the progenitor massive static star is located in
the Galactic plane, is strongly influenced by the distribution of the magnetic
field in the ambient medium. The contrasting evolution and instabilities
observed in the magnetized and unmagnetized cases emphasize the significant
role played by the interstellar medium’s magnetic field in shaping the
dynamics and morphology of the PWN. These findings underscore the importance
of considering the magnetic field effects when studying the evolution of PWN
and their interaction with the surrounding environment.

Figure 2:  Normalised radio synchrotron emission map of the plerionic
supernova remnants with an inclination angle of $\theta_{\rm obs}=45^{\circ}$
between the observer’s line of sight and the nebula’s symmetry axis. The left-
hand panels correspond to the hydrodynamical model (Run-35-HD-0.79-PWN), and
right-hand panel to the model with magnetised ISM (Run-35-MHD-0.79-PWN). The
top figures show the remnants at time $25\,\rm kyr$ and the bottom figures
display them at time $45\,\rm kyr$.

##  4 Discussion

In this section, we will discuss the applications and limitations, of our
model. We will also examine the non-thermal characteristics of the simulated
pulsar wind nebulae and compare our findings to existing observational data.
By doing so, we aim to provide a comprehensive analysis of our model’s
strengths and weaknesses and assess its compatibility with the observed
properties of pulsar wind nebulae.

###  4.1 Model limitations

Let us first consider four aspects central to the model. First, the
simulations conducted in this study are two-dimensional, assuming axisymmetry
and not accounting for variations in the supernova progenitor or the pulsar’s
spin. While this approach offers computational efficiency and valuable
insights, it’s essential to acknowledge that a fully three-dimensional
treatment is not only important to capture the realistic properties of the ISM
but also crucial for a comprehensive understanding of the pulsar wind nebula
and the supernova remnant. A 3D model would better represent the complex
interactions of the PWN and supernova remnant with the surrounding medium,
including the realistic behaviour of magnetic fields. Moreover, the
magnetization of the pulsar wind is a fundamental parameter that plays a
significant role in the evolution of the PWN and supernova remnant. While we
have considered a weak magnetization of the pulsar wind in this study, it’s
essential to discuss its implications thoroughly. State-of-the-art simulations
in both 2D and 3D have shown that the strength and longitudinal variation of
magnetization are subjects of debate (Coroniti, 2017; Olmi et al., 2016).
Future investigations will explore the influence of higher magnetization on
the evolution of the PWN in its interaction with supernova remnant and
circumstellar medium.

Furthermore, we acknowledge that our modelling of the pulsar wind involves
simplified assumptions. A more realistic modelling approach should also
involve a better physical description of the wind properties, including its
relativistic speed and composition. Addressing these aspects will be crucial
in future research for a more comprehensive understanding of the system’s
dynamics and morphology. Another aspect to consider is the absence of pulsar
motion in the simulations. Incorporating pulsar motion would introduce
additional complexities and offer a more realistic representation of the
interaction between the pulsar wind and the surrounding medium. Furthermore,
accounting for the oblique rotation of the pulsar’s magnetic axis would allow
for a more accurate reproduction of the observed characteristics of the PWN.
These are important considerations for future research. The chosen two-
dimensional setup and static pulsar position provide valuable insights into
general behaviour and trends. However, future investigations can explore the
impact of three-dimensional effects, pulsar motion, higher magnetization, and
improved modelling of the pulsar wind to obtain a more comprehensive
characterization of the system’s dynamics and morphology.

###  4.2 Non-thermal emission

To enhance the comparison between our MHD models of the PWN embedded in an
elongated circumstellar medium and the available observational data, we
performed radiative transfer calculations to generate synthetic images that
accurately capture the non-thermal emissions, particularly synchrotron
emissions in the radio band. These calculations were specifically carried out
at the different evolution stages of the PWN that were previously discussed.
The synchrotron radio emission was calculated by considering a non-thermal
electron spectrum described by the expression,

$N(E)=KE^{-s}\propto nE^{-2},$ (23)

where here, $n$ represents the gas number density, $s$ is the spectral index
and $E$ denotes the energy of the non-thermal electrons in the post-shock
region of the advancing blast wave. The emission coefficient is given by:

$j_{\rm sync}(\nu)\propto nB_{\perp}^{(s+1)/2}\nu^{-(s-1)/2},$ (24)

being $\nu$ the observed frequency and $B_{\perp}$ the component of the
magnetic field perpendicular to the observer’s line-of-sight.

Intensity maps were obtained by performing the projection given by,

$I=\int_{\rm SNR}j_{\rm sync}\left(\theta_{\rm obs}\right)dl,$ (25)

where $\theta_{\rm obs}$ denote the inclination angle of the remnant with
respect to the sky plane. These calculations were conducted using the
radiative transfer code RADMC-3D444https://www.ita.uni-
heidelberg.de/dullemond/software/radmc-3d/, and the methodology described in
detail by Meyer et al. (2022a). Note that since the investigated numerical
simulations are non-relativistic, they do not account for the beaming effect,
and this issue will be addressed in our upcoming work.

Figure 2 illustrates the normalized emission maps representing our numerical
simulations, specifically Run-35-HD-0.79-PWN (left-hand panels) and
Run-35-MHD-0.79-PWN (right-hand panels), showcasing the non-thermal
synchrotron emissions in the radio waveband. The top panels correspond to
$25\rm kyr$, while the bottom panels depict a time of $45\,\rm kyr$. The
intensity is plotted assuming an observer angle ($\theta_{\rm obs}$) of
$45^{\circ}$, representing the angle between the plane of the sky and the
plane of symmetry of the supernova remnant. Figure 2a displays the pulsar wind
nebula at the age of $25\,\rm kyr$ within an unmagnetized ISM. As highlighted
in Meyer & Meliani (2022), no trace of the circumstellar medium is visible in
the emission maps because of the absence of the ISM magnetic field. Indeed,
the emission map focuses on the pulsar wind and its associated nebula. The
image reveals an ovoidal shape, with slightly brighter regions observed at the
polar and dimmer regions in the equatorial plane. This brightness variation
can be attributed to the toroidal component of the pulsar wind, which applies
lateral pressure on the pulsar wind material, causing it to be displaced
sideways in the equatorial plane.

At a later evolution time, with a pulsar age of $45\,\rm kyr$, the radio
synchrotron map of the PWN in an unmagnetized ISM is shown in Fig. 2c. The PWN
exhibits a jet-torus-like shape, with brighter regions observed at the polar
zones. These bright regions result from the strong interaction between the
pulsar wind and the supernova ejecta along the pulsar’s rotational axis. On
the other hand, in the equatorial plane, the strong pulsar wind, driven by the
centrifugal force and toroidal magnetic field pressure (Komissarov &
Lyubarsky, 2004), extends outward. The gas is more diluted in this region,
which explains why the equatorial plane is not the brightest region in the
hydrodynamical plerion model Run-35-HD-0.79-PWN. In the case of a magnetized
ISM, significant changes are observed in the synthetic radio image. The
corresponding image is shown at $25\,\rm kyr$ in Fig. 2b. It reveals the
presence of two bright arcs parallel to the direction of the ISM magnetic
field. These arcs, observed in our axisymmetric setup and aligned with the
pulsar’s rotation axis, are formed as a result of the interaction between the
supernova ejecta and the contact discontinuity between the stellar wind and
the magnetized ISM within the elongated cavity (Meyer et al., 2022a). The
influence of the ISM magnetic field plays a crucial role in shaping these
arcs, ultimately forming a PWN enclosed within a rectangular supernova
remnant.

Fig. 2d depicts the older remnant within a magnetized ambient medium,
showcasing characteristics of both a supernova shock wave that has interacted
with the cavity’s border and the growing pulsar wind nebula inside it. The
presence of the pulsar wind prevents the reverberation of the supernova shock
wave towards the centre of the explosion, as described in Meyer et al.
(2022a), resulting in an empty region near the rotating neutron star. The
overall morphology of the plerionic remnant still exhibits features of a
rectangularly reflected supernova shock wave, with the pulsar wind distributed
as an elongated structure. The brightest regions are observed as two polar
spots located beyond the termination shock of the pulsar wind.

###  4.3 Comparison with observations

The models presented in this study focus on the evolution of the circumstellar
medium surrounding static high-mass stellar objects that eventually undergo
supernova explosions, leaving behind a static pulsar. We aim to investigate
the formation of elongated pulsar wind nebulae, similar to those observed in
Igoshev (2020). It is important to note that these elongated PWN, where the
leptonic wind is channelled into the cavity created by the stellar wind shaped
by the organized ISM magnetic field, should not be confused with the long
tails observed behind the bow shocks of runaway pulsars (e.g., Bucciantini,
2002, 2018; De Luca et al., 2013; Barkov et al., 2019a). The class of
torus/jet-like pulsar wind nebulae, as classified in the catalogue based on
Chandra X-ray data, provides strong support for the conclusions drawn from our
model. These objects naturally exhibit both an equatorial structure and a
jet/counter-jet system, as observed in studies such as Kargaltsev & Pavlov
(2010); Kargaltsev et al. (2012) and references therein. Notable examples
include the famous Crab nebula with its twisted double jet (Mignone et al.,
2013) and the Vela supernova remnant. Magneto-hydrodynamical models have
successfully reproduced such structures without considering the stellar wind
or supernova ejecta as initial conditions, as demonstrated in Klingler et al.
(2014).

The influence of the environment on the morphology of pulsar wind tails/jets
has been demonstrated in cases such as the Geminga pulsar wind nebula, which
exhibits two curved antennae representing its jets/counter-jet that bend under
the influence of the bow shock formed due to the interaction between the fast
pulsar motion and the surrounding medium (Posselt et al., 2017). Similar
effects have been observed in the case of B0355+54 (Klingler et al., 2014). We
propose that the pre-supernova environment plays a similar role, and further
modelling efforts are highly desirable, as discussed in Meyer & Meliani
(2022). The peculiar morphology of certain pulsar wind nebulae, which cannot
be classified as either torus/jet-like objects or bow shock/tail systems, may
result from their interaction with a particularly complex surrounding medium.
This medium could be shaped by the asymmetric stellar wind during the evolved
phases of the progenitor’s pre-supernova life, which influences the forward
shock of the ejecta and causes aspherical propagation (Velázquez et al., 2023;
Villagran et al., 2023).

##  5 Conclusion

This paper presents a study on the modelling of PWN in core-collapse supernova
remnants associated with static massive stars in the warm phase of a
magnetized spiral arm of the Milky Way. By utilizing 2.5-dimensional
simulations, we demonstrate that the reflection of the supernova blast wave
against the elongated contact discontinuity between the stellar wind and
magnetise ISM of the magnetically elongated stellar wind cavity in the
progenitor’s circumstellar medium has a significant impact on the morphology
of the resulting PWN. This phenomenon might be responsible for forming
rectangular supernova remnants, such as Puppis A, as described in Meyer et al.
(2022a). The reverberation of the shock wave leads to the compression of the
pulsar wind and imposes a preferred expansion direction perpendicular to the
plane of the pulsar’s spin. As a result, the PWN within the rectangular
supernova remnant becomes elongated rather than adopting the jet-torus-like
shape typically observed in previous studies, as described by Komissarov &
Lyubarsky (2004).

The radio synchrotron emission maps of plerionic supernova remnants exhibit a
complex morphology that evolves over time. Initially, the morphology is
characterized by a young, growing, ovoidal PWN combined with the rectangular
shape produced by the interaction between the supernova ejecta and the walls
of the unshocked stellar wind cavity of the progenitor star. This interaction
gives rise to the rectangular appearance observed in Puppis A, as discussed in
Meyer et al. (2022a). As time progresses, the influence of the ISM magnetic
field becomes more prominent in shaping the remnant’s morphology. The
channelling effect of the pulsar wind into the elongated circumstellar wind
cavity of the progenitor extends along the pulsar’s rotation axis.
Instabilities at the interface between the pulsar wind and the ejecta result
in a knotty nebula, manifesting as bright spots within the plerion. The
irregular shapes observed in many pulsar wind nebulae may indicate the complex
nature of the surrounding environment, influenced by both the distribution of
material in the ambient medium and the stellar wind history of the supernova
progenitor. In this complex environment, the interaction between the supernova
ejecta and the pulsar wind gives rise to observed irregular morphologies.

## Acknowledgements

The authors acknowledge the referee for numerous constructive comments. DMAM
thanks O. Petruk for his comments on the manuscript. The authors acknowledge
the North-German Supercomputing Alliance (HLRN) for providing HPC resources
that have contributed to the research results reported in this paper. PFV
acknowledges financial support from the PAPIIT IG100422 grant. DFT has been
supported by grants PID2021-124581OB-I00 funded by
MCIN/AEI/10.13039/501100011033, 2021SGR00426 of the Generalitat de Catalunya,
by the program Unidad de Excelencia María de Maeztu CEX2020-001058-M, and by
MCIN with funding from European Union NextGeneration EU (PRTR-C17.I1).

## Data availability

This research made use of the pluto code developed at the University of Torino
by A. Mignone (http://plutocode.ph.unito.it/). and of the radmc-3d code
developed at the University of Heidelberg by C. Dullemond
(https://www.ita.uni-heidelberg.de/$\sim$dullemond/software/radmc-3d/). 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

  * Abdalla et al. (2021) Abdalla H., et al., 2021, ApJ, 917, 6
  * Acero et al. (2023) Acero F., et al., 2023, Astroparticle Physics, 150, 102850
  * Aharonian et al. (2006) Aharonian F., et al., 2006, A&A, 457, 899
  * Arias et al. (2019) Arias M., et al., 2019, A&A, 622, A6
  * Aschenbach & Leahy (1999) Aschenbach B., Leahy D. A., 1999, A&A, 341, 602
  * 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
  * Bandiera et al. (2020) Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2020, MNRAS, 499, 2051
  * Bandiera et al. (2021) Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2021, MNRAS, 508, 3194
  * Bandiera et al. (2023a) Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2023a, MNRAS,
  * Bandiera et al. (2023b) Bandiera R., Bucciantini N., Olmi B., Torres D. F., 2023b, MNRAS, 525, 2839
  * Barkov et al. (2019a) Barkov M. V., Lyutikov M., Khangulyan D., 2019a, MNRAS, 484, 4760
  * Barkov et al. (2019b) Barkov M. V., Lyutikov M., Klingler N., Bordas P., 2019b, MNRAS, 485, 2041
  * Barkov et al. (2020) Barkov M. V., Lyutikov M., Khangulyan D., 2020, MNRAS, 497, 2605
  * Begelman (1998) Begelman M. C., 1998, ApJ, 493, 291
  * Begelman & Li (1992) Begelman M. C., Li Z.-Y., 1992, ApJ, 397, 187
  * Blondin & Chevalier (2017) Blondin J. M., Chevalier R. A., 2017, ApJ, 845, 139
  * Blondin et al. (2001) Blondin J. M., Chevalier R. A., Frierson D. M., 2001, ApJ, 563, 806
  * Bock et al. (1998) Bock D. C. J., Turtle A. J., Green A. J., 1998, AJ, 116, 1886
  * Borkowski et al. (2016) Borkowski K. J., Reynolds S. P., Roberts M. S. E., 2016, ApJ, 819, 160
  * Bucciantini (2002) Bucciantini N., 2002, A&A, 387, 1066
  * Bucciantini (2011) Bucciantini N., 2011, in High-Energy Emission from Pulsars and their Systems. p. 473 (arXiv:1005.4781), doi:10.1007/978-3-642-17251-9_39
  * Bucciantini (2018) Bucciantini N., 2018, MNRAS, 478, 2074
  * Bucciantini & Bandiera (2001) Bucciantini N., Bandiera R., 2001, A&A, 375, 1032
  * Bucciantini et al. (2004) Bucciantini N., Bandiera R., Blondin J. M., Amato E., Del Zanna L., 2004, A&A, 422, 609
  * Burns et al. (2023) Burns R. A., et al., 2023, Nature Astronomy,
  * Camus et al. (2009) Camus N. F., Komissarov S. S., Bucciantini N., Hughes P. A., 2009, MNRAS, 400, 1241
  * 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 & Luo (1994) Chevalier R. A., Luo D., 1994, ApJ, 421, 225
  * Chiotellis et al. (2019) Chiotellis A., Boumis P., Derlopa S., Steffen W., 2019, arXiv e-prints, p. arXiv:1909.08947
  * Comerón & Kaper (1998) Comerón F., Kaper L., 1998, A&A, 338, 273
  * Coroniti (1990) Coroniti F. V., 1990, ApJ, 349, 538
  * Coroniti (2017) Coroniti F. V., 2017, ApJ, 850, 184
  * Cox et al. (1999) Cox D. P., Shelton R. L., Maciejewski W., Smith R. K., Plewa T., Pawl A., Różyczka M., 1999, ApJ, 524, 179
  * Cox et al. (2012) Cox N. L. J., Kerschbaum F., van Marle A. J., Decin L., Ladjal D., Mayer A., 2012, A&A, 543, C1
  * Crawford et al. (2001) Crawford F., Gaensler B. M., Kaspi V. M., Manchester R. N., Camilo F., Lyne A. G., Pivovaroff M. J., 2001, ApJ, 554, 152
  * Das et al. (2022) Das S., Brose R., Meyer D. M. A., Pohl M., Sushch I., Plotko P., 2022, A&A, 661, A128
  * De Luca et al. (2013) De Luca A., et al., 2013, ApJ, 765, L19
  * Decin (2012) Decin L., 2012, Advances in Space Research, 50, 843
  * Decin et al. (2012) Decin L. N. L. J., et al., 2012, A&A, 548, A113
  * Del Zanna et al. (2004) Del Zanna L., Amato E., Bucciantini N., 2004, A&A, 421, 1063
  * Del Zanna et al. (2006) Del Zanna L., Volpi D., Amato E., Bucciantini N., 2006, A&A, 453, 621
  * Derlopa et al. (2019) Derlopa S., Boumis P., Chiotellis A., Steffen W., Akras S., 2019, arXiv e-prints, p. arXiv:1909.06131
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium
  * Driessen et al. (2018) Driessen L. N., Domček V., Vink J., Hessels J. W. T., Arias M., Gelfand J. D., 2018, ApJ, 860, 133
  * 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
  * Elbakyan et al. (2023) Elbakyan V. G., Nayakshin S., Meyer D. M. A., Vorobyov E. I., 2023, MNRAS, 518, 791
  * Eldridge et al. (2006) Eldridge J. J., Genet F., Daigne F., Mochkovitch R., 2006, MNRAS, 367, 186
  * Fossati et al. (2015) Fossati L., et al., 2015, A&A, 574, A20
  * 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
  * Gaensler & Slane (2006) Gaensler B. M., Slane P. O., 2006, ARA&A, 44, 17
  * Gull & Sofia (1979) Gull T. R., Sofia S., 1979, ApJ, 230, 782
  * Henney et al. (2009) Henney W. J., Arthur S. J., de Colle F., Mellema G., 2009, MNRAS, 398, 157
  * Herbst et al. (2020) Herbst K., et al., 2020, ApJ, 897, L27
  * Hester (2008) Hester J. J., 2008, ARA&A, 46, 127
  * Hummer (1994) Hummer D. G., 1994, MNRAS, 268, 109
  * Igoshev (2020) Igoshev A. P., 2020, MNRAS, 494, 3663
  * Kargaltsev & Pavlov (2010) Kargaltsev O., Pavlov G. G., 2010, in Comastri A., Angelini L., Cappi M., eds, American Institute of Physics Conference Series Vol. 1248, X-ray Astronomy 2009; Present Status, Multi-Wavelength Approach and Future Perspectives. pp 25–28 (arXiv:1002.0885), doi:10.1063/1.3475228
  * Kargaltsev et al. (2012) Kargaltsev O., Schmitt B. M., Pavlov G. G., Misanovic Z., 2012, ApJ, 745, 99
  * Kargaltsev et al. (2017) Kargaltsev O., Pavlov G. G., Klingler N., Rangelov B., 2017, Journal of Plasma Physics, 83, 635830501
  * Katsuda et al. (2018) Katsuda S., Takiwaki T., Tominaga N., Moriya T. J., Nakamura K., 2018, ApJ, 863, 127
  * Kennel & Coroniti (1984) Kennel C. F., Coroniti F. V., 1984, ApJ, 283, 710
  * Kervella et al. (2018) Kervella P., et al., 2018, A&A, 609, A67
  * Klingler et al. (2014) Klingler N., Rangelov B., Kargaltsev O., Pavlov G. G., Romani R. W., Slane P. O., XVP PWN Collaboration 2014, in American Astronomical Society Meeting Abstracts #223. p. 153.28
  * Kolb et al. (2017) Kolb C., Blondin J., Slane P., Temim T., 2017, ApJ, 844, 1
  * Komissarov (2006) Komissarov S. S., 2006, MNRAS, 367, 19
  * Komissarov & Lyubarsky (2003) Komissarov S. S., Lyubarsky Y. E., 2003, MNRAS, 344, L93
  * Komissarov & Lyubarsky (2004) Komissarov S. S., Lyubarsky Y. E., 2004, MNRAS, 349, 779
  * Komissarov & Lyutikov (2011) Komissarov S. S., Lyutikov M., 2011, MNRAS, 414, 2017
  * Lasky (2015) Lasky P. D., 2015, Publ. Astron. Soc. Australia, 32, e034
  * Lhaaso Collaboration et al. (2021) Lhaaso Collaboration et al., 2021, Science, 373, 425
  * Liu et al. (2020) Liu S.-Y., Su Y.-N., Zinchenko I., Wang K.-S., Meyer D. M. A., Wang Y., Hsieh I. T., 2020, ApJ, 904, 181
  * Matheson et al. (2013) Matheson H., Safi-Harb S., Kothes R., 2013, ApJ, 774, 33
  * Mestel et al. (1985) Mestel L., Robertson J. A., Wang Y. M., Westfold K. C., 1985, MNRAS, 217, 443
  * 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. (2020) Meyer D. M. A., Petrov M., Pohl M., 2020, MNRAS, 493, 3548
  * Meyer et al. (2021a) Meyer D. M. A., Pohl M., Petrov M., Oskinova L., 2021a, MNRAS, 502, 5340
  * Meyer et al. (2021b) Meyer D. M. A., Mignone A., Petrov M., Scherer K., Velázquez P. F., Boumis P., 2021b, MNRAS, 506, 5170
  * Meyer et al. (2022a) Meyer D. M. A., et al., 2022a, MNRAS, 515, 594
  * Meyer et al. (2022b) Meyer D. M. A., Vorobyov E. I., Elbakyan V. G., Kraus S., Liu S. Y., Nayakshin S., Sobolev A. M., 2022b, MNRAS, 517, 4795
  * Meyer et al. (2023) Meyer D. M.-A., Pohl M., Petrov M., Egberts K., 2023, arXiv e-prints, p. arXiv:2303.12579
  * 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
  * Mignone et al. (2013) Mignone A., Striani E., Tavani M., Ferrari A., 2013, MNRAS, 436, 1102
  * Olmi & Bucciantini (2019) Olmi B., Bucciantini N., 2019, MNRAS, 488, 5690
  * Olmi & Bucciantini (2023) Olmi B., Bucciantini N., 2023, Publ. Astron. Soc. Australia, 40, e007
  * Olmi et al. (2014) Olmi B., Del Zanna L., Amato E., Bandiera R., Bucciantini N., 2014, MNRAS, 438, 1518
  * Olmi et al. (2016) Olmi B., Del Zanna L., Amato E., Bucciantini N., Mignone A., 2016, Journal of Plasma Physics, 82, 635820601
  * Osterbrock & Bochkarev (1989) Osterbrock D. E., Bochkarev N. G., 1989, Soviet Ast., 33, 694
  * Pacini (1967) Pacini F., 1967, Nature, 216, 567
  * Parker (1958) Parker E. N., 1958, ApJ, 128, 664
  * Pavan et al. (2016) Pavan L., et al., 2016, A&A, 591, A91
  * Pétri (2022) Pétri J., 2022, A&A, 666, A5
  * Pnigouras (2019) Pnigouras P., 2019, Phys. Rev. D, 100, 063016
  * Pnigouras & Kokkotas (2015) Pnigouras P., Kokkotas K. D., 2015, Phys. Rev. D, 92, 084018
  * Pnigouras & Kokkotas (2016) Pnigouras P., Kokkotas K. D., 2016, Phys. Rev. D, 94, 024053
  * 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
  * Popov et al. (2019) Popov M. V., Andrianov A. S., Burgin M. S., Zuga V. A., Rudnitskii A. G., Smirnova T. V., Soglasnov V. A., Fadeev E. N., 2019, Astronomy Reports, 63, 391
  * Porth et al. (2014) Porth O., Komissarov S. S., Keppens R., 2014, MNRAS, 438, 278
  * Posselt et al. (2017) Posselt B., et al., 2017, ApJ, 835, 66
  * 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
  * Rees & Gunn (1974) Rees M. J., Gunn J. E., 1974, MNRAS, 167, 1
  * Reynolds et al. (2017) Reynolds S. P., Pavlov G. G., Kargaltsev O., Klingler N., Renaud M., Mereghetti S., 2017, Space Sci. Rev., 207, 175
  * Rozyczka & Franco (1996) Rozyczka M., Franco J., 1996, ApJ, 469, L127
  * 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
  * Slane (2017) Slane P., 2017, in Alsabti A. W., Murdin P., eds, , Handbook of Supernovae. p. 2159, doi:10.1007/978-3-319-21846-5_95
  * Steiner et al. (2005) Steiner A. W., Prakash M., Lattimer J. M., Ellis P. J., 2005, Phys. Rep., 411, 325
  * Stock & Barlow (2010) Stock D. J., Barlow M. J., 2010, MNRAS, 409, 1429
  * Sun et al. (1999) Sun M., Wang Z.-r., Chen Y., 1999, ApJ, 511, 274
  * Temim et al. (2017) Temim T., Slane P., Plucinsky P. P., Gelfand J., Castro D., Kolb C., 2017, ApJ, 851, 128
  * Temim et al. (2022) Temim T., Slane P., Raymond J. C., Patnaude D., Murray E., Ghavamian P., Renzo M., Jacovich T., 2022, ApJ, 932, 26
  * Toropina et al. (2019) Toropina O. D., Romanova M. M., Lovelace R. V. E., 2019, MNRAS, 484, 1475
  * Torres & Lin (2018) Torres D. F., Lin T., 2018, ApJ, 864, L2
  * Torres et al. (2014) Torres D., Cillis A., Martín J., de Oña Wilhelmi E., 2014, Journal of High Energy Astrophysics, 1-2, 31
  * Truelove & McKee (1999) Truelove J. K., McKee C. F., 1999, ApJS, 120, 299
  * Vaidya et al. (2018) Vaidya B., Mignone A., Bodo G., Rossi P., Massaglia S., 2018, ApJ, 865, 144
  * 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
  * Villagran et al. (2023) Villagran M. A., et al., 2023, MNRAS,
  * 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
  * Weber (1999) Weber F., 1999, Journal of Physics G Nuclear Physics, 25, R195
  * Weber & Davis (1967) Weber E. J., Davis Leverett J., 1967, ApJ, 148, 217
  * 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 (1996a) Wilkin F. P., 1996a, ApJ, 459, L31
  * Wilkin (1996b) Wilkin F. P., 1996b, ApJ, 459, L31
  * Williams et al. (2008) Williams B. J., et al., 2008, ApJ, 687, 1054
  * Wolfire et al. (2003) Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ, 587, 278
  * Yusef-Zadeh et al. (2003) Yusef-Zadeh F., Wardle M., Rho J., Sakano M., 2003, ApJ, 585, 319
  * de Vries & Romani (2020) de Vries M., Romani R. W., 2020, ApJ, 896, L7
  * de Vries et al. (2021) de Vries M., et al., 2021, ApJ, 908, 50
  * van Marle et al. (2015a) van Marle A. J., Meliani Z., Marcowith A., 2015a, A&A, 584, A49
  * van Marle et al. (2015b) van Marle A. J., Meliani Z., Marcowith A., 2015b, 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
  * van der Swaluw (2003) van der Swaluw E., 2003, A&A, 404, 939
  * van der Swaluw et al. (2001) van der Swaluw E., Achterberg A., Gallant Y. A., Tóth G., 2001, A&A, 380, 309
  * van der Swaluw et al. (2003) van der Swaluw E., Achterberg A., Gallant Y. A., Downes T. P., Keppens R., 2003, A&A, 397, 913
  * van der Swaluw et al. (2004) van der Swaluw E., Downes T. P., Keegan R., 2004, A&A, 420, 937