On the plerionic rectangular supernova remnants of static progenitors

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

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 phe-nomena, fluid dynamics, general relativity, and nuclear physics, see, e.g., Weber (1999); Bucciantini (2011); Steiner et al. (2005); Lasky (2015); Pnigouras & Kokkotas (2015, 2016); Pnigouras (2019).Pulsars have very powerful magnetospheres with strong magnetic fields on the order of kilogauss (kG), which play a crucial role in the evolution see e.g., Mestel et al. (1985).The rotating magnetospheres extract energy from the pulsar and generate a powerful wind, see, e.g., Pétri (2022).The interaction of the pulsar wind with the ambient medium produces the so-called pulsar nebulae, which can be located inside or outside the supernova remnant of the progenitor star, depending on whether the supernova explosion had kicked off the pulsar.Various observations have triggered investigations into such phenomena, see, amongst others, the studies of Pavan et al. (2016); Kargaltsev et al. (2017); de Vries & Romani (2020); Igoshev (2020); de Vries et al. (2021).
A particular class of supernova remnants containing a pulsar exhibit a succession of structured shocks powered by the pulsar's magnetic wind, producing multi-wavelength polarized non-thermal emission.Examples of such plerions include the Crab Nebula (Hester 2008), as well as the Geminga pulsar residing within the Vela supernova remnant (Bock et al. 1998;Popov et al. 2019).Additionally, one can observe youthful supernova remnants hosting both a pulsar and a pulsar-wind nebula, such as B0540-693 (Williams et al. 2008) and G11.2-0.3 (Borkowski et al. 2016).
The modelling of PWN has been a long-standing challenge for several reasons.First, the physics involved in PWN is inherently complex, involving the interaction between the pulsar's relativistic wind and the surrounding medium.This requires a multi-disciplinary approach.Second, the environment in which the pulsar wind is launched is often structured, as it depends on the supernova remnant's properties, the progenitor star's circumstellar medium, or the interstellar medium in which the pulsar resides.The properties of the surrounding medium can significantly affect the dynamics and emission of the PWN.These factors together make the modelling of PWN a complex and multi-faceted problem, requiring sophisticated theoretical models and numerical simulations to understand the physics at play fully.The Crab Nebula stands out as a prominent example of a plerion.Extensive research, both theoretical (Kennel & Coroniti 1984;Coroniti 1990;Begelman & Li 1992;Begelman 1998) and numerical, has been dedicated to studying PWN like the Crab Nebula.These investigations encompass relativistic axisymmetric 2D simulations (e.g., Komissarov & Lyubarsky 2003, 2004;Komissarov 2006;Del Zanna et al. 2006;Camus et al. 2009;Komissarov & Lyutikov 2011;Olmi et al. 2014) as well as relativistic 3D simulations (e.g., Mignone et al. 2013;Porth et al. 2014;Olmi et al. 2016).
Moreover, the Crab Nebula is pivotal in advancing our understanding of pulsar physics and their interactions with supernova remnants.Notably, the dynamics and morphology of pulsar wind nebulae experience significant transformations as they expand within a supernova remnant.This influence can be even more pronounced when the pulsar receives a kick during a supernova explosion, as observed in the case of PWN CTB 87 (Matheson et al. 2013).Extensive studies have been conducted on the interaction between PWN and supernova remnants, focusing on no-moving pulsars in 1D and 2D scenarios (van der Swaluw et al. 2001;Blondin & Chevalier 2017), including a mock complex surrounding the neutron star (van der Swaluw 2003; Blondin et al. 2001).
These studies have also been extended to moving pulsars inside the supernova ejecta, revealing the development of strong bow-shocks between the pulsar wind and supernova remnant (van der Swaluw et al. 2003(van der Swaluw et al. , 2004;;Temim et al. 2017;Kolb et al. 2017;Temim et al. 2022) and the study was extended using relativistic MHD (Bucciantini et al. 2004).In some cases, the strong interaction between the PWN and the reverse shock of the supernova remnant can result in compression.This interaction phase is referred to as reverberation, during and after which the morphology, spectrum, and dynamics of PWN could undergo significant changes, see the recent studies by Torres & Lin (2018); Bandiera et al. (2020Bandiera et al. ( , 2023a,b),b).
In later phases, as the moving pul sar leaves the supernova remnant and begins to interact with the interstellar medium, a bow shock nebula forms around the runaway pulsar.This intriguing phenomenon has been extensively studied analytically (Bucciantini & Bandiera 2001) and through numerical simulations in two dimensions with relativistic considerations (Bucciantini 2018).Furthermore, research into this phenomenon has delved into three-dimensional simulations, encompassing non-relativistic pulsar winds (Toropina et al. 2019) and relativistic pulsar winds (Barkov et al. 2019a(Barkov et al. ,b, 2020;;Olmi & Bucciantini 2019).Despite the high complexity of these simulations and the numerous questions they leave unanswered (Olmi & Bucciantini 2023), all of these studies still neglect the effects of the circumstellar medium of the defunct star in which the supernova remnant and the PWN expand during the initial phases.
The circumstellar medium of a defunct star is formed through the interaction between the star's wind and luminosity with the sur-rounding interstellar medium (ISM).The shape and properties of the circumstellar medium depend on various factors, including the evolution of the star, such as its mass, age, and stage of evolution, as well as the characteristics of the surrounding ISM, such as density, temperature, and magnetic field van Marle et al. (2015a); Meyer et al. (2022a).In the context of massive stars, the circumstellar medium undergoes successive structural changes.During its early life, it forms an accretion disc (Liu et al. 2020;Meyer et al. 2022b;Elbakyan et al. 2023;Burns et al. 2023).In its main sequence phase, it expands into a wind bubble (Weaver et al. 1977;Gull & Sofia 1979;Wilkin 1996a).Later on, it evolves into expanding shells (Stock & Barlow 2010;Cox et al. 2012;Decin 2012;Decin et al. 2012).If a supernova explosion occurs, it leaves behind an expanding remnant shell (Aschenbach & Leahy 1999;Yusef-Zadeh et al. 2003;Katsuda et al. 2018;Arias et al. 2019;Chiotellis et al. 2019;Derlopa et al. 2019).
Once the pulsar emits a relativistic and powerful wind, it initially interacts with the surrounding supernova ejecta (Cox et al. 1999;Sun et al. 1999;Crawford et al. 2001;Olmi & Bucciantini 2023).As the PWN passes through the supernova ejecta, it subsequently interacts with the circumstellar medium of the defunct star.The distribution of ejecta, stellar wind, and ISM gas acts as a matrix that channels the expansion of the pulsar wind (Kolb et al. 2017;Temim et al. 2022).This is particularly important when the supernova progenitor is a runaway star, as the bow shock created by its surrounding stellar wind can influence the subsequent evolution and emission of the supernova ejecta and the PWN (Meyer & Meliani 2022).
This study aims to investigate how a magnetised ambient medium influences the dynamics, morphologies, and emission properties of PWN with static massive star progenitors.The multi-dimensional magnetohydrodynamical (MHD) simulations conducted by van Marle et al. (2015b) have revealed that the circumstellar medium of high-mass stars is significantly influenced by the organized magnetic field of its ambient medium.This finding has profound implications on the understanding of stellar wind bubbles around massive stars, as previously studied by Freyer et al. (2003); Dwarkadas (2005); Freyer et al. (2006); Dwarkadas (2007).The presence of a magnetic field can cause expanding stellar wind bubbles to become elongated and adopt an oblong morphology along the direction of the magnetic field lines.Our previous work (Meyer et al. 2022a) has shown that such asymmetric pre-supernova environments can result in a peculiar reflection of the supernova shock wave, forming rectangular-shaped remnants like Puppis A. In this study, we further investigate the effects of the reflection of the supernova blastwave in asymmetric, magnetized wind bubbles that are generated by a static, rotating star in the warm phase of the Galactic plane and how it may impact the evolution of plerionic pulsar wind nebulae.
The paper is structured as follows.In Section 2, we present the modelling methods used in this study.This includes the description of the numerical simulations of pulsar wind nebulae, which are detailed in Section 3. We then discuss the outcomes of our study in Section 4, and present our conclusions in Section 5.

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

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  ISM = 0.79 cm −3 , while the magnetic field of the ISM is uniform and structured, with a strength of  ISM = 7 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 et al. (2009), which is suitable for a solar metallicity environment (Asplund et al. 2009).The cooling law accounts for hydrogen and helium as the primary coolants at temperatures  < 10 6 K, and various metals' emission lines at temperatures  ⩾ 10 6 K.The cooling curve is further enhanced with [Oiii] , 5007 collisionally excited forbidden lines, as described in Asplund et al. (2009) and Henney et al. (2009).This paper presents a model that captures the evolution of the circumstellar medium surrounding a static massive star with an initial mass of 35 M ⊙ at the zero-age main sequence.The star is considered to be rotating with an angular velocity ratio of Ω ★ /Ω K = 0.1, where Ω ★ represents the star's initial angular frequency and Ω K is its equatorial Keplerian angular velocity.Consequently, the equatorial velocity of the star can be expressed as,  rot () = Ω ★ () ★ (). (1) Here,  ★ () 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 onedimensional stellar evolution model provided by the geneva library, as described in Ekström et al. (2012) 1 .Specifically, we extract the mass-loss rate  () and the effective temperature  eff () of the star at each stage of evolution from this database, and derive from it the wind density, In this equation,  represents the radial distance from the star, and  () corresponds to the mass loss rate of the star at time .
1 https://www.unige.ch/sciences/astro/evolution/en/database/syclist/ The terminal velocity of the stellar wind, denoted as  w (), is calculated based on the escape velocity  esc ().The escape velocity depends on the star's effective temperature  eff and is determined using the conversion law, where  represents the gravitational constant, and () is a the normalisation factor introduced by Eldridge et al. (2006).
We adopt the time-dependent evolution of the surface magnetic field  ★ of the supernova progenitor as derived in Meyer et al. (2023) where the magnetic field strength at the surface of the star are scaled to that of the Sun, as described in Scherer et al. (2020); Herbst et al. (2020);Baalmann et al. (2020Baalmann et al. ( , 2021)); Meyer et al. (2021b).Specifically, we assume a magnetic field strength at star surface of  ★ = 500 G during the main-sequence phase (Fossati et al. 2015;Castro et al. 2015;Przybilla et al. 2016;Castro et al. 2017), to a Betelgeuse-like field of  ★ = 0.2 G for the red supergiant phase (Vlemmings et al. 2002(Vlemmings et al. , 2005;;Kervella et al. 2018) and  ★ = 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, and a toroidal component, respectively, with, being the latitude-dependent surface velocity of the rotating massive star (Parker 1958;Weber & Davis 1967;Pogorelov & Semenov 1997;Pogorelov & Matsuda 2000;Chevalier & Luo 1994;Rozyczka & Franco 1996).
At the end of a star's evolution, it enters the supernova phase, during which we model the expanding supernova ejecta as a spherically symmetric distribution within a radius of  max .The ejecta has a total energy of  SN = 10 51 erg and a mass of  SN , which takes into account the star's mass loss throughout its entire evolution until the immediate pre-supernova time  psn , as well as the mass  NS of the neutron star that forms at the centre.Specifically, we set, with  psn and  NS = 1.4 M ⊙ (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  core , where  core represents the core radius.In this region, the density decreases with time following a power-law relationship of  −3 , where  denotes the time after the explosion.The second region is the outer edge, extending from  core to  max , where  max corresponds to the maximum radius.In this region, the density decreases steeply with radius, following a power-law relationship of  ∝  − , with the exponent  set to 11.Additionally, the density in the outer edge region decreases with time as  − (3+) .These density profiles can be expressed as follows: and, respectively.These density profiles are commonly used for corecollapse supernovae (Chevalier 1982).
For the velocity, we utilize a homologous radial profile for the supernova ejecta, given by  = /, across all regions from 0 to  max at time  max .The characteristics of the supernova ejecta profile are computed following the methodology outlined in Truelove & McKee (1999) and Whalen et al. (2008).
The velocity at the core radius, denoted as  core ( core ), is determined as, where  sn represents the total energy of the supernova ejecta and  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  max , is set to: This choice of  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, , according to the following equation: with  o is the initial spin-down of the pulsar, defined as, where the initial spin period of the pulsar is set to  o = 0.3 s, and its time derivative is set to  o = 10 −17 s s −1 .The braking index is assumed to be  = 3, which corresponds to magnetic dipole spindown, as outlined in Pacini (1967).Furthermore, we assume that the pulsar's wind maintains a constant speed of  psw = 10 −2 , where  denotes the speed of light in vacuum.It's important to acknowledge that this speed is significantly lower than the realistic pulsar wind speeds, which can reach , corresponding to a Lorentz factor of 10 6 (as demonstrated in Kennel & Coroniti 1984).This decision to employ a reduced pulsar wind speed can lead to noticeable alterations in the properties of its termination shock.These changes encompass compression rates, speeds, and, subsequently, shock positions and influence the development of associated instabilities.It is crucial to emphasize that our paper's primary objective is to replicate the overall evolution of the PWN accurately.This evolution is predominantly governed by the wind's momentum flux (Wilkin 1996b).In terms of magnetization, we have opted for a low value of  = 10 −3 in this study, a choice in line with descriptions found in Rees & Gunn (1974), Kennel & Coroniti (1984), Slane (2017), Begelman & Li (1992) and Torres et al. (2014).This magnetization value implies that a significant portion of the magnetic field is converted into kinetic energy.It is worth noting that recent multi-dimensional simulations have demonstrated that larger magnetization values, such as  = 0.01 in 2D (e.g., Komissarov & Lyubarsky 2003, 2004;Del Zanna et al. 2004, 2006) and even  > 1 in 3D (Porth et al. 2014;Barkov et al. 2019a), can accurately reproduce the features of termination shocks of PWN.
However, it is essential to recognize that the value of pulsar wind magnetization remains a topic of debate, as it significantly influences PWN termination shock strength and, consequently, particle acceleration.Moreover, the magnetization of the equatorial wind zone may decrease, leading to lower magnetization values due to the annihilation of equatorial wind magnetic stripes (Coroniti 2017).By selecting such a low magnetization, as in Bucciantini et al. (2004), the Pulsar Wind Nebula tends to expand more in the equatorial plane, resulting in a stronger termination shock.
Furthermore, Komissarov and Lyubarsky in 2003-2004and Del Zanna et al. 2004, showed that the properties of the inner nebula can only be recovered with 2D simulations if the injected magnetization is larger than 0.01.
The magnetic field is assumed to have only a toroidal component.The total kinetic energy, magnetic field strength, and kinetic energy are functions of the radial distance  and polar angle , as described in Komissarov & Lyubarsky (2004).
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.

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. (2015Meyer et al. ( , 2020) ) and that we extended after to PWN in Meyer & Meliani (2022).The magneto-hydrodynamical simulations are conducted in a 2.5-dimensional, axisymmetric cylindrical coordinate system.The simulation box extends over the range [;  max ] × [ min ;  max ] and is discretized using a uniform grid of  R ×  z cells.Consequently, the spatial resolution is consistent along both directions, with each grid cell having a size of Δ =  max / R .We employ two different spatial resolutions throughout the evolutionary process.During the progenitor star wind phases, the circumstellar medium is resolved using a grid resolution of  R = 2000 and  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Δ, 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  R = 3000 and  z = 6000 cells.The supernova ejecta is confined within a central sphere of radius  max , as described in section 2.1.Simultaneously, the pulsar wind is imposed within a sphere of radius  ns wind = 20Δ, 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 ⊙ 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-PWNand Run-35-MHD-0.79-PWN.In the magnetised external medium case, the adopted strength of the background magnetic field is set to that measured in the spiral arms of the Galaxy, with an average strength of  ISM = 7 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(Mignone et al. , 2012;;Vaidya et al. 2018)2 and we solve the following set of equations, and with the gas density , velocity , momentum  =  and magnetic field , as well as the the total pressure   and the energy of the gas The sound speed of the medium reads where the adiabatic index is  = 5/3.Last, radiative cooling by optically-thin processes and photo-heating are included into the equations via the term Φ(, ), with the gas temperature , accounting for the prescriptions of Meyer et al. (2014).Regarding the cooling/heating processes of the gas, we assume the gas to be optically thin throughout the entire progenitor's life.After this point, with the launch of the pulsar wind, the cooling and heating mechanisms are disabled.We employ a Godunov-type numerical scheme with the Harten-Lax-van Leer approximate Riemann solver (HLL) and utilize the 8-waves magnetic field formulation (Powell 1997).For time integration, a third-order Runge-Kutta scheme is employed, with the time step controlled by the Courant-Friedrichs-Lewy (CFL) number.The numerical simulations are performed at the North-German Supercomputing Alliance (HLRN3 ) using the LISE cluster in Berlin, which is equipped with Cray XC40/30 processors.

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

Model with unmagnetised ISM
In Figure 1, the density contour is shown for the unmagnetized case Run-35-HD-0.79-PWN(left panels) and the magnetized case Run-35-MHD-0.79-PWN(right panels) at different evolution times, from to top to the bottom.The density contour is represented in logarithmic scale in cm −3 units.In both cases, the red contour marks indicate the region of the plerionic supernova remnants where the contribution of the pulsar wind reaches 50 times of the number density.
In Figure 1a (top-left), we present the pre-supernova circumstellar medium.At this stage, it forms a large-scale quasi-spherical stellar bubble (Weaver et al. 1977), and its spherical forward shock extends to distances of approximately 90 pc.Throughout the star's evolution, the stellar wind interacts strongly with the ambient medium.Each phase of evolution contributes to the formation of successive shock structures, which appear in order from the farthest to the nearest region to the star.The thick and dense shell located farthest from the star, with a radial extent of ⩾ 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 et al. 2003;Dwarkadas 2005;Freyer et al. 2006;Dwarkadas 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 due to the interaction between the red-giant wind and the main-sequence wind.Additionally, the main-sequence wind interacts with the surrounding ambient medium, forming an external dense shell that is limited by the contact discontinuity surface.
It is worth noting that the contact discontinuity, which marks the interface between the wind and the ISM, exhibits a slightly aspherical morphology, particularly in the region close to the symmetry axis.This aspherical shape is influenced by the presence of the magnetic field and the star's rotation.The variations between the bubbles depicted in Fig. 1a and Fig. 1 of Meyer et al. (2022a) highlight this effect.Furthermore, the grid's proximity to the near-polar axis amplifies this asymmetry.Moving on to Fig. 1c (middle-left), we can observe the supernova remnant at 25 kyr after the explosion.The expanding shock wave from the supernova remnant propagates outward, sweeping up and pushing away all the previously formed dense shells associated with the successive stellar winds.As the shock wave reaches the contact discontinuity surface between the main-sequence stellar wind and the ISM, it interacts with this surface, causing reflection, as described in Meyer et al. (2015Meyer et al. ( , 2021a)); Meyer & Meliani (2022).This interaction and reverberation contribute to the observed structure and morphology of the supernova remnant.
After the supernova explosion, a pulsar wind with high initial mechanical luminosity,  psr,0 = 10 38 ergs −1 , is launched.However, this luminosity decreases over time according to Eq. 21 of Pacini (1967).This pulsar wind interacts with the dense supernova ejecta (van der Swaluw et al. 2004), resulting in the formation of a complex structure as described in Meyer & Meliani (2022) for a runaway progenitor star with a zero-age main-mass of 20 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 interaction of the pulsar wind with the expanding supernova remnant.A pulsar wind termination shock is formed at the interface between the unperturbed pulsar wind and the shocked pulsar wind.The outermost region of the pulsar wind nebula behind the termination shock contains the contact discontinuity.This contact discontinuity marks the interface between the supernova ejecta and the shocked pulsar wind (depicted by the red contour in Fig. 1).Beyond the contact discontinuity, a transmitted pulsar wind forward shock propagates through the still unshocked supernova ejecta and further travels into the surrounding medium.
The pulsar wind contact discontinuity undergoes expansion to larger radii due to the fast rotation of the magnetized neutron star.This expansion leads to the characteristic shape with an equatorial torus and an elongated polar jet, as found by Komissarov & Lyubarsky (2004); Del Zanna et al. (2006); Porth et al. (2014); Olmi et al. (2016).However, it's important to note that due to limitations in the numerical scheme applied to the 2D symmetry axis, the jet along the polar axis may appear more elongated than it would in a full 3D simulation.Nevertheless, despite these limitations, the general behaviors of the PWN remains accurate.This shape can be observed at a later time, specifically 45 kyr after the explosion, as shown in Fig. 1e.As the contact discontinuity surface expands, it encounters Rayleigh-Taylor instabilities due to the significant differences in density and velocity between the pulsar wind and the supernova ejecta.These instabilities are further amplified by the reverberation of the reverse shock from the supernova ejecta, as illustrated in Fig. 1e.

Model with magnetized ISM
During the main-sequence phase of a massive star, the influence of the ISM magnetic field becomes particularly significant.During this phase, the interaction between the stellar wind and the magnetized ISM carves out a large-scale circumstellar wind bubble.This wind bubble plays a crucial role in shaping the propagation of the supernova forward shock.Additionally, the wind bubble's presence influences the pulsar wind's dynamics, further highlighting the interplay between the stellar wind, the ISM magnetic field, and the subsequent evolution of the system.We will describe it in detail in the following.In Fig. 1b (top-right), we can observe the circumstellar medium surrounding the massive star in the presence of a magnetized ISM, as represented in the model Run-35-MHD-0.79-PWN.The black arrows indicate the magnetic field lines of the ISM, which are initially aligned with the polar axis.The overall structure of the circumstellar medium in the presence of the magnetized ISM remains similar to the unmagnetized model (Run-35-HD-0.79-PWN,Fig. 1a).However, the morphology of the shocked shells within the low-density cavity, up to the contact discontinuity between the shocked stellar wind and the shocked ISM, appears to be more elongated along the polar axis due to the influence of the ISM magnetic field.
Indeed, as the expanding stellar bubble interacts with the magnetized ISM, it compresses the magnetic field lines, increasing magnetic pressure and tension along the polar axis.This phenomenon has been extensively studied and described in detail in van Marle et al. (2015b).During the last evolution phase, when the Wolf-Rayet wind material reaches the main-sequence termination shock, it undergoes reflection near the equator.This anisotropic reflection causes a change in the direction of propagation of the shocked material, resulting in the loss of the initially spherical shape of the shocked shell from the Wolf-Rayet wind.The interaction with the magnetized ISM further influences the shape and dynamics of the shocked shell, leading to the observed rectangular morphology of the resulting supernova ejecta.Furthermore, as the expanding supernova blast wave propagates within the elongated cavity (as shown in the left panel of Fig. 1), it interacts with the reflected dense shells resulting from the Wolf-Rayet wind and the elongated contact discontinuity.These interactions lead to anisotropic reverberation at the contact discontinuity of the supernova ejecta.As a result, the shape of the supernova ejecta becomes rectangular, reflecting the influence of the asymmetric interactions with the elongated structure induced by the magnetized circumstellar medium.This mechanism is specifically described within the context of the remnant Puppis A in Meyer et al. (2022a).
In Fig. 1d and f, the influence of the magnetized ISM on the shaping of the PWN can be observed.The ISM magnetic field, which plays a significant role in determining the morphology of the circumstellar medium and supernova blastwave, also affects the confinement and shape of the pulsar wind.Under the influence of the ISM magnetic field, the reflected and the supernova blastwave adopts a rectangular morphology along the direction perpendicular to the magnetic field.This happens because the ram pressure of the supernova ejecta is directed towards the polar axis, causing compression and confinement of the pulsar wind in that direction.In contrast, in the direction parallel to the magnetic field, the pressure exerted by the supernova ejecta is lower, resulting in a more extended shape of the PWN.This interplay between the magnetic field of the ISM, the reverse shock, and the pulsar wind contributes to the complex and asymmetric morphology observed in the PWN, as depicted in Fig. 1d and f.Indeed, the presence of a magnetized ISM influences the propagation of the PWN, resulting in distinct behavior compared to an unmagnetized ISM.In the magnetized ISM model (Run-35-MHD-0.79-PWN), the expansion of the PWN is less pronounced in the equatorial plane compared to the hydrodynamical simulation (Run-35-HD-0.79-PWN),as illustrated in Fig. 1c and d.As time progresses, at a later evolution time of 45 kyr as depicted in Fig. 1f, the pulsar wind continues to be channeled along the direction of the ISM's magnetic field, leading to the formation of a stretched PWN.The presence of the ISM magnetic field affects the dynamics of the PWN and leads to enhanced instabilities at the termination shock of the pulsar wind.These instabilities, which arise from the interaction between the pulsar wind and the magnetized ISM, are more pronounced in the magnetized ISM model (Run-35-MHD-0.79-PWN)compared to the hydrodynamical simulation (Run-35-HD-0.79-PWN).
Our models provide compelling evidence that the morphology of the PWN inside a subsequent supernova, when the progenitor massive static star is located in the Galactic plane, is strongly influenced by the distribution of the magnetic field in the ambient medium.The contrasting evolution and instabilities observed in the magnetized and unmagnetized cases emphasize the significant role played by the interstellar medium's magnetic field in shaping the dynamics and morphology of the PWN.These findings underscore the importance of considering the magnetic field effects when studying the evolution of PWN and their interaction with the surrounding environment.

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

Model limitations
Let us first consider four aspects central to the model.First, the simulations conducted in this study are two-dimensional, assuming axisymmetry and not accounting for variations in the supernova progenitor or the pulsar's spin.While this approach offers computational efficiency and valuable insights, it's essential to acknowledge that a fully three-dimensional treatment is not only important to capture the realistic properties of the ISM but also crucial for a comprehensive understanding of the pulsar wind nebula and the supernova remnant.A 3D model would better represent the complex interactions of the PWN and supernova remnant with the surrounding medium, including the realistic behaviour of magnetic fields.Moreover, the magnetization of the pulsar wind is a fundamental parameter that plays a significant role in the evolution of the PWN and supernova remnant.While we have considered a weak magnetization of the pulsar wind in this study, it's essential to discuss its implications thoroughly.State-of-the-art simulations in both 2D and 3D have shown that the strength and longitudinal variation of magnetization are subjects of debate (Coroniti 2017;Olmi et al. 2016).Future investigations will explore the influence of higher magnetization on the evolution of the PWN in its interaction with supernova remnant and circumstellar medium.
Furthermore, we acknowledge that our modelling of the pulsar wind involves simplified assumptions.A more realistic modelling approach should also involve a better physical description of the wind properties, including its relativistic speed and composition.Addressing these aspects will be crucial in future research for a more comprehensive understanding of the system's dynamics and morphology.Another aspect to consider is the absence of pulsar motion in the simulations.Incorporating pulsar motion would introduce additional complexities and offer a more realistic representation of the interaction between the pulsar wind and the surrounding medium.Furthermore, accounting for the oblique rotation of the pulsar's magnetic axis would allow for a more accurate reproduction of the observed characteristics of the PWN.These are important considerations for future research.The chosen two-dimensional setup and static pulsar position provide valuable insights into general behaviour and trends.However, future investigations can explore the impact of three-dimensional effects, pulsar motion, higher magnetization, and improved modelling of the pulsar wind to obtain a more comprehensive characterization of the system's dynamics and morphology.

Non-thermal emission
To enhance the comparison between our MHD models of the PWN embedded in an elongated circumstellar medium and the available observational data, we performed radiative transfer calculations to gen-erate 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, where here,  represents the gas number density,  is the spectral index and  denotes the energy of the non-thermal electrons in the postshock region of the advancing blast wave.The emission coefficient is given by: being  the observed frequency and  ⊥ the component of the magnetic field perpendicular to the observer's line-of-sight.Intensity maps were obtained by performing the projection given by, where  obs denote the inclination angle of the remnant with respect to the sky plane.These calculations were conducted using the radiative transfer code RADMC-3D4 , and the methodology described in detail by Meyer et al. (2022a).Note that since the investigated numerical simulations are non-relativistic, they do not account for the beaming effect, and this issue will be addressed in our upcoming work.
Figure 2 illustrates the normalized emission maps representing our numerical simulations, specifically Run-35-HD-0.79-PWN(left-hand panels) and Run-35-MHD-0.79-PWN(right-hand panels), showcasing the non-thermal synchrotron emissions in the radio waveband.The top panels correspond to 25kyr, while the bottom panels depict a time of 45 kyr.The intensity is plotted assuming an observer angle ( obs ) of 45 • , representing the angle between the plane of the sky and the plane of symmetry of the supernova remnant.Figure 2a displays the pulsar wind nebula at the age of 25 kyr within an unmagnetized ISM.As highlighted in Meyer & Meliani (2022), no trace of the circumstellar medium is visible in the emission maps because of the absence of the ISM magnetic field.Indeed, the emission map focuses on the pulsar wind and its associated nebula.The image reveals an ovoidal shape, with slightly brighter regions observed at the polar and dimmer regions in the equatorial plane.This brightness variation can be attributed to the toroidal component of the pulsar wind, which applies lateral pressure on the pulsar wind material, causing it to be displaced sideways in the equatorial plane.
At a later evolution time, with a pulsar age of 45 kyr, the radio synchrotron map of the PWN in an unmagnetized ISM is shown in Fig. 2c.The PWN exhibits a jet-torus-like shape, with brighter regions observed at the polar zones.These bright regions result from the strong interaction between the pulsar wind and the supernova ejecta along the pulsar's rotational axis.On the other hand, in the equatorial plane, the strong pulsar wind, driven by the centrifugal force and toroidal magnetic field pressure (Komissarov & Lyubarsky 2004), extends outward.The gas is more diluted in this region, which explains why the equatorial plane is not the brightest region in the hydrodynamical plerion model Run-35-HD-0.79-PWN.In the case of a magnetized ISM, significant changes are observed in the synthetic radio image.The corresponding image is shown at 25 kyr in Fig. 2b.It reveals the presence of two bright arcs parallel to the direction of the ISM magnetic field.These arcs, observed in our axisymmetric setup and aligned with the pulsar's rotation axis, are formed as a result of the interaction between the supernova ejecta and the contact discontinuity between the stellar wind and the magnetized ISM within the elongated cavity (Meyer et al. 2022a).The influence of the ISM magnetic field plays a crucial role in shaping these arcs, ultimately forming a PWN enclosed within a rectangular supernova remnant.
Fig. 2d depicts the older remnant within a magnetized ambient medium, showcasing characteristics of both a supernova shock wave that has interacted with the cavity's border and the growing pulsar wind nebula inside it.The presence of the pulsar wind prevents the reverberation of the supernova shock wave towards the centre of the explosion, as described in Meyer et al. (2022a), resulting in an empty region near the rotating neutron star.The overall morphology of the plerionic remnant still exhibits features of a rectangularly reflected supernova shock wave, with the pulsar wind distributed as an elongated structure.The brightest regions are observed as two polar spots located beyond the termination shock of the pulsar wind.

Comparison with observations
The models presented in this study focus on the evolution of the circumstellar medium surrounding static high-mass stellar objects that eventually undergo supernova explosions, leaving behind a static pulsar.We aim to investigate the formation of elongated pulsar wind nebulae, similar to those observed in Igoshev (2020).It is important to note that these elongated PWN, where the leptonic wind is channelled into the cavity created by the stellar wind shaped by the organized ISM magnetic field, should not be confused with the long tails observed behind the bow shocks of runaway pulsars (e.g., Bucciantini 2002Bucciantini , 2018;;De Luca et al. 2013;Barkov et al. 2019a).
The class of torus/jet-like pulsar wind nebulae, as classified in the catalogue based on Chandra X-ray data, provides strong support for the conclusions drawn from our model.These objects naturally exhibit both an equatorial structure and a jet/counter-jet system, as observed in studies such as Kargaltsev & Pavlov (2010); Kargaltsev et al. (2012) and references therein.Notable examples include the famous Crab nebula with its twisted double jet (Mignone et al. 2013) and the Vela supernova remnant.Magneto-hydrodynamical models have successfully reproduced such structures without considering the stellar wind or supernova ejecta as initial conditions, as demonstrated in Klingler et al. (2014).
The influence of the environment on the morphology of pulsar wind tails/jets has been demonstrated in cases such as the Geminga pulsar wind nebula, which exhibits two curved antennae representing its jets/counter-jet that bend under the influence of the bow shock formed due to the interaction between the fast pulsar motion and the surrounding medium (Posselt et al. 2017).Similar effects have been observed in the case of B0355+54 (Klingler et al. 2014).We propose that the pre-supernova environment plays a similar role, and further modelling efforts are highly desirable, as discussed in Meyer & Meliani (2022).The peculiar morphology of certain pulsar wind nebulae, which cannot be classified as either torus/jet-like objects or bow shock/tail systems, may result from their interaction with a particularly complex surrounding medium.This medium could be shaped by the asymmetric stellar wind during the evolved phases of the progenitor's pre-supernova life, which influences the forward shock of the ejecta and causes aspherical propagation (Velázquez et al. 2023;Villagran et al. 2023).

CONCLUSION
This paper presents a study on the modelling of PWN in core-collapse supernova remnants associated with static massive stars in the warm phase of a magnetized spiral arm of the Milky Way.By utilizing 2.5-dimensional simulations, we demonstrate that the reflection of the supernova blast wave against the elongated contact discontinuity between the stellar wind and magnetise ISM of the magnetically elongated stellar wind cavity in the progenitor's circumstellar medium has a significant impact on the morphology of the resulting PWN.This phenomenon might be responsible for forming rectangular supernova remnants, such as Puppis A, as described in Meyer et al. (2022a).The reverberation of the shock wave leads to the compression of the pulsar wind and imposes a preferred expansion direction perpendicular to the plane of the pulsar's spin.As a result, the PWN within the rectangular supernova remnant becomes elongated rather than adopt-ing the jet-torus-like shape typically observed in previous studies, as described by Komissarov & Lyubarsky (2004).
The radio synchrotron emission maps of plerionic supernova remnants exhibit a complex morphology that evolves over time.Initially, the morphology is characterized by a young, growing, ovoidal PWN combined with the rectangular shape produced by the interaction between the supernova ejecta and the walls of the unshocked stellar wind cavity of the progenitor star.This interaction gives rise to the rectangular appearance observed in Puppis A, as discussed in Meyer et al. (2022a).As time progresses, the influence of the ISM magnetic field becomes more prominent in shaping the remnant's morphology.The channelling effect of the pulsar wind into the elongated circumstellar wind cavity of the progenitor extends along the pulsar's rotation axis.Instabilities at the interface between the pulsar wind and the ejecta result in a knotty nebula, manifesting as bright spots within the plerion.The irregular shapes observed in many pulsar wind nebulae may indicate the complex nature of the surrounding environment, influenced by both the distribution of material in the ambient medium and the stellar wind history of the supernova progenitor.In this complex environment, the interaction between the supernova ejecta and the pulsar wind gives rise to observed irregular morphologies.

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

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

Table 1 .
List of models in this study.All simulations assume a rotating static massive star of mass  ★ at solar metallicity, in a medium of number density  ISM and organised magnetic field strength  ISM .The initial rotation rate of the central massive star is Ω ★ /Ω K = 0.1.