<figure><img src="image_1.png" /></figure>

HAL
open science

# On the plerionic rectangular supernova remnants of static progenitors

D.M.-A Meyer, Z Meliani, Velazquez P.F., M Pohl, D.F Torres

► To cite this version:

D.M.-A Meyer, Z Meliani, Velazquez P.F., M Pohl, D.F Torres. On the plerionic rectangular supernova remnants of static progenitors. Monthly Notices of the Royal Astronomical Society, 2023, 527 (3), pp.5514-5524. 10.1093/mnras/stad3495. hal-04324160

HAL Id: hal-04324160

https://hal.science/hal-04324160

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the plerionic rectangular supernova remnants of static progenitors

D.M.-A Meyer, Z Meliani, Velazquez P.F., M Pohl, D.F Torres

► To cite this version:

D.M.-A Meyer, Z Meliani, Velazquez P.F., M Pohl, D.F Torres. On the plerionic rectangular supernova remnants of static progenitors. Monthly Notices of the Royal Astronomical Society, 2023, 527 (3), pp.5514-5524. 10.1093/mnras/stad3495. hal-04324160

HAL Id: hal-04324160

https://hal.science/hal-04324160v1

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# On the plerionic rectangular supernova remnants of static progenitors

D. M.-A. Meyer<sup>1,★</sup> Z. Meliani<sup>2</sup> P. F. Velázquez<sup>3</sup> M. Pohl<sup>1,4</sup> and D. F. Torres<sup>5,6,7</sup>

Accepted 2023 November 7. Received 2023 October 16; in original form 2023 May 30

## ABSTRACT

Pulsar wind nebulae are a possible final stage of the circumstellar evolution of massive stars, where a fast-rotating, magnetized 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 magnetohydrodynamic simulations. Our models explore the impact of the interstellar medium’s (ISM) 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 ISM. 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.

**Key words:** MHD – stars: evolution – stars: massive – pulsars: general – ISM: supernova remnants.

## 1 INTRODUCTION

Pulsar wind nebulae (PWNe) are intriguing astrophysical phenomena powered by pulsars, rotating neutron stars with high magnetic field. 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; Steiner et al. 2005; Lasky 2015; Pnigouras & Kokkotas 2015, 2016; Pnigouras 2019).

Pulsars have very powerful magnetospheres with strong magnetic fields of 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 pulsar had been kicked off by the supernova explosion (Barkov et al. 2019b; Igoshev 2020; de Vries et al. 2021).

A particular class of supernova remnants that contain a pulsar exhibits a succession of structured shocks powered by the pulsar’s magnetic wind, which in turn produces multiwavelength polarized non-thermal emission. Notable examples of such plerions include the Crab Nebula (Kennel & Coroniti 1984; Gallant & Arons 1994; Aharonian et al. 2006) and the Geminga pulsar in the Vela supernova remnant (Bock, Turtle & Green 1998; Popov et al. 2019).

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 multidisciplinary approach. Secondly, the environment in which the pulsar wind is launched is often structured, as it depends on the properties of the supernova remnant, the circumstellar medium of the progenitor star, or the interstellar medium (ISM) 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 multifaceted problem, requiring sophisticated theoretical models and numerical simulations to understand the physics at play in full.

The Crab Nebula is one of the most well-known examples of a plerion. Extensive theoretical (Kennel & Coroniti 1984; Coroniti 1990; Begelman & Li 1992; Begelman 1998) and numerical studies have been conducted on the Crab Nebula in both 2.5D (e.g. Komissarov &

<footnote n="1"><sup>1</sup>Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, D-14476 Potsdam, Germany</footnote>

<footnote n="2"><sup>2</sup>Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, Université de Paris, CNRS, F-92190 Meudon, France</footnote>

<footnote n="3"><sup>3</sup>Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, CP 04510 Mexico City, Mexico</footnote>

<footnote n="4"><sup>4</sup>Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, D-15738 Zeuthen, Germany</footnote>

<footnote n="5"><sup>5</sup>Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans s/n, E-08193 Barcelona, Spain</footnote>

<footnote n="6"><sup>6</sup>Institut d'Estudis Espacials de Catalunya (IEEC), Gran Capità 2-4, E-08034 Barcelona, Spain</footnote>

<footnote n="7"><sup>7</sup>Institució Catalana de Recerca i Estudis Avançats (ICREA), E-08010 Barcelona, Spain</footnote>

<footnote>*E-mail: dmameyer.astro@gmail.com</footnote>
# On the plerionic rectangular supernova remnants of static progenitors

D. M.-A. Meyer<sup>1,★</sup> Z. Meliani<sup>2</sup> P. F. Velázquez<sup>3</sup> M. Pohl<sup>1,4</sup> and D. F. Torres<sup>5,6,7</sup>

Accepted 2023 November 7. Received 2023 October 16; in original form 2023 May 30

## ABSTRACT

Pulsar wind nebulae are a possible final stage of the circumstellar evolution of massive stars, where a fast-rotating, magnetized 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 magnetohydrodynamic simulations. Our models explore the impact of the interstellar medium’s (ISM) 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 ISM. 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.

**Key words:** MHD – stars: evolution – stars: massive – pulsars: general – ISM: supernova remnants.

## 1 INTRODUCTION

Pulsar wind nebulae (PWNe) are intriguing astrophysical phenomena powered by pulsars, rotating neutron stars with high magnetic field. 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; Steiner et al. 2005; Lasky 2015; Pnigouras & Kokkotas 2015, 2016; Pnigouras 2019).

Pulsars have very powerful magnetospheres with strong magnetic fields of 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 pulsar had been kicked off by the supernova explosion (Barkov et al. 2019b; Igoshev 2020; de Vries et al. 2021).

A particular class of supernova remnants that contain a pulsar exhibits a succession of structured shocks powered by the pulsar’s magnetic wind, which in turn produces multiwavelength polarized non-thermal emission. Notable examples of such plerions include the Crab Nebula (Kennel & Coroniti 1984; Gallant & Arons 1994; Aharonian et al. 2006) and the Geminga pulsar in the Vela supernova remnant (Bock, Turtle & Green 1998; Popov et al. 2019).

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 multidisciplinary approach. Secondly, the environment in which the pulsar wind is launched is often structured, as it depends on the properties of the supernova remnant, the circumstellar medium of the progenitor star, or the interstellar medium (ISM) 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 multifaceted problem, requiring sophisticated theoretical models and numerical simulations to understand the physics at play in full.

The Crab Nebula is one of the most well-known examples of a plerion. Extensive theoretical (Kennel & Coroniti 1984; Coroniti 1990; Begelman & Li 1992; Begelman 1998) and numerical studies have been conducted on the Crab Nebula in both 2.5D (e.g. Komissarov &

<footnote n="1"><sup>1</sup>Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, D-14476 Potsdam, Germany</footnote>

<footnote n="2"><sup>2</sup>Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, Université de Paris, CNRS, F-92190 Meudon, France</footnote>

<footnote n="3"><sup>3</sup>Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, CP 04510 Mexico City, Mexico</footnote>

<footnote n="4"><sup>4</sup>Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, D-15738 Zeuthen, Germany</footnote>

<footnote n="5"><sup>5</sup>Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans s/n, E-08193 Barcelona, Spain</footnote>

<footnote n="6"><sup>6</sup>Institut d'Estudis Espacials de Catalunya (IEEC), Gran Capità 2-4, E-08034 Barcelona, Spain</footnote>

<footnote n="7"><sup>7</sup>Institució Catalana de Recerca i Estudis Avançats (ICREA), E-08010 Barcelona, Spain</footnote>

<footnote>*E-mail: dmameyer.astro@gmail.com</footnote>
Lyubarsky 2003, 2004; Komissarov 2006; Komissarov & Lyutikov 2011) and 3D simulations (e.g. Porth, Komissarov & Keppens 2014; Olmi & Bucciantini 2019). Furthermore, the Crab Nebula provides important insights into the physics of pulsars and their interaction with supernova remnants.

The dynamics and morphology of PWNe undergo significant changes when they expand inside 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, Safi-Harb & Kothes 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 (Blondin, Chevalier & Frierson 2001; van der Swaluw 2003). These studies have also been extended to moving neutron stars inside the supernova ejecta, revealing the development of strong shocks between the pulsar wind and supernova remnant (van der Swaluw et al. 2003; van der Swaluw, Downes & Keegan 2004; Kolb et al. 2017; Temim et al. 2017, 2022).

In some cases, the strong interaction between the PWN and the reverse shock of the supernova remnant can result in a 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) and Bandiera et al. (2020, 2023)].

For later phases, when the pulsar leaves the supernova remnant and starts interacting with the ISM, a bow shock nebula forms around the runaway pulsar. This phenomenon has been extensively studied using 3D simulations (Bucciantini & Bandiera 2001; Bucciantini 2002a,b, 2018; Bucciantini et al. 2004; Barkov, Lyutikov & Khangulyan 2019a, 2020) and further detailed models have been tailored to the Crab Nebula (Mignone et al. 2013). 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 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, Meliani & Marcowith 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; Burns et al. 2023; Elbakyan et al. 2023). In its main-sequence phase, it expands into a wind bubble (Weaver et al. 1977; Gull & Sofia 1979; Wilkin 1996). 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. 2020). Once the pulsar emits a relativistic and powerful wind, it initially interacts with the surrounding supernova ejecta (Cox et al. 1999; Sun, Wang & Chen 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 the dynamics, morphologies, and emission properties of PWNe with static massive star progenitors are influenced by a magnetized ambient medium. The multidimensional magnetohydrodynamic (MHD) simulations conducted by van Marle, Meliani & Marcowith (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, Hensler & Yorke (2003, 2006) and Dwarkadas (2005, 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, leading to the formation of rectangular-shaped remnants like Puppis A. In this study, we further investigate the effects of the reflection of the supernova blast wave 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 PWNe.

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 PWNe, 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 set-up used in this study to generate models of PWNe 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_{ISM} = 0.79 \text{ cm}^{-3}$, while the magnetic field of the ISM is uniform and structured, with a strength of $B_{ISM} = 7 \mu\text{G}$. In these models, we assume that the ionized gas has a temperature of 8000 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, Schaye & Smith (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$ K, and various metals' emission lines at temperatures $T \ge 10^6$ K. The cooling curve is further enhanced with [O III] $\lambda$ 5007 collisionally excited forbidden lines, as described in Asplund et al. (2009) and Henney et al. (2009).

In this paper, we present a model that captures the evolution of the circumstellar medium surrounding a static massive star with an 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).
Lyubarsky 2003, 2004; Komissarov 2006; Komissarov & Lyutikov 2011) and 3D simulations (e.g. Porth, Komissarov & Keppens 2014; Olmi & Bucciantini 2019). Furthermore, the Crab Nebula provides important insights into the physics of pulsars and their interaction with supernova remnants.

The dynamics and morphology of PWNe undergo significant changes when they expand inside 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, Safi-Harb & Kothes 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 (Blondin, Chevalier & Frierson 2001; van der Swaluw 2003). These studies have also been extended to moving neutron stars inside the supernova ejecta, revealing the development of strong shocks between the pulsar wind and supernova remnant (van der Swaluw et al. 2003; van der Swaluw, Downes & Keegan 2004; Kolb et al. 2017; Temim et al. 2017, 2022).

In some cases, the strong interaction between the PWN and the reverse shock of the supernova remnant can result in a 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) and Bandiera et al. (2020, 2023)].

For later phases, when the pulsar leaves the supernova remnant and starts interacting with the ISM, a bow shock nebula forms around the runaway pulsar. This phenomenon has been extensively studied using 3D simulations (Bucciantini & Bandiera 2001; Bucciantini 2002a,b, 2018; Bucciantini et al. 2004; Barkov, Lyutikov & Khangulyan 2019a, 2020) and further detailed models have been tailored to the Crab Nebula (Mignone et al. 2013). 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 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, Meliani & Marcowith 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; Burns et al. 2023; Elbakyan et al. 2023). In its main-sequence phase, it expands into a wind bubble (Weaver et al. 1977; Gull & Sofia 1979; Wilkin 1996). 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. 2020). Once the pulsar emits a relativistic and powerful wind, it initially interacts with the surrounding supernova ejecta (Cox et al. 1999; Sun, Wang & Chen 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 the dynamics, morphologies, and emission properties of PWNe with static massive star progenitors are influenced by a magnetized ambient medium. The multidimensional magnetohydrodynamic (MHD) simulations conducted by van Marle, Meliani & Marcowith (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, Hensler & Yorke (2003, 2006) and Dwarkadas (2005, 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, leading to the formation of rectangular-shaped remnants like Puppis A. In this study, we further investigate the effects of the reflection of the supernova blast wave 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 PWNe.

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 PWNe, 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 set-up used in this study to generate models of PWNe 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_{ISM} = 0.79 \text{ cm}^{-3}$, while the magnetic field of the ISM is uniform and structured, with a strength of $B_{ISM} = 7 \mu\text{G}$. In these models, we assume that the ionized gas has a temperature of 8000 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, Schaye & Smith (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$ K, and various metals' emission lines at temperatures $T \ge 10^6$ K. The cooling curve is further enhanced with [O III] $\lambda$ 5007 collisionally excited forbidden lines, as described in Asplund et al. (2009) and Henney et al. (2009).

In this paper, we present a model that captures the evolution of the circumstellar medium surrounding a static massive star with an 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).
<table><caption>**Table 1.** List of models in this study. All simulations assume a rotating static massive star of mass $M_*$ at solar metallicity, in a medium of number density $n_{ISM}$ and organized magnetic field strength $B_{ISM}$. The initial rotation rate of the central massive star is $\Omega_*/\Omega_K = 0.1$.</caption><thead><tr><th>Model</th><th><i>M</i><sub>*</sub> (M<sub>☉</sub>)</th><th><i>n</i><sub>ISM</sub> (cm<sup>-1</sup>)</th><th><i>B</i><sub>ISM</sub> (μG)</th><th>Description</th></tr></thead><tbody><tr><td>Run-35-HD-0.79-PWN</td><td>35</td><td>0.79</td><td>0</td><td>35 M<sub>☉</sub> induced PWN in unmagnetized warm ISM</td></tr><tr><td>Run-35-MHD-0.79-PWN</td><td>35</td><td>0.79</td><td>7</td><td>35 M<sub>☉</sub> induced PWN in magnetized warm ISM</td></tr></tbody></table>

initial mass of 35 M☉ at the zero-age main sequence. The star is considered to be rotating with an angular velocity ratio of $\Omega_*/\Omega_K = 0.1$, where $\Omega_*$ represents the star's initial angular frequency and $\Omega_K$ is its equatorial Keplerian angular velocity. Consequently, the equatorial velocity of the star can be expressed as

$$v_{\text{rot}}(t) = \Omega_{*}(t)R_{*}(t). \qquad (1)$$

Here, $R_*(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 1D stellar evolution model provided by the GENEVA library, as described in Ekström et al. (2012).<fnref n="1" /> Specifically, we extract the mass-loss rate $\dot{M}(t)$ and the effective temperature $T_{\text{eff}}(t)$ of the star at each stage of evolution from this data base, and derive from it the wind density,

$$\rho_w(r, t) = \frac{\dot{M}(t)}{4\pi r^2 v_w(t)}. \qquad (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_w(t)$, is calculated based on the escape velocity $v_{\text{esc}}(t)$. The escape velocity depends on the star's effective temperature $T_{\text{eff}}$ and is determined using the conversion law,

$$v_w(t) = \sqrt{\beta(T)} v_{\text{esc}}(t) = \sqrt{\beta(T) \frac{2GM_*(t)}{R_*(t)}}, \qquad (3)$$

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

We adopt the time-dependent evolution of the surface magnetic field $B_*$ of the supernova progenitor as derived in Meyer et al. (2023) where the magnetic field strength at the surface of the star is scaled to that of the Sun, as described in Baalmann et al. (2020, 2021), Herbst et al. (2020), Scherer et al. (2020), and Meyer et al. (2021). Specifically, we assume a magnetic field strength at star surface of $B_* = 500$ G during the main-sequence phase (Castro et al. 2015, 2017; Fossati et al. 2015; Przybilla et al. 2016), to a Betelgeuse-like field of $B_* = 0.2$ G for the red supergiant phase (Vlemmings, Diamond & van Langevelde 2002; Vlemmings, van Langevelde & Diamond 2005; Kervella et al. 2018) and $B_* = 100$ 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_r(r, t) = B_*(t) \left( \frac{R_*(t)}{r} \right)^2, \qquad (4)$$

and a toroidal component,

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

respectively, with

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

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

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_{\text{max}}$. The ejecta has a total energy of $E_{\text{sn}} = 10^{51}$ erg and a mass of $M_{\text{sn}}$, which takes into account the star's mass loss throughout its entire evolution until the immediate pre-supernova time $t_{\text{psn}}$, as well as the mass $M_{\text{NS}}$ of the neutron star that forms at the centre. Specifically, we set

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

with $t_{\text{psn}}$ and $M_{\text{NS}} = 1.4 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_{\text{core}}$, where $r_{\text{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_{\text{core}}$ to $r_{\text{max}}$, where $r_{\text{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_{\text{core}} = \frac{1}{4\pi n} \frac{(10E_{\text{sn}}^{n-5})^{-3/2}}{(3M_{\text{sn}}^{n-3})^{-5/2}} \frac{1}{t_{\text{max}}^3} \qquad (8)$$

and

$$\rho_{\text{max}}(r) = \frac{1}{4\pi n} \frac{(10E_{\text{sn}}^{n-5})^{(n-3)/2}}{(3M_{\text{sn}}^{n-3})^{(n-5)/2}} \frac{1}{t_{\text{max}}^3} \left(\frac{r}{t_{\text{max}}}\right)^{-n}, \qquad (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_{\text{max}}$ at time $t_{\text{max}}$. The characteristics of the supernova ejecta profile are computed following the methodology outlined in Truelove & McKee (1999) and Whalen et al. (2008).

<footnote n="1">¹https://www.unige.ch/sciences/astro/evolution/en/database/syclist/</footnote>
<table><caption>**Table 1.** List of models in this study. All simulations assume a rotating static massive star of mass $M_*$ at solar metallicity, in a medium of number density $n_{ISM}$ and organized magnetic field strength $B_{ISM}$. The initial rotation rate of the central massive star is $\Omega_*/\Omega_K = 0.1$.</caption><thead><tr><th>Model</th><th><i>M</i><sub>*</sub> (M<sub>☉</sub>)</th><th><i>n</i><sub>ISM</sub> (cm<sup>-1</sup>)</th><th><i>B</i><sub>ISM</sub> (µG)</th><th>Description</th></tr></thead><tbody><tr><td>Run-35-HD-0.79-PWN</td><td>35</td><td>0.79</td><td>0</td><td>35 M<sub>☉</sub> induced PWN in unmagnetized warm ISM</td></tr><tr><td>Run-35-MHD-0.79-PWN</td><td>35</td><td>0.79</td><td>7</td><td>35 M<sub>☉</sub> induced PWN in magnetized warm ISM</td></tr></tbody></table>

initial mass of 35 M☉ at the zero-age main sequence. The star is considered to be rotating with an angular velocity ratio of $\Omega_*/\Omega_K = 0.1$, where $\Omega_*$ represents the star's initial angular frequency and $\Omega_K$ is its equatorial Keplerian angular velocity. Consequently, the equatorial velocity of the star can be expressed as

$$v_{\text{rot}}(t) = \Omega_{*}(t)R_{*}(t). \qquad (1)$$

Here, $R_*(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 1D stellar evolution model provided by the GENEVA library, as described in Ekström et al. (2012).<fnref n="1" /> Specifically, we extract the mass-loss rate $\dot{M}(t)$ and the effective temperature $T_{\text{eff}}(t)$ of the star at each stage of evolution from this data base, and derive from it the wind density,

$$\rho_w(r, t) = \frac{\dot{M}(t)}{4\pi r^2 v_w(t)}. \qquad (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_w(t)$, is calculated based on the escape velocity $v_{\text{esc}}(t)$. The escape velocity depends on the star's effective temperature $T_{\text{eff}}$ and is determined using the conversion law,

$$v_w(t) = \sqrt{\beta(T)} v_{\text{esc}}(t) = \sqrt{\beta(T) \frac{2GM_*(t)}{R_*(t)}}, \qquad (3)$$

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

We adopt the time-dependent evolution of the surface magnetic field $B_*$ of the supernova progenitor as derived in Meyer et al. (2023) where the magnetic field strength at the surface of the star is scaled to that of the Sun, as described in Baalmann et al. (2020, 2021), Herbst et al. (2020), Scherer et al. (2020), and Meyer et al. (2021). Specifically, we assume a magnetic field strength at star surface of $B_* = 500$ G during the main-sequence phase (Castro et al. 2015, 2017; Fossati et al. 2015; Przybilla et al. 2016), to a Betelgeuse-like field of $B_* = 0.2$ G for the red supergiant phase (Vlemmings, Diamond & van Langevelde 2002; Vlemmings, van Langevelde & Diamond 2005; Kervella et al. 2018) and $B_* = 100$ 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_r(r, t) = B_*(t) \left( \frac{R_*(t)}{r} \right)^2, \qquad (4)$$

and a toroidal component,

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

respectively, with

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

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

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_{\text{max}}$. The ejecta has a total energy of $E_{\text{sn}} = 10^{51}$ erg and a mass of $M_{\text{sn}}$, which takes into account the star's mass loss throughout its entire evolution until the immediate pre-supernova time $t_{\text{psn}}$, as well as the mass $M_{\text{NS}}$ of the neutron star that forms at the centre. Specifically, we set

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

with $t_{\text{psn}}$ and $M_{\text{NS}} = 1.4 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_{\text{core}}$, where $r_{\text{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_{\text{core}}$ to $r_{\text{max}}$, where $r_{\text{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_{\text{core}} = \frac{1}{4\pi n} \frac{(10E_{\text{sn}}^{n-5})^{-3/2}}{(3M_{\text{sn}}^{n-3})^{-5/2}} \frac{1}{t_{\text{max}}^3} \qquad (8)$$

and

$$\rho_{\text{max}}(r) = \frac{1}{4\pi n} \frac{(10E_{\text{sn}}^{n-5})^{(n-3)/2}}{(3M_{\text{sn}}^{n-3})^{(n-5)/2}} \frac{1}{t_{\text{max}}^3} \left(\frac{r}{t_{\text{max}}}\right)^{-n}, \qquad (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_{\text{max}}$ at time $t_{\text{max}}$. The characteristics of the supernova ejecta profile are computed following the methodology outlined in Truelove & McKee (1999) and Whalen et al. (2008).

<footnote n="1">¹https://www.unige.ch/sciences/astro/evolution/en/database/syclist/</footnote>
The velocity at the core radius, denoted as $v_{\text{core}}(r_{\text{core}})$, is determined as

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

where $E_{\text{sn}}$ represents the total energy of the supernova ejecta and $M_{\text{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_{\text{max}}$, is set to

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

This choice of $v_{\text{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}_{\text{psw}} = \dot{E}_{\text{psw},0} \left(1 + \frac{t}{\tau_0}\right)^{-\frac{n+1}{n-1}}. \quad (12)$$

Here, $\tau_0$ is the initial spin-down of the pulsar, given by

$$\tau_0 = \frac{P_o}{(n-1)\dot{P}_o}. \quad (13)$$

The initial total power of the pulsar wind is set to $\dot{E}_{\text{psw},0} = 10^{38} \text{ erg s}^{-1}$. The initial spin period of the pulsar is set to $P_o = 0.3 \text{ s}$, and its time derivative is set to $\dot{P}_o = 10^{-17} \text{ s s}^{-1}$. The braking index is assumed to be $n=3$, which corresponds to magnetic dipole spin-down, as outlined in Slane (2017).

Moreover, the pulsar's wind is assumed to have a constant speed of $v_{\text{psw}} = 10^{-2}c$ ($c$ is the speed of light) and magnetization of $\sigma = 10^{-3}$, as outlined in Slane (2017). 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).

$$\dot{E}_{\text{total}} = \frac{\dot{E}_{\text{psw}}}{r^2} (\sin^2 \theta + 1/\sigma), \quad (14)$$

$$B = \sqrt{\frac{4\pi}{c}} \frac{1}{r} \sin\theta (1 - 2\theta/\pi), \quad (15)$$

$$\dot{E}_{\text{kinetic}} = \dot{E}_{\text{total}} - \frac{B^2}{4\pi} c. \quad (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) and Meyer, Petrov & Pohl (2020), and that we extended after to PWN in Meyer & Meliani (2022). The MHD simulations are conducted in a 2.5D, axisymmetric cylindrical coordinate system. The simulation box extends over the range $[0; R_{\text{max}}] \times [z_{\text{min}}; z_{\text{max}}]$ and is discretized using a uniform grid of $N_R \times N_z$ cells. Consequently, the spatial resolution is consistent along both directions, with each grid cell having a size of $\Delta = R_{\text{max}}/N_R$. We employ two different spatial resolutions throughout the evolutionary process. During the

phases involving the progenitor star wind, the circumstellar medium is resolved using a grid resolution of $N_R = 2000$ and $N_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_R = 3000$ and $N_z = 6000$ cells. The supernova ejecta is confined within a central sphere of radius $r_{\text{max}}$, as described in Section 2.1. Simultaneously, the pulsar wind is imposed within a sphere of radius $r_{\text{nswind}} = 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 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 magnetized 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_{\text{ISM}} = 7 \mu\text{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.

The numerical simulations are conducted using the PLUTO code (Mignone et al. 2007, 2012; Vaidya et al. 2018)<fnref n="2" /> and we solve the following set of equations,

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \quad (17)$$

$$\frac{\partial \mathbf{m}}{\partial t} + \nabla \cdot (\mathbf{m} \otimes \mathbf{v} - \mathbf{B} \otimes \mathbf{B} + \hat{\mathbf{I}} p_t) \mathbf{0}, \quad (18)$$

$$\frac{\partial E}{\partial t} + \nabla \cdot \left( (E + p_t)\mathbf{v} - \mathbf{B}(\mathbf{v} \cdot \mathbf{B}) \right) = \Phi(T, \rho), \quad (19)$$

and

$$\frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{v} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{v}) = \mathbf{0}, \quad (20)$$

with the gas density $\rho$, velocity $v$, momentum $m = \rho v$, and magnetic field $\mathbf{B}$, as well as the total pressure $p_t$ and the energy of the gas

$$E = \frac{p}{(\gamma - 1)} + \frac{\mathbf{m} \cdot \mathbf{m}}{2\rho} + \frac{\mathbf{B} \cdot \mathbf{B}}{2}. \quad (21)$$

The sound speed of the medium reads

$$c_s = \sqrt{\frac{\gamma p}{\rho}}, \quad (22)$$

where the adiabatic index is $\gamma = 5/3$. Lastly, radiative cooling by optically thin processes and photoheating 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 and utilize

<footnote n="2">²http://plutocode.ph.unito.it/</footnote>
The velocity at the core radius, denoted as $v_{\text{core}}(r_{\text{core}})$, is determined as

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

where $E_{\text{sn}}$ represents the total energy of the supernova ejecta and $M_{\text{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_{\text{max}}$, is set to

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

This choice of $v_{\text{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}_{\text{psw}} = \dot{E}_{\text{psw},0} \left(1 + \frac{t}{\tau_0}\right)^{-\frac{n+1}{n-1}}. \qquad (12)$$

Here, $\tau_0$ is the initial spin-down of the pulsar, given by

$$\tau_0 = \frac{P_o}{(n-1)\dot{P}_o}. \qquad (13)$$

The initial total power of the pulsar wind is set to $\dot{E}_{\text{psw},0} = 10^{38} \text{ erg s}^{-1}$. The initial spin period of the pulsar is set to $P_o = 0.3 \text{ s}$, and its time derivative is set to $\dot{P}_o = 10^{-17} \text{ s s}^{-1}$. The braking index is assumed to be $n=3$, which corresponds to magnetic dipole spin-down, as outlined in Slane (2017).

Moreover, the pulsar's wind is assumed to have a constant speed of $v_{\text{psw}} = 10^{-2}c$ ($c$ is the speed of light) and magnetization of $\sigma = 10^{-3}$, as outlined in Slane (2017). 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).

$$\dot{E}_{\text{total}} = \frac{\dot{E}_{\text{psw}}}{r^2} (\sin^2 \theta + 1/\sigma), \qquad (14)$$

$$B = \sqrt{\frac{4\pi}{c} \frac{1}{r} \sin\theta (1 - 2\theta/\pi)}, \qquad (15)$$

$$\dot{E}_{\text{kinetic}} = \dot{E}_{\text{total}} - \frac{B^2}{4\pi} c. \qquad (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) and Meyer, Petrov & Pohl (2020), and that we extended after to PWN in Meyer & Meliani (2022). The MHD simulations are conducted in a 2.5D, axisymmetric cylindrical coordinate system. The simulation box extends over the range $[0; R_{\text{max}}] \times [z_{\text{min}}; z_{\text{max}}]$ and is discretized using a uniform grid of $N_R \times N_z$ cells. Consequently, the spatial resolution is consistent along both directions, with each grid cell having a size of $\Delta = R_{\text{max}}/N_R$. We employ two different spatial resolutions throughout the evolutionary process. During the

phases involving the progenitor star wind, the circumstellar medium is resolved using a grid resolution of $N_R = 2000$ and $N_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_R = 3000$ and $N_z = 6000$ cells. The supernova ejecta is confined within a central sphere of radius $r_{\text{max}}$, as described in Section 2.1. Simultaneously, the pulsar wind is imposed within a sphere of radius $r_{\text{nswind}} = 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 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 magnetized 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_{\text{ISM}} = 7 \mu\text{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.

The numerical simulations are conducted using the PLUTO code (Mignone et al. 2007, 2012; Vaidya et al. 2018)$^2$ and we solve the following set of equations,

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (17)$$

$$\frac{\partial \mathbf{m}}{\partial t} + \nabla \cdot (\mathbf{m} \otimes \mathbf{v} - \mathbf{B} \otimes \mathbf{B} + \hat{\mathbf{I}} p_t) \mathbf{0}, \qquad (18)$$

$$\frac{\partial E}{\partial t} + \nabla \cdot \left( (E + p_t)\mathbf{v} - \mathbf{B}(\mathbf{v} \cdot \mathbf{B}) \right) = \Phi(T, \rho), \qquad (19)$$

and

$$\frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{v} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{v}) = \mathbf{0}, \qquad (20)$$

with the gas density $\rho$, velocity $v$, momentum $\mathbf{m} = \rho\mathbf{v}$, and magnetic field $\mathbf{B}$, as well as the total pressure $p_t$ and the energy of the gas

$$E = \frac{p}{(\gamma - 1)} + \frac{\mathbf{m} \cdot \mathbf{m}}{2\rho} + \frac{\mathbf{B} \cdot \mathbf{B}}{2}. \qquad (21)$$

The sound speed of the medium reads

$$c_s = \sqrt{\frac{\gamma p}{\rho}}, \qquad (22)$$

where the adiabatic index is $\gamma = 5/3$. Lastly, radiative cooling by optically thin processes and photoheating 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 and utilize

<footnote>$^2$http://plutocode.ph.unito.it/</footnote>
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 number.

The numerical simulations are performed at the North-German Supercomputing Alliance (HLRN<fnref n="3" />) using the LISE cluster in Berlin, which is equipped with Cray XC40/30 processors.

# 3 RESULTS

In this section, we will analyse 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 unmagnetized ISM

In Fig. 1, the density contour is shown for the unmagnetized case Run-35-HD-0.79-PWN (left-hand panels) and the magnetized case Run-35-MHD-0.79-PWN (right-hand panels) at different evolution times, from top to the bottom. The density contour is represented in logarithmic scale in units of cm⁻<fnref n="3" />. 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 Fig. 1(a) (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 pc. Throughout the evolution of the star, 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 ≥50 pc, is the result of the interaction between the stellar wind and the ISM, and it occurs mainly during the main-sequence phase (Freyer, Hensler & Yorke 2003, 2006; Dwarkadas 2005, 2007). In the central region, within a radius of less than 20 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 pc, is dense and exhibits unstable behaviour. Subsequently, a second, less dense shell is formed as a result of the interaction between the red-giant wind and the main-sequence wind. Additionally, the main-sequence wind interacts with the surrounding ambient medium, leading to the formation of 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 rotation of the star. The variations between the bubbles depicted in Fig. 1(a) 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. 1(c) (middle-left), we can observe the supernova remnant at 25 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 partial reverberation. This interaction and reverberation contribute to the observed structure and morphology of the supernova remnant.

After the supernova explosion, a powerful pulsar wind is launched with a mechanical luminosity of $E = 10^{38}$ erg s⁻¹. This pulsar wind interacts with the dense supernova ejecta (van der Swaluw, Downes & Keegan 2004), resulting in the formation of a complex structure as described in Meyer & Meliani (2022) for a runaway progenitor star with a mass of $20 M_⊙$. 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 reverberation of the pulsar wind. At the interface between the unperturbed pulsar wind and the shocked pulsar wind, a pulsar wind termination shock is formed. The outermost region of the PWN is characterized by a contact discontinuity, which marks the interface between the ejecta and the pulsar wind (depicted by the red contour in the figure). Beyond the contact discontinuity, a transmitted shock propagates through the still unshocked ejecta, further travelling 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 diamond-like shape, as found by Komissarov & Lyubarsky (2004). This shape can be observed at a later time, specifically 45 kyr after the explosion, as shown in Fig. 1(e). 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. 1(e).

## 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. It is during this phase that 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 presence of the wind bubble influences the dynamics of the pulsar wind, 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. 1(b) (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. 1(a)]. 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, leading to an increase in magnetic pressure and tension along the polar axis. This phenomenon has been extensively studied and described in

<footnote n="3">³https://www.hlrn.de/</footnote>
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 number.

The numerical simulations are performed at the North-German Supercomputing Alliance (HLRN<fnref n="3" />) using the LISE cluster in Berlin, which is equipped with Cray XC40/30 processors.

# 3 RESULTS

In this section, we will analyse 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 unmagnetized ISM

In Fig. 1, the density contour is shown for the unmagnetized case Run-35-HD-0.79-PWN (left-hand panels) and the magnetized case Run-35-MHD-0.79-PWN (right-hand panels) at different evolution times, from top to the bottom. The density contour is represented in logarithmic scale in units of cm⁻<fnref n="3" />. 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 Fig. 1(a) (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 pc. Throughout the evolution of the star, 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 ≥50 pc, is the result of the interaction between the stellar wind and the ISM, and it occurs mainly during the main-sequence phase (Freyer, Hensler & Yorke 2003, 2006; Dwarkadas 2005, 2007). In the central region, within a radius of less than 20 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 pc, is dense and exhibits unstable behaviour. Subsequently, a second, less dense shell is formed as a result of the interaction between the red-giant wind and the main-sequence wind. Additionally, the main-sequence wind interacts with the surrounding ambient medium, leading to the formation of 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 rotation of the star. The variations between the bubbles depicted in Fig. 1(a) 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. 1(c) (middle-left), we can observe the supernova remnant at 25 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 partial reverberation. This interaction and reverberation contribute to the observed structure and morphology of the supernova remnant.

After the supernova explosion, a powerful pulsar wind is launched with a mechanical luminosity of $E = 10^{38}$ erg s⁻¹. This pulsar wind interacts with the dense supernova ejecta (van der Swaluw, Downes & Keegan 2004), resulting in the formation of a complex structure as described in Meyer & Meliani (2022) for a runaway progenitor star with a mass of $20 M_⊙$. 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 reverberation of the pulsar wind. At the interface between the unperturbed pulsar wind and the shocked pulsar wind, a pulsar wind termination shock is formed. The outermost region of the PWN is characterized by a contact discontinuity, which marks the interface between the ejecta and the pulsar wind (depicted by the red contour in the figure). Beyond the contact discontinuity, a transmitted shock propagates through the still unshocked ejecta, further travelling 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 diamond-like shape, as found by Komissarov & Lyubarsky (2004). This shape can be observed at a later time, specifically 45 kyr after the explosion, as shown in Fig. 1(e). 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. 1(e).

## 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. It is during this phase that 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 presence of the wind bubble influences the dynamics of the pulsar wind, 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. 1(b) (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. 1(a)]. 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, leading to an increase in magnetic pressure and tension along the polar axis. This phenomenon has been extensively studied and described in

<footnote n="3">³https://www.hlrn.de/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Number density fields in our MHD simulation of the PWN forming in the supernova remnant of a static 35 M☉ star rotating with $\Omega_*/\Omega_K = 0.1$ in an unmagnetized (left) and magnetized (right) ISM. The red contours highlight the region with a 50 per cent contribution of pulsar wind material. The streamlines in the right-hand side of panels (b), (d), and (f) mark the ISM magnetic field lines.</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Number density fields in our MHD simulation of the PWN forming in the supernova remnant of a static 35 M☉ star rotating with $\Omega_*/\Omega_K = 0.1$ in an unmagnetized (left) and magnetized (right) ISM. The red contours highlight the region with a 50 per cent contribution of pulsar wind material. The streamlines in the right-hand side of panels (b), (d), and (f) mark the ISM magnetic field lines.</figcaption></figure>
detail in van Marle, Meliani & Marcowith (2015b). During the last evolution phase, when the Wolf–Rayet wind material reaches the main-sequence terminal 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-hand 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. 1(d) 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 blast wave, also affects the confinement and shape of the pulsar wind. The reflected and the supernova blast wave, under the influence of the ISM magnetic field, 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 to 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. 1(d) and (f).

Indeed, the presence of a magnetized ISM influences the propagation of the PWN, resulting in distinct behaviour 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. 1(c) and (d). As time progresses, at a later evolution time of 45 kyr as depicted in Fig. 1(f), the pulsar wind continues to be channelled 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 ISM’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 PWNe and their interaction with the surrounding environment.

# 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 PWNe and compare our findings to existing observational data. By doing so, we aim to provide a comprehensive analysis of the strengths and weaknesses of our model and to assess its compatibility with the observed properties of PWNe.

## 4.1 Model limitations

Let us first consider two aspects central to the model. First, the simulations conducted in this study are 2D, which means they assume axisymmetry and do not account for variations in the supernova progenitor or the pulsar’s spin. While this approach allows for computational efficiency and provides valuable insights, it should be acknowledged that a fully 3D treatment would be necessary to capture the complete realism of the ISM. In particular, the inclusion of magnetic turbulence and the multiphased structure of the ISM could influence the development of the stellar wind bubble and the expansion of the supernova blast wave. Thus, future investigations employing a 3D approach would offer a more comprehensive understanding of these phenomena.

Another aspect to consider is the absence of pulsar motion in the simulations. Incorporating pulsar motion would introduce additional complexities and provide 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 PWN.

These are avenues for future research. The chosen 2D set-up and static pulsar position provide valuable insights into the general behaviour and trends, but future investigations can explore the impact of 3D effects and pulsar motion to obtain a more comprehensive characterization of the system’s dynamics and morphology.

## 4.2 Non-thermal emission

In order 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}, \qquad (23)$$

where 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_{\text{sync}}(\nu) \propto n B_{\perp}^{(s+1)/2} \nu^{-(s-1)/2}, \qquad (24)$$

$\nu$ being 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_{\text{SNR}} j_{\text{sync}}(\theta_{\text{obs}}) dl, \qquad (25)$$
detail in van Marle, Meliani & Marcowith (2015b). During the last evolution phase, when the Wolf–Rayet wind material reaches the main-sequence terminal 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-hand 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. 1(d) 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 blast wave, also affects the confinement and shape of the pulsar wind. The reflected and the supernova blast wave, under the influence of the ISM magnetic field, 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 to 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. 1(d) and (f).

Indeed, the presence of a magnetized ISM influences the propagation of the PWN, resulting in distinct behaviour 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. 1(c) and (d). As time progresses, at a later evolution time of 45 kyr as depicted in Fig. 1(f), the pulsar wind continues to be channelled 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 ISM’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 PWNe and their interaction with the surrounding environment.

# 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 PWNe and compare our findings to existing observational data. By doing so, we aim to provide a comprehensive analysis of the strengths and weaknesses of our model and to assess its compatibility with the observed properties of PWNe.

## 4.1 Model limitations

Let us first consider two aspects central to the model. First, the simulations conducted in this study are 2D, which means they assume axisymmetry and do not account for variations in the supernova progenitor or the pulsar’s spin. While this approach allows for computational efficiency and provides valuable insights, it should be acknowledged that a fully 3D treatment would be necessary to capture the complete realism of the ISM. In particular, the inclusion of magnetic turbulence and the multiphased structure of the ISM could influence the development of the stellar wind bubble and the expansion of the supernova blast wave. Thus, future investigations employing a 3D approach would offer a more comprehensive understanding of these phenomena.

Another aspect to consider is the absence of pulsar motion in the simulations. Incorporating pulsar motion would introduce additional complexities and provide 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 PWN.

These are avenues for future research. The chosen 2D set-up and static pulsar position provide valuable insights into the general behaviour and trends, but future investigations can explore the impact of 3D effects and pulsar motion to obtain a more comprehensive characterization of the system’s dynamics and morphology.

## 4.2 Non-thermal emission

In order 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}, \qquad (23)$$

where 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_{\text{sync}}(\nu) \propto n B_{\perp}^{(s+1)/2} \nu^{-(s-1)/2}, \qquad (24)$$

ν being the observed frequency and $B_⊥$ 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_{\text{SNR}} j_{\text{sync}}(\theta_{\text{obs}}) dl, \qquad (25)$$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Normalized radio synchrotron emission map of the plerionic supernova remnants with an inclination angle of $\theta_{\text{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 magnetized ISM (Run-35-MHD-0.79-PWN). The top figures show the remnants at time 25 kyr and the bottom figures display them at time 45 kyr.</figcaption></figure>

where $\theta_{\text{obs}}$ denotes the inclination angle of the remnant with respect to the sky plane. These calculations were conducted using the radiative transfer code RADMC-3D,<fnref n="4" /> and the methodology described in detail by Meyer et al. (2022a).

Fig. 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 a time of 25 kyr, while the bottom panels depict a time of 45 kyr. The intensity is plotted assuming an observer angle ($\theta_{\text{obs}}$) of 45°, representing the angle between the plane of the sky and the plane of symmetry of the supernova remnant.

Fig. 2(a) displays the PWN at an age of 25 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

<footnote n="4">⁴https://www.ita.uni-heidelberg.de/dullemond/software/radmc-3d/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Normalized radio synchrotron emission map of the plerionic supernova remnants with an inclination angle of $\theta_{\text{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 magnetized ISM (Run-35-MHD-0.79-PWN). The top figures show the remnants at time 25 kyr and the bottom figures display them at time 45 kyr.</figcaption></figure>

where $\theta_{\text{obs}}$ denotes the inclination angle of the remnant with respect to the sky plane. These calculations were conducted using the radiative transfer code RADMC-3D,<fnref n="4" /> and the methodology described in detail by Meyer et al. (2022a).

Fig. 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 a time of 25 kyr, while the bottom panels depict a time of 45 kyr. The intensity is plotted assuming an observer angle ($\theta_{\text{obs}}$) of 45°, representing the angle between the plane of the sky and the plane of symmetry of the supernova remnant.

Fig. 2(a) displays the PWN at an age of 25 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

<footnote n="4">⁴https://www.ita.uni-heidelberg.de/dullemond/software/radmc-3d/</footnote>
the polar regions 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 kyr, the radio synchrotron map of the PWN in an unmagnetized ISM is shown in Fig. 2(c). The PWN exhibits a diamond-like shape, with brighter regions observed at the polar zones. These bright regions are a result of 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. In this region, the gas is more diluted, 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. At an age of 25 kyr, the corresponding image is shown in Fig. 2(b). It reveals the presence of two bright arcs parallel to the direction of the ISM magnetic field. These arcs, observed in our axisymmetric set-up 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 leading to the formation of a PWN enclosed within a rectangular supernova remnant.

Fig. 2(d) 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 PWN 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

### 4.3.1 Generic discussion

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. Our objective is to investigate the formation of elongated PWNe, similar to those observed in Igoshev (2020). It is important to note that these elongated PWNe, 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 (De Luca et al. 2013).

The class of torus/jet-like PWNe, 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. MHD 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 PWN, which exhibits two curved antennae representing its jets/counter-jets 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 PWNe, which cannot be classified as either torus/jet-like objects or bow shock/tail systems, may be the result of 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).

### 4.3.2 The case of the Boomerang nebula

The plerionic supernova remnant G106.3+2.7, which hosts the Boomerang PWN associated with the pulsar J2229+6114, is worth further discussion in the context of our numerical models.

G106.3+2.7 is an old supernova remnant with a size of approximately 200 pc, assuming a distance of 12 kpc. Its age has been constrained to 1.26 Myr, as reported by Pineault & Joncas (2000), which presents continuum millimetre observations conducted by the Dominion Radio Astrophysical Observatory Synthesis Telescope. The remnant is in its final expansion stage, with a forward shock speed of the order of 15 km s⁻¹, interacting with a complex environment consisting of hot gas, warm material, cold molecular material, and a nearby H II region (Kothes, Uyaniker & Pineault 2001).

The Boomerang PWN associated with J2229+6114 is characterized by a bright synchrotron arc that borders and distorts the outermost shock of the expanding supernova remnant. The pulsar wind passes through the shock, leaving the remnant and entering the unperturbed ISM (Kothes, Uyaniker & Pineault 2001). This observation suggests that the off-centring of the pulsar within the remnant is not solely due to the typical high velocity of a neutron star kicked during a supernova explosion (Kothes, Reich & Uyaniker 2006).

## 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.5D simulations, we demonstrate that the reflection of the supernova blast wave against the elongated contact discontinuity between the stellar wind and magnetized 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 the formation of 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 diamond-like shape typically observed in previous studies, as described by Komissarov & Lyubarsky (2004).
the polar regions 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 kyr, the radio synchrotron map of the PWN in an unmagnetized ISM is shown in Fig. 2(c). The PWN exhibits a diamond-like shape, with brighter regions observed at the polar zones. These bright regions are a result of 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. In this region, the gas is more diluted, 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. At an age of 25 kyr, the corresponding image is shown in Fig. 2(b). It reveals the presence of two bright arcs parallel to the direction of the ISM magnetic field. These arcs, observed in our axisymmetric set-up 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 leading to the formation of a PWN enclosed within a rectangular supernova remnant.

Fig. 2(d) 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 PWN 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

### 4.3.1 Generic discussion

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. Our objective is to investigate the formation of elongated PWNe, similar to those observed in Igoshev (2020). It is important to note that these elongated PWNe, 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 (De Luca et al. 2013).

The class of torus/jet-like PWNe, 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. MHD 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 PWN, which exhibits two curved antennae representing its jets/counter-jets 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 PWNe, which cannot be classified as either torus/jet-like objects or bow shock/tail systems, may be the result of 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).

### 4.3.2 The case of the Boomerang nebula

The plerionic supernova remnant G106.3+2.7, which hosts the Boomerang PWN associated with the pulsar J2229+6114, is worth further discussion in the context of our numerical models.

G106.3+2.7 is an old supernova remnant with a size of approximately 200 pc, assuming a distance of 12 kpc. Its age has been constrained to 1.26 Myr, as reported by Pineault & Joncas (2000), which presents continuum millimetre observations conducted by the Dominion Radio Astrophysical Observatory Synthesis Telescope. The remnant is in its final expansion stage, with a forward shock speed of the order of 15 km s⁻¹, interacting with a complex environment consisting of hot gas, warm material, cold molecular material, and a nearby H II region (Kothes, Uyaniker & Pineault 2001).

The Boomerang PWN associated with J2229+6114 is characterized by a bright synchrotron arc that borders and distorts the outermost shock of the expanding supernova remnant. The pulsar wind passes through the shock, leaving the remnant and entering the unperturbed ISM (Kothes, Uyaniker & Pineault 2001). This observation suggests that the off-centring of the pulsar within the remnant is not solely due to the typical high velocity of a neutron star kicked during a supernova explosion (Kothes, Reich & Uyaniker 2006).

## 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.5D simulations, we demonstrate that the reflection of the supernova blast wave against the elongated contact discontinuity between the stellar wind and magnetized 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 the formation of 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 diamond-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 PWNe 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. The interaction between the supernova ejecta and the pulsar wind with this complex environment gives rise to the observed irregular morphologies.

In future research, the numerical model presented in this study will be refined and customized to specific supernova remnants hosting PWNe, such as the case of the old plerion G106.3+2.7 and its central pulsar J2229+6114, which drives the Boomerang synchrotron shell. By tailoring the model to these specific systems, researchers aim to achieve a more detailed and accurate understanding of the morphologies and dynamics of these objects.

## ACKNOWLEDGEMENTS

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 grant PID2021-124581OB-I00 funded by MCIN/AEI/10.13039/501100011033, 2021SGR00426 of the Generalitat de Catalunya, by the programme Unidad de Excelencia María de Maestu CEX2020-001058-M, and by MCIN with funding from European Union NextGeneration EU (PRTR-C17.11).

## 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.itai.uni-heidelberg.de/~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

* Aharonian F. et al., 2006, A&A, 457, 899

* Arias M. et al., 2019, A&A, 622, A6

* Aschenbach B., Leahy D. A., 1999, A&A, 341, 602

* Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

* Baalmann L. R., Scherer K., Fichtner H., Kleimann J., Bomans D. J., Weis K., 2020, A&A, 634, A67

* Baalmann L. R., Scherer K., Kleimann J., Fichtner H., Bomans D. J., Weis K., 2021, A&A, 650, A36

* Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2020, MNRAS, 499, 2051

* Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2021, MNRAS, 508, 3194

* Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2023, MNRAS, 520, 2451

* Barkov M. V., Lyutikov M., Khangulyan D., 2019a, MNRAS, 484, 4760

* Barkov M. V., Lyutikov M., Klingler N., Bordas P., 2019b, MNRAS, 485, 2041

* Barkov M. V., Lyutikov M., Khangulyan D., 2020, MNRAS, 497, 2605

* Begelman M. C., 1998, ApJ, 493, 291

* Begelman M. C., Li Z.-Y., 1992, ApJ, 397, 187

* Blondin J. M., Chevalier R. A., 2017, ApJ, 845, 139

* Blondin J. M., Chevalier R. A., Frierson D. M., 2001, ApJ, 563, 806

* Bock D. C. J., Turtle A. J., Green A. J., 1998, AJ, 116, 1886

* Bucciantini N., 2002a, A&A, 387, 1066

* Bucciantini N., 2002b, A&A, 393, 629

* Bucciantini N., 2018, MNRAS, 478, 2074

* Bucciantini N., Bandiera R., 2001, A&A, 375, 1032

* Bucciantini N., Bandiera R., Blondin J. M., Amato E., Del Zanna L., 2004, A&A, 422, 609

* Burns R. A. et al., 2023, Nat. Astron., 7, 557

* Castro N. et al., 2015, A&A, 581, A81

* Castro N. et al., 2017, A&A, 597, L6

* Chevalier R. A., 1982, ApJ, 258, 790

* Chevalier R. A., Luo D., 1994, ApJ, 421, 225

* Chiotellis A., Boumis P., Derlopa S., Steffen W., 2019, Proceedings of the conference "Supernova Remnants II: An Odyssey in Space after Stellar Death", preprint (arXiv:1909.08947)

* Comerón F., Kaper L., 1998, A&A, 338, 273

* Coroniti F. V., 1990, ApJ, 349, 538

* Cox D. P., Shelton R. L., Maciejewski W., Smith R. K., Plewa T., Pawl A., Różyczka M., 1999, ApJ, 524, 179

* Cox N. L. J., Kerschbaum F., van Marle A. J., Decin L., Ladjal D., Mayer A., 2012, A&A, 543, C1

* Crawford F., Gaensler B. M., Kaspi V. M., Manchester R. N., Camilo F., Lyne A. G., Pivovaroff M. J., 2001, ApJ, 554, 152

* Das S., Brose R., Meyer D. M. A., Pohl M., Sushch I., Plotko P., 2022, A&A, 661, A128

* Decin L., 2012, Adv. Space Res., 50, 843

* Decin L. N. L. J. et al., 2012, A&A, 548, A113

* De Luca A. et al., 2013, ApJ, 765, L19

* Derlopa S., Boumis P., Chiotellis A., Steffen W., Akras S., 2020, MNRAS, 499, 5410

* de Vries M. et al., 2021, ApJ, 908, 50

* Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium. Princeton Univ. Press, Princeton, NJ

* Dwarkadas V. V., 2005, ApJ, 630, 892

* Dwarkadas V. V., 2007, ApJ, 667, 226

* Ekström S. et al., 2012, A&A, 537, A146

* Elbakyan V. G., Nayakshin S., Meyer D. M. A., Vorobyov E. I., 2023, MNRAS, 518, 791

* Eldridge J. J., Genet F., Daigne F., Mochkovitch R., 2006, MNRAS, 367, 186

* Fossati L. et al., 2015, A&A, 574, A20

* Freyer T., Hensler G., Yorke H. W., 2003, ApJ, 594, 888

* Freyer T., Hensler G., Yorke H. W., 2006, ApJ, 638, 262

* Gallant Y. A., Arons J., 1994, ApJ, 435, 230

* Gull T. R., Sofia S., 1979, ApJ, 230, 782

* Henney W. J., Arthur S. J., de Colle F., Mellema G., 2009, MNRAS, 398, 157

* Herbst K. et al., 2020, ApJ, 897, L27

* Hummer D. G., 1994, MNRAS, 268, 109

* Igoshev A. P., 2020, MNRAS, 494, 3663

* Kargaltsev O., Pavlov G. G., 2010, in Comastri A., Angelini L., Cappi M., eds, AIP Conf. Proc. Vol. 1248, X-ray Astronomy 2009; Present Status, Multi-Wavelength Approach and Future Perspectives. Am. Inst. Phys., New York, p. 25
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 PWNe 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. The interaction between the supernova ejecta and the pulsar wind with this complex environment gives rise to the observed irregular morphologies.

In future research, the numerical model presented in this study will be refined and customized to specific supernova remnants hosting PWNe, such as the case of the old plerion G106.3+2.7 and its central pulsar J2229+6114, which drives the Boomerang synchrotron shell. By tailoring the model to these specific systems, researchers aim to achieve a more detailed and accurate understanding of the morphologies and dynamics of these objects.

## ACKNOWLEDGEMENTS

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 grant PID2021-124581OB-I00 funded by MCIN/AEI/10.13039/501100011033, 2021SGR00426 of the Generalitat de Catalunya, by the programme Unidad de Excelencia María de Maestu CEX2020-001058-M, and by MCIN with funding from European Union NextGeneration EU (PRTR-C17.11).

## 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.itai.uni-heidelberg.de/~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

* Aharonian F. et al., 2006, A&A, 457, 899

* Arias M. et al., 2019, A&A, 622, A6

* Aschenbach B., Leahy D. A., 1999, A&A, 341, 602

* Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

* Baalmann L. R., Scherer K., Fichtner H., Kleimann J., Bomans D. J., Weis K., 2020, A&A, 634, A67

* Baalmann L. R., Scherer K., Kleimann J., Fichtner H., Bomans D. J., Weis K., 2021, A&A, 650, A36

* Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2020, MNRAS, 499, 2051

* Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2021, MNRAS, 508, 3194

* Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 2023, MNRAS, 520, 2451

* Barkov M. V., Lyutikov M., Khangulyan D., 2019a, MNRAS, 484, 4760

* Barkov M. V., Lyutikov M., Klingler N., Bordas P., 2019b, MNRAS, 485, 2041

* Barkov M. V., Lyutikov M., Khangulyan D., 2020, MNRAS, 497, 2605

* Begelman M. C., 1998, ApJ, 493, 291

* Begelman M. C., Li Z.-Y., 1992, ApJ, 397, 187

* Blondin J. M., Chevalier R. A., 2017, ApJ, 845, 139

* Blondin J. M., Chevalier R. A., Frierson D. M., 2001, ApJ, 563, 806

* Bock D. C. J., Turtle A. J., Green A. J., 1998, AJ, 116, 1886

* Bucciantini N., 2002a, A&A, 387, 1066

* Bucciantini N., 2002b, A&A, 393, 629

* Bucciantini N., 2018, MNRAS, 478, 2074

* Bucciantini N., Bandiera R., 2001, A&A, 375, 1032

* Bucciantini N., Bandiera R., Blondin J. M., Amato E., Del Zanna L., 2004, A&A, 422, 609

* Burns R. A. et al., 2023, Nat. Astron., 7, 557

* Castro N. et al., 2015, A&A, 581, A81

* Castro N. et al., 2017, A&A, 597, L6

* Chevalier R. A., 1982, ApJ, 258, 790

* Chevalier R. A., Luo D., 1994, ApJ, 421, 225

* Chiotellis A., Boumis P., Derlopa S., Steffen W., 2019, Proceedings of the conference "Supernova Remnants II: An Odyssey in Space after Stellar Death". preprint (arXiv:1909.08947)

* Comerón F., Kaper L., 1998, A&A, 338, 273

* Coroniti F. V., 1990, ApJ, 349, 538

* Cox D. P., Shelton R. L., Maciejewski W., Smith R. K., Plewa T., Pawl A., Różyczka M., 1999, ApJ, 524, 179

* Cox N. L. J., Kerschbaum F., van Marle A. J., Decin L., Ladjal D., Mayer A., 2012, A&A, 543, C1

* Crawford F., Gaensler B. M., Kaspi V. M., Manchester R. N., Camilo F., Lyne A. G., Pivovaroff M. J., 2001, ApJ, 554, 152

* Das S., Brose R., Meyer D. M. A., Pohl M., Sushch I., Plotko P., 2022, A&A, 661, A128

* Decin L., 2012, Adv. Space Res., 50, 843

* Decin L. N. L. J. et al., 2012, A&A, 548, A113

* De Luca A. et al., 2013, ApJ, 765, L19

* Derlopa S., Boumis P., Chiotellis A., Steffen W., Akras S., 2020, MNRAS, 499, 5410

* de Vries M. et al., 2021, ApJ, 908, 50

* Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium. Princeton Univ. Press, Princeton, NJ

* Dwarkadas V. V., 2005, ApJ, 630, 892

* Dwarkadas V. V., 2007, ApJ, 667, 226

* Ekström S. et al., 2012, A&A, 537, A146

* Elbakyan V. G., Nayakshin S., Meyer D. M. A., Vorobyov E. I., 2023, MNRAS, 518, 791

* Eldridge J. J., Genet F., Daigne F., Mochkovitch R., 2006, MNRAS, 367, 186

* Fossati L. et al., 2015, A&A, 574, A20

* Freyer T., Hensler G., Yorke H. W., 2003, ApJ, 594, 888

* Freyer T., Hensler G., Yorke H. W., 2006, ApJ, 638, 262

* Gallant Y. A., Arons J., 1994, ApJ, 435, 230

* Gull T. R., Sofia S., 1979, ApJ, 230, 782

* Henney W. J., Arthur S. J., de Colle F., Mellema G., 2009, MNRAS, 398, 157

* Herbst K. et al., 2020, ApJ, 897, L27

* Hummer D. G., 1994, MNRAS, 268, 109

* Igoshev A. P., 2020, MNRAS, 494, 3663

* Kargaltsev O., Pavlov G. G., 2010, in Comastri A., Angelini L., Cappi M., eds, AIP Conf. Proc. Vol. 1248, X-ray Astronomy 2009; Present Status, Multi-Wavelength Approach and Future Perspectives. Am. Inst. Phys., New York, p. 25
Kargaltsev O., Schmitt B. M., Pavlov G. G., Misanovic Z., 2012, ApJ, 745, 99

Katsuda S., Takiwaki T., Tominaga N., Moriya T. J., Nakamura K., 2018, ApJ, 863, 127

Kennel C. F., Coroniti F. V., 1984, ApJ, 283, 710

Kervella P. et al., 2018, A&A, 609, A67

Klingler N., Rangelov B., Kargaltsev O., Pavlov G. G., Romani R. W., Slane P. O., XVP PWN Collaboration, 2014, AAS Meeting Abstracts, #223

Kolb C., Blondin J., Slane P., Temim T., 2017, ApJ, 844, 1

Komissarov S. S., 2006, MNRAS, 367, 19

Komissarov S. S., Lyubarsky Y. E., 2003, MNRAS, 344, L93

Komissarov S. S., Lyubarsky Y. E., 2004, MNRAS, 349, 779

Komissarov S. S., Lyutikov M., 2011, MNRAS, 414, 2017

Klothes R., Uyaniker B., Pineault S., 2001, ApJ, 560, 236

Klothes R., Reich W., Uyaniker B., 2006, ApJ, 638, 225

Lasky P. D., 2015, Publ. Astron. Soc. Aust., 32, e034

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 H., Safi-Harb S., Kothes R., 2013, ApJ, 774, 33

Mestel L., Robertson J. A., Wang Y. M., Westfold K. C., 1985, MNRAS, 217, 443

Meyer D. M. A., 2021, MNRAS, 507, 4697

Meyer D. M. A., Meliani Z., 2022, MNRAS, 515, L29

Meyer D. M.-A., Mackey J., Langer N., Gvaramadze V. V., Mignone A., Izzard R. G., Kaper L., 2014, MNRAS, 444, 2754

Meyer D. M.-A., Langer N., Mackey J., Velázquez P. F., Gusdorf A., 2015, MNRAS, 450, 3080

Meyer D. M. A., Petrov M., Pohl M., 2020, MNRAS, 493, 3548

Meyer D. M. A., Mignone A., Petrov M., Scherer K., Velázquez P. F., Boumis P., 2021, MNRAS, 506, 5170

Meyer D. M. A. et al., 2022a, MNRAS, 515, 594

Meyer D. M. A., Vorobiyov E. I., Elbakyan V. G., Kraus S., Liu S. Y., Nayakshin S., Sobolev A. M., 2022b, MNRAS, 517, 4795

Meyer D. M.-A., Pohl M., Petrov M., Egberts K., 2023, MNRAS, 521, 5354

Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228

Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G., 2012, ApJS, 198, 7

Mignone A., Striani E., Tavani M., Ferrari A., 2013, MNRAS, 436, 1102

Olmi B., Bucciantini N., 2019, MNRAS, 488, 5690

Olmi B., Bucciantini N., 2023, Publ. Astron. Soc. Aust., 40, e007

Osterbrock D. E., Bochkarev N. G., 1989, SvA, 33, 694

Parker E. N., 1958, ApJ, 128, 664

Pétri J., 2022, A&A, 666, A5

Pineault S., Joncas G., 2000, AJ, 120, 3218

PNigouras P., 2019, Phys. Rev. D, 100, 063016

PNigouras P., Kokkotas K. D., 2015, Phys. Rev. D, 92, 084018

PNigouras P., Kokkotas K. D., 2016, Phys. Rev. D, 94, 024053

Pogorelov N. V., Matsuda T., 2000, A&A, 354, 697

Pogorelov N. V., Semenov A. Y., 1997, A&A, 321, 330

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, Astron. Rep., 63, 391

Porth O., Komissarov S. S., Keppens R., 2014, MNRAS, 438, 278

Posselt B. et al., 2017, ApJ, 835, 66

Powell K. G., 1997, in Hussaini M.Y., van Leer B., Van Rosendale J., eds, An Approximate Riemann Solver for Magnetohydrodynamics. Springer-Verlag, Berlin, p. 570

Przybilla N. et al., 2016, A&A, 587, A7

Rozyczka M., Franco J., 1996, ApJ, 469, L127

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 P., 2017, in Alsabti A. W., Murdin P., eds, Handbook of Supernovae. Springer, p. 2159

Steiner A. W., Prakash M., Lattimer J. M., Ellis P. J., 2005, Phys. Rep., 411, 325

Stock D. J., Barlow M. J., 2010, MNRAS, 409, 1429

Sun M., Wang Z.-r., Chen Y., 1999, ApJ, 511, 274

Temim T., Slane P., Plucinsky P. P., Gelfand J., Castro D., Kolb C., 2017, ApJ, 851, 128

Temim T., Slane P., Raymond J. C., Patnaude D., Murray E., Ghavamian P., Renzo M., Jacovich T., 2022, ApJ, 932, 26

Torres D. F., Lin T., 2018, ApJ, 864, L2

Truelove J. K., McKee C. F., 1999, ApJS, 120, 299

Vaidya B., Mignone A., Bodo G., Rossi P., Massaglia S., 2018, ApJ, 865, 144

van der Swaluw E., 2003, A&A, 404, 939

van der Swaluw E., Achterberg A., Gallant Y. A., Tóth G., 2001, A&A, 380, 309

van der Swaluw E., Achterberg A., Gallant Y. A., Downes T. P., Keppens R., 2003, A&A, 397, 913

van der Swaluw E., Downes T. P., Keegan R., 2004, A&A, 420, 937

van Marle A. J., Meliani Z., Marcowith A., 2015a, A&A, 584, A49

van Marle A. J., Meliani Z., Marcowith A., 2015b, A&A, 584, A49

van Veelen B., Langer N., Vink J., García-Segura G., van Marle A. J., 2009, A&A, 503, 495

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

Vlemmings W. H. T., Diamond P. J., van Langevelde H. J., 2002, A&A, 394, 589

Vlemmings W. H. T., van Langevelde H. J., Diamond P. J., 2005, A&A, 434, 1029

Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377

Weber F., 1999, J. Phys. G: Nucl. Phys., 25, R195

Weber E. J., Davis Leverett J., 1967, ApJ, 148, 217

Whalen D., van Veelen B., O'Shea B. W., Norman M. L., 2008, ApJ, 682, 49

Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99

Wilkin F. P., 1996, ApJ, 459, L31

Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ, 587, 278

Yusef-Zadeh F., Wardle M., Rho J., Sakano M., 2003, ApJ, 585, 319

This paper has been typeset from a *TEX/LATEX* file prepared by the author.
Kargaltsev O., Schmitt B. M., Pavlov G. G., Misanovic Z., 2012, ApJ, 745, 99

Katsuda S., Takiwaki T., Tominaga N., Moriya T. J., Nakamura K., 2018, ApJ, 863, 127

Kennel C. F., Coroniti F. V., 1984, ApJ, 283, 710

Kervella P. et al., 2018, A&A, 609, A67

Klingler N., Rangelov B., Kargaltsev O., Pavlov G. G., Romani R. W., Slane P. O., XVP PWN Collaboration, 2014, AAS Meeting Abstracts, #223

Kolb C., Blondin J., Slane P., Temim T., 2017, ApJ, 844, 1

Komissarov S. S., 2006, MNRAS, 367, 19

Komissarov S. S., Lyubarsky Y. E., 2003, MNRAS, 344, L93

Komissarov S. S., Lyubarsky Y. E., 2004, MNRAS, 349, 779

Komissarov S. S., Lyutikov M., 2011, MNRAS, 414, 2017

Klothes R., Uyaniker B., Pineault S., 2001, ApJ, 560, 236

Klothes R., Reich W., Uyaniker B., 2006, ApJ, 638, 225

Lasky P. D., 2015, Publ. Astron. Soc. Aust., 32, e034

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 H., Safi-Harb S., Kothes R., 2013, ApJ, 774, 33

Mestel L., Robertson J. A., Wang Y. M., Westfold K. C., 1985, MNRAS, 217, 443

Meyer D. M. A., 2021, MNRAS, 507, 4697

Meyer D. M. A., Meliani Z., 2022, MNRAS, 515, L29

Meyer D. M.-A., Mackey J., Langer N., Gvaramadze V. V., Mignone A., Izzard R. G., Kaper L., 2014, MNRAS, 444, 2754

Meyer D. M.-A., Langer N., Mackey J., Velázquez P. F., Gusdorf A., 2015, MNRAS, 450, 3080

Meyer D. M. A., Petrov M., Pohl M., 2020, MNRAS, 493, 3548

Meyer D. M. A., Mignone A., Petrov M., Scherer K., Velázquez P. F., Boumis P., 2021, MNRAS, 506, 5170

Meyer D. M. A. et al., 2022a, MNRAS, 515, 594

Meyer D. M. A., Vorobiyov E. I., Elbakyan V. G., Kraus S., Liu S. Y., Nayakshin S., Sobolev A. M., 2022b, MNRAS, 517, 4795

Meyer D. M.-A., Pohl M., Petrov M., Egberts K., 2023, MNRAS, 521, 5354

Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228

Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G., 2012, ApJS, 198, 7

Mignone A., Striani E., Tavani M., Ferrari A., 2013, MNRAS, 436, 1102

Olmi B., Bucciantini N., 2019, MNRAS, 488, 5690

Olmi B., Bucciantini N., 2023, Publ. Astron. Soc. Aust., 40, e007

Osterbrock D. E., Bochkarev N. G., 1989, SvA, 33, 694

Parker E. N., 1958, ApJ, 128, 664

Pétri J., 2022, A&A, 666, A5

Pineault S., Joncas G., 2000, AJ, 120, 3218

PNigouras P., 2019, Phys. Rev. D, 100, 063016

PNigouras P., Kokkotas K. D., 2015, Phys. Rev. D, 92, 084018

PNigouras P., Kokkotas K. D., 2016, Phys. Rev. D, 94, 024053

Pogorelov N. V., Matsuda T., 2000, A&A, 354, 697

Pogorelov N. V., Semenov A. Y., 1997, A&A, 321, 330

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, Astron. Rep., 63, 391

Porth O., Komissarov S. S., Keppens R., 2014, MNRAS, 438, 278

Posselt B. et al., 2017, ApJ, 835, 66

Powell K. G., 1997, in Hussaini M.Y., van Leer B., Van Rosendale J., eds, An Approximate Riemann Solver for Magnetohydrodynamics. Springer-Verlag, Berlin, p. 570

Przybilla N. et al., 2016, A&A, 587, A7

Rozyczka M., Franco J., 1996, ApJ, 469, L127

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 P., 2017, in Alsabti A. W., Murdin P., eds, Handbook of Supernovae. Springer, p. 2159

Steiner A. W., Prakash M., Lattimer J. M., Ellis P. J., 2005, Phys. Rep., 411, 325

Stock D. J., Barlow M. J., 2010, MNRAS, 409, 1429

Sun M., Wang Z.-r., Chen Y., 1999, ApJ, 511, 274

Temim T., Slane P., Plucinsky P. P., Gelfand J., Castro D., Kolb C., 2017, ApJ, 851, 128

Temim T., Slane P., Raymond J. C., Patnaude D., Murray E., Ghavamian P., Renzo M., Jacovich T., 2022, ApJ, 932, 26

Torres D. F., Lin T., 2018, ApJ, 864, L2

Truelove J. K., McKee C. F., 1999, ApJS, 120, 299

Vaidya B., Mignone A., Bodo G., Rossi P., Massaglia S., 2018, ApJ, 865, 144

van der Swaluw E., 2003, A&A, 404, 939

van der Swaluw E., Achterberg A., Gallant Y. A., Tóth G., 2001, A&A, 380, 309

van der Swaluw E., Achterberg A., Gallant Y. A., Downes T. P., Keppens R., 2003, A&A, 397, 913

van der Swaluw E., Downes T. P., Keegan R., 2004, A&A, 420, 937

van Marle A. J., Meliani Z., Marcowith A., 2015a, A&A, 584, A49

van Marle A. J., Meliani Z., Marcowith A., 2015b, A&A, 584, A49

van Veelen B., Langer N., Vink J., García-Segura G., van Marle A. J., 2009, A&A, 503, 495

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

Vlemmings W. H. T., Diamond P. J., van Langevelde H. J., 2002, A&A, 394, 589

Vlemmings W. H. T., van Langevelde H. J., Diamond P. J., 2005, A&A, 434, 1029

Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377

Weber F., 1999, J. Phys. G: Nucl. Phys., 25, R195

Weber E. J., Davis Leverett J., 1967, ApJ, 148, 217

Whalen D., van Veelen B., O'Shea B. W., Norman M. L., 2008, ApJ, 682, 49

Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99

Wilkin F. P., 1996, ApJ, 459, L31

Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ, 587, 278

Yusef-Zadeh F., Wardle M., Rho J., Sakano M., 2003, ApJ, 585, 319

This paper has been typeset from a *TEX/LATEX* file prepared by the author.