Core-collapse supernova inside the core of a young massive star cluster: 3D MHD simulations

Young massive stars in compact stellar clusters could end their evolution as core-collapse supernovae a few million years after the cluster was built. The blast wave of a supernova propagates through the inner cluster region with multiple stellar winds of young luminous stars. We present the results of 3D magnetohydrodynamic simulations of the plasma flows produced by a supernova event inside a cluster with a population of massive stars similar to that in Westerlund 1. We followed its evolution over a few thousand years (i.e. a few shock crossing times). The plasma temperature, density and magnetic field, which are highly disturbed by supernova event, relax to values close to the initial over the studied period. The relaxation time of a cluster is a few thousand years, which is a sizeable fraction of the period between the successive supernova events for a massive cluster of a few million years age. The spectra of the cluster diffuse X-ray emission simulated here should be representative for the galactic and extragalactic young massive clusters. The resultant magnetic fields are highly intermittent, so we derived the volume filling factors for a set of magnetic field ranges. Highly amplified magnetic fields of magnitude well above 100 $\mu$G fill in a few per cent of the cluster volume, but still dominate the magnetic energy. The structure of the magnetic fields and high velocity plasma flows with shocks in the system are favorable for both proton and electron acceleration to energies well above TeV.


INTRODUCTION
Young massive star clusters (YMSCs) host rich and dense populations of luminous massive stars which have short lifespans of a few Myr.Some of these stars end their lives as core-collapse supernovae (SNe) releasing ∼ 10 51 erg of kinetic energy into the surrounding medium.The interaction of powerful winds from luminous stars can create a complex and highly magnetized medium with numerous shocks.When the forward shock of a supernova remnant (SNR) propagates through the circumstellar medium of massive stars it is expected to further compress the gas and amplify the magnetic field.In this paper, we aimed to investigate the dynamics of a supernova remnant in the context of a YMSC and its impact on the inner cluster environment.
YMSCs play an important role in star-forming regions being the powerful sources of ionising radiation, kinetic energy and momentum, which affect their parent clouds (see for review Krumholz et al. 2019;Krause et al. 2020) and the interstellar medium (ISM) (e.g.Müller et al. 2020).The burst mode of episodic accretion can explain the high-mass star formation processes (see e.g.Meyer et al. 2019Meyer et al. , 2021)).In the galaxies with high star formation rate a substantial fraction of massive stars are born in YMSCs (Adamo et al. 2020).Moreover, they were suggested recently as likely places for the formation of binary black holes which are the sources of gravitational waves (e.g.Rastello et al. 2023).Both stellar and diffuse emission ★ E-mail: danirbadmaev.astro@gmail.comfrom some clusters was thoroughly studied over the whole range of the electromagnetic waves including the high energy photons.
Diffuse X-ray emission from hot plasma produced by colliding stellar winds (SWs) and SNe was revealed in YMSCs (see Muno et al. 2006;Townsley et al. 2019;Kavanagh 2020;Sasaki et al. 2022;Bykov et al. 2023).Some clusters were also detected in very high energy gamma-rays (e.g.Aharonian et al. 2007Aharonian et al. , 2022;;Ohm et al. 2010;Yang et al. 2018) which is a clear indication of efficient cosmic ray (CR) acceleration processes in these sources which can be studied in detail with the forthcoming Cherenkov Telescope Array observatory (see e.g.Acero et al. 2023;Acharyya et al. 2023).CRs accelerated in the clusters can penetrate deep into the cores of nearby molecular clouds providing the gas heating and ionization.This effect is known to be important to both star formation process and interstellar chemistry (e.g.Draine 2011;Padovani et al. 2020).Indeed, due to the high extinction the UV and soft X-ray photon fluxes drop drastically inside the dense cloud clumps and the low-energy CRs are the main ionization agent.Being localized inside the molecular clouds and rather long-lived (up to some Myr) the low-energy CR sources associated with YMSCs can be an important feedback agent in the star formation process in galaxies.On the other hand, the high-energy protons accelerated in YMSCs interacting with the surrounding matter can contribute to the observed high-energy neutrino fluxes (see e.g.Bykov et al. 2015).An adequate interpretation of the high-energy diffuse emission requires detailed modeling of the energy partition processes in the complex plasma flows, produced by the colliding winds of young massive stars and SN events.Models of particle acceleration in YMSCs at the evolution phase, dominated by the powerful winds of OB and Wolf Rayet (WR) stars and SNe (it is discussed in Bykov 2014; Gupta et al. 2020;Morlino et al. 2021;Vieu & Reville 2023;Gabici 2023), rely on the dynamics of the shocks of different strengths and scales, as well as on the magnetic fields in clusters and their vicinity.The hydrodynamic modeling of the structure of YMSCs and their impact on the galactic environment was performed extensively (see e.g.Chevalier & Clegg 1985;Stevens & Hartwell 2003;Torres & Domingo-Santamaría 2007;Rogers & Pittard 2013;Wünsch et al. 2017;Gallegos-Garcia et al. 2020).To study the structure and evolution of magnetic fields in YMSCs at the stage dominated by powerful stellar winds we performed a 3D magnetohydrodynamic (MHD) simulation of the quasi-stationary medium inside a YMSC core (see Paper I: Badmaev et al. 2022).In Paper I we discussed the formation of cluster-scaled magnetic field filaments, where the amplitude of the field reaches > ∼ 100 G .This model did not include possible supernova events and thus it is only applicable for a period of time (< 10 kyr) when there are no strong and rapid ( dyn ≲ 1 kyr) perturbations as the stellar system evolves.Chevalier (1982) and Nadezhin (1985) proposed a self-similar solution for an SNR freely expanding at the earliest stages of evolution.This solution considers the expansion of the supernova ejecta into the surrounding medium and can be applied to different types (Ia, Ib/c-II) of SNe occurring in either a uniform or wind-blown circumstrellar medium.Hydrodynamic modeling of supernova remnant evolution inside a centrally symmetric wind-blown bubble was performed, considering both stationary (Blinnikov et al. 1982;Chevalier & Liang 1989;Tenorio-Tagle et al. 1990) and moving massive progenitor (Rozyczka et al. 1993), demonstrating complex behaviour of the ejected material and velocity fields.The origin of strong magnetic fields in young SNRs has been investigated through 2D MHD simulations by Jun & Norman (1996).Models with a smooth transition between the free-expansion and Sedov-Taylor phases have also been introduced (see Truelove & McKee 1999).The evolution of a SN shock in the bubble of the massive progenitor star which accounted for different mass-loss regimes at the evolution phases of O star, red supergiant, and WR was simulated by Dwarkadas (2007).MHD models have also been applied by Marcowith et al. (2018) to the problem of CR acceleration.The effect of magnetic fields in the wind of the progenitor star of SN1987a has been studied in detail by Orlando et al. (2019).Recent developments include MHD models for the dynamics of core-collapsed SNe (Petruk et al. 2021), connecting the early phases of supernovae to the observed morphology of SNRs in multi-wavelength observations (see e,g.Vink 2020;Ferrand et al. 2021).Some of the latest 3D HD/MHD models demonstrate the connection between the dynamical properties of SNRs and the internal structure of their progenitor stars, as well as the circumstellar medium (see Orlando et al. 2020Orlando et al. , 2022)).Non-spherical morphology of SNRs propagating in anisotropic wind-blown bubbles shaped by interstellar magnetic fields has also been investigated by Meyer et al. (2022).Recently, Meyer et al. (2023) have studied the mixing of materials in magnetized core-collapse SNRs moving through the ISM over a period of ∼ 10 kyr.Additionally, collisions between supernova remnants and stellar winds from nearby massive stars have been simulated using both hydrodynamic and MHD models (e.g.Badmaev & Bykov 2021;Velázquez et al. 2003Velázquez et al. , 2023)).
Given the number of massive stars ( > 25  ⊙ ) in a YMSC core ∼ 100 and considering YMSCs of ages below 5 Myr, one could expect the average rate of supernovae (SNe) to be as high as ∼ 0.1 kyr −1 as it was estimated by Muno et al. (2006) for a rich galactic stellar cluster Westerlund 1 at its current evolution stage.In reality this implies that the system withstands recurrent perturbation and relaxation phases caused by the propagating SNe shock fronts.We performed 3D MHD modeling of SNR evolution inside YMSC in order to study the energy partitioning, thermal and magnetic structure of a cluster right after the SN event.The metal rich matter ejected by SN is expected to change the plasma composition in a cluster for ∼ 1000 years and thus may affect the X-ray spectrum of the hot cluster interior.Shocks and magnetic field structure are the key ingredients of non-thermal particle acceleration and radiation models in YMSCs.
The paper is organized as follows.The 3D MHD model using the PLUTO code which includes the numerical scheme, stellar cluster setup and the supernova initialization are discussed in §2.The simulated plasma density, temperature and magnetic fields distributions in YMSC perturbed by supernova events are presented in §3 where 3D rendering images and plane maps both for the central and peripheral locations of SN are considered.The temporal evolution of the SN, energy partitions, statistical distributions of the cluster magnetic fields and SN ejected mass are presented in §4.In this section we also illustrate the X-ray spectra of the hot optically thin plasma and discuss in brief the non-thermal components in the simulated cluster.

Governing equations
The simulations were performed using the well-proven open source code pluto (Mignone et al. 2007(Mignone et al. , 2012;;Vaidya et al. 2018) based on the Godunov method, and created specifically for problems of computational astrophysics.According to our problem, the code integrates the following set of non-ideal magnetohydrodynamic equations: where  =  represents the momentum density vector of a control volume,  is the magnetic field vector,  is the identity matrix,  tot =  +  • /2 is the total pressure.The total energy density of the systems reads, and the sound speed,  s = √︁  /, closes the above system of equations, where  = 5/3 is the adiabatic index.The source term Φ (, ) on the right-hand side of the total energy conservation equation represents optically thin radiative losses and heating.The plasma heat flux is determined by the vector  c .
We took into account the gains and losses by optically thin radiative processes following the recipe from Meyer et al. (2014); Green et al. (2019) for the case of photo-ionization equilibrium: where Γ () and Λ () are the radiative heating and cooling rates, respectively, and  H is the hydrogen number density.The efficiency of thermal conduction in the collisionless turbulent plasma is still under debate.The recent work of Meinecke et al. (2022) revealed strong suppression of the thermal conduction in the turbulent, weakly collisional, magnetized plasma by more than two orders of magnitude comparing with the Spitzer's thermal conduction in collisional plasma.Strong supernova shocks are expected to produce a highly enhanced level of plasma turbulence (at scales well below of that resolved in the MHD simulations) which would suppress the thermal conduction on the time scales studied in the paper.The effect of the Bale et al. (2013) solar wind thermal conduction at the appropriate time and length scales on the plasma flows in SSC was studied in Paper I.
This system of equations is solved by use of the unsplit secondorder Runge-Kutta algorithm with linear reconstruction of the variables between adjacent cells together with the combination of HLLD (Miyoshi & Kusano 2005;Mignone et al. 2007) and HLL (Harten et al. 1983) Riemann solvers.The divergence-free condition for the magnetic field is ensured by the Hyperbolic Divergence Cleaning algorithm (Dedner et al. 2002) in the whole computational domain.The time-marching algorithm is controlled by the standard Courant-Friedrichs-Lewy parameter that we set to  cfl = 0.2.

Supernova remnant initialization
For the simulations of the SN event inside the star cluster environment we used a mapping strategy (e.g.Meyer et al. 2015).Assuming the SNR shock to be isotropic, which is sufficient for our scope, we calculated its early expansion in 1D with high resolution.This simulation accurately tracks the evolution of an SNR starting from the time  0 ≈ 10 −2 yr after the core-collapse up to the time  map ≈ 30 yr as it propagates through the wind of a WR-star progenitor (i.e.SN Type Ib/c).Then, we mapped the obtained 1D solution into the 3D domain, which contains the pre-simulated ISM of an SSC from Paper I with randomly oriented stellar spins (see it for the details on the interacting winds model), replacing the fixed wind injection region of the progenitor star.
For the SNR model we followed a standard procedure (e.g.Chevalier & Liang 1989;Truelove & McKee 1999;Whalen et al. 2008;Telezhinsky et al. 2013;Petruk et al. 2021) that leads to the classical self-similar type solution of Chevalier (1982) and Nadezhin (1985), see Fig. 1.It was assumed that at the time  0 =  0 / 0 shortly after the SN event the ejecta expands freely at the velocity  0 = 30000 km s −1 , and consists of an inner (core) and outer layers with density profiles are being, respectively, uniform and steep: where the index  = 9 corresponds to a Type Ib/c remnant (see Dwarkadas 2005), and  = /4  is a progenitor's wind parameter.Here,  and  c =  c / 0 are normalization constants that must be determined from the required ejecta mass  ej = 10 M ⊙ and SN event energy  ej = 10 51 erg.In general, in order to find these constants one has to follow the two-step numerical procedure described in Whalen et al. (2008), or, if the inequality  c ≪  0 deliberately holds, use the following analytic expressions: In our case the circumstellar environment around the expanding SNR is a free WR wind that inherits its parameters from Paper I, hence:  = 6.50 × 10 −5 M ⊙ yr −1 , and  w = 1600 km s −1 .We did not consider the magnetic field in this 1D SNR solution, as (1) it is not dynamically important at the earliest free-expansion SNR phase (see Das et al. 2022), (2) it is assumed to contribute a quantitatively negligible fraction to the interstellar magnetic fields that will be compressed by the SNR forward shock.
We used a simplified description of the wind-driven circumstellar medium around the supernova progenitor star which is specific for a SN in a YMSC.More complex model could account for evolutionary changes in the stellar mass loss rate (see e.g.Dwarkadas 2022) which is important for isolated massive stars.In the case of SN in a YMSC of a few Myr age, which we are discussing here, the surroundings of the progenitor star are rather determined by the colliding winds of massive stars on sub-parsec scale sizes.Also, the possible previous SNe could sweep out the matter from the cluster and thus erase the stellar evolution history to a significant extent.
Once the initial 1D SNR solution was mapped into 3D MHD domain the pluto code solved the equations (1-4) in Cartesian grid (, , ).The computational domain was extended in the intervals of [−2.16; 2.16] pc in all directions and covered with a uniform grid of 540 3 and 270 3 cells for the first 1500 yr and the last 6000 yr (downscaled relaxation part, see § 4.1) of integration time, respectively.At the domain borders we used a modified 'free outflow' boundary condition that prohibits any possible backflow of the gas.We performed two simulations for the two cases of the SNR progenitor position inside the cluster core: near the cluster center,  dist ≈ 0.6 pc ('central'), and at the cluster periphery,  dist ≈ 2 pc ('peripheral').In the latter case it was expected that the SNR shock propagating through the whole cluster volume would sweep and compress more material flowing both towards (in the first half of the cluster volume) and away (in the second half) from the SNR forward shock, possibly resulting in the more effective Axford-Cranfill-type magnetic field amplification.
The ejecta material was traced using the scalar marker  which was described by the linear advection equation: This marker was passively advected with the fluid, allowing us to distinguish between the ejecta and interstellar material.We set the value () to be negative for the ejecta and positive for the SWs and interstellar material, here  is the vector position of a given cell in the simulation domain.

RESULTS
Here we present the results of two simulations from  snr = 100 yr up to  snr = 500 and  snr = 700 yr for the centrally located and peripheral SN events inside the YMSC core, respectively.These moments in time track the expansion of the SNR until the moment when its forward shock reaches the computational domain borders.
We use two forms of presentation: 2D maps (cross-sections of a 3D domain) and 3D volumetric renders of scalar and vector fields.
For both simulations the integration times up to  snr = 7500 yr are thoroughly analysed in §4.

Overall dynamics and structure of flows
We have obtained a 3D MHD data on the passage of the SNR shock wave through the Wd1-like YMSC core at a resolution of 8 × 10 −3 pc cell −1 over the dynamical time  dyn = / prop ∼ 1200 yr, where  prop ∼ 3000 km s −1 is the SNR shock propagation speed and  ≈ 4 pc is the length scale.The results are presented in 2D maps of density, temperature, and magnetic field amplitude for two cases of SN location: in the center of the cluster core (see Fig. 2) and on its periphery (see Fig. 3).The figures show cross-sections perpendicular to the -axis passing through the points where the SN progenitors were located.Flow streamlines were plotted over the density maps using velocity vector field data.Since the initialized remnant is spherically symmetric, the morphology of the expanded SNR shell depends on the non-uniform and intermittent collective flows inside the simulated cluster.Therefore, it is important to note that different cross-sections may reveal slightly different structures of the SNR shell in detail.
The local structure and morphology of the expanding SNR shell are determined by the winds of neighboring stars.The forward shock crushes through the SWs within the cluster core, creating numerous bow shocks with various geometries depending on the wind type.The size and geometry of a bow shock depend on the balance between the ram pressure of the SNR and the kinetic power of the wind (see Wilkin 1996).The widest and thickest bow shocks are formed around cold and dense CSG-winds, as shown in Fig. 2. In the areas of impact, SWs form cometary structures that penetrate through the supernova remnant shell and leave heated up 'dents' on its surface, see Fig. 4. Magnetic fields are amplified to approximately 100 G at the bow shocks and tails of these cometary structures.The overall spherical structure of the remnant is preserved during expansion.When the SNR shell leaves the volume of the domain, the SWs relax to their roughly unperturbed sizes in ∼ 3000 yr, as shown in Fig. 5.

The SNR expansion
Just a hundred years after the SN event one could clearly see the remnant with nearly undistorted initial structure, see Figs. 2-3.At this stage the SNR expands freely and it has a typical structure with two shock waves (reverse and forward) and a contact discontinuity (CD) (see e.g.Truelove & McKee 1999).
The latter can be easily identified by the Rayleigh-Taylor (RT) instability modes at the surface of the inner dense ( ∼ 30 cm −3 ) and thin (∼ 0.1 pc) shell shown at the top of Fig. 4. In the immediate vicinity of the CD, there is a reverse shock that heats a thick but dense layer of SN ejecta.The lack of a reverse shock moving towards the center, as explained in (see Petruk et al. 2021), can be attributed to the high amount of ejected mass compared to the mass of the cluster diffuse gas and the short time scales.As a result, no extensive interior heating is observed, and most of the SNR volume cools adiabatically.However, this shock effectively heats the inner dense SNR shell to the temperatures ∼ 10 7 K.On the other hand, the forward shock is responsible for heating the swept-up interstellar medium plasma to temperatures up to a few 10 8 K, see the temperature data in Fig. 4. It is also responsible for effective amplification of interstellar magnetic fields well above 200 G in some filaments and thin isolated regions, see the bottom of Fig. 4.This amplified field is probably carried into the bow shock tails, this is clearly seen in Fig. 3.The thickness of the shocked interstellar gas layer is ∼ 0.5 pc with the density  ∼ 4 cm −3 .Some minor differences in thickness of the shocked gas layer and dense RT-unstable shell are the consequence of different initial positions of the SNRs.No significant quantitative differences in terms of densities, temperatures, and magnetic field strength are detected for the case of peripheral SN over the central one.

Intracluster medium relaxation
For a stellar cluster similar to Wd1 the most massive stars are likely to produce the first SN events when the cluster is ∼ 3 Myr old, and then, following the estimation by Muno et al. (2006), one should expect ∼ 1 SN every 7-13 kyr for the time interval of duration ≳ 1 Myr.It is important therefore to derive the relaxation times of the intracluster plasma after the SN event.These times are somewhat different for the dynamics, energy partition, and ejected material admixture.The cluster crossing time by the SN shock is ≲ 1000 yr depending on the progenitor star position.The simulation presented above shows that the violent disturbance of the intracluster medium produced by the core-collapse SN shock flow relaxes to the quasi-stationary state 4000-5000 yr after the SN event.This is consistent with what we see in Fig. 5, where the general structure of flows is recovered at ∼ 3500 yr and during the next several thousand years the small-scale density perturbations are being suppressed.In terms of kinetic and thermal energy partitioning, as it is apparent in Fig. 6, it takes ∼ 4000 yr for the system to reach the relaxed state.Yet one should note that in the case of the peripheral SN event the sub-dominant magnetic energy part, while slowly reaching some plateau at 7000 yr, still struggles to recover back to its base level at ≈ 3 × 10 47 erg.This could be both the effect of resolution down-scaling after  snr = 1500 yr and the consequence of slow magnetic field reorganisation.If we track the ejecta admixture (see Fig. 7), the relaxation takes place at ∼ 4000 yr, when highly diluted leftovers of the ejecta finally leave the domain volume, regardless of the SN event initial position.
The energy partition and total mass of the gas presented in Figs. 6 and 9 as functions of time were integrated explicitly: where     () represents the scalar density field,  cell is the domain control volume, and  cell = 270.For the mass loss rate estimate we calculated the integral of the flow     ()    () through the faces of the cubic computational domain.

Magnetic field evolution
Magnetic fields in the simulated core of a YMSC are highly intermittent and span a broad range of magnitudes from a few to hundreds of G with various volume filling factors (cf.Badmaev et al. 2022).The energetic SN event severely disturbs the quasi-stationary magnetic field configuration as well as the other important characteristics inside the cluster (see §4.1).The 3D render of the magnetic field structure for  = 500 yr after the centrally located SN event is shown in Fig. 4.
The magnetic energy of the simulated cluster core is dominated by the magnetic fields of amplitudes ≳ 100 G with a few per cent volumetric filling factor.This is consistent with the results of Inoue et al. (2009) 2D simulations of a young supernova remnant in a turbulent environment.One can see that the magnetic field magnitude reaches the values ≲ 750 G in some highly compressed regions of the filamentary structures, see Fig. 4. The evolution of the magnetic structures over time can be observed in Fig. 8.It shows the changes in the filling factors of magnetic fields of different magnitudes for both the central and peripheral supernova events.The temporal evolution differs between central and peripheral supernovae, but in both cases, a quasi-stationary distribution is reached after ∼ 4000 yr.
The decrease in volumetric filling factors observed in high amplitude magnetic field bands can be attributed to the crushing and sweeping of the initial magnetic field configuration which was previously formed by the colliding winds of massive stars.The strong forward shock driven by the SNR travels through highly inhomogeneous cluster medium.It compresses the interstellar magnetic field immediately behind the collisionless shock and possibly amplifies the fluctuating fields downstream through dynamo-like effects.The thin shell of high magnetic field accompanies the SNR forward shock and the downstream flows, overrunning powerful stellar winds as depicted by magnetic field data in Figs.2-3 and 4.
The total magnetic energy of the cluster is ∼ 10 47 erg and exhibits some apparent temporal evolution especially during the central SN event (see Fig. 6).The magnetic energy has a prominent minimum that occurs when the SNR shock front exits the domain volume and lasts until the fast colliding SWs start restoring the magnetic fields in the cluster core, before the next possible SN event.It is evident that during the propagation of the SNR shock through the cluster core volume, the range of magnetic field magnitudes || between 3 and 30 G becomes spatially dominant.The total energy in magnetic fields is mainly contributed by fields with magnitudes well above 30 G, with a significant contribution from regions with magnitudes higher than 100 G, as seen in Figs.6-8.The presence of these high magnetic fields can be tested by searching for synchrotron X-rays from multi-TeV leptons, which are expected to be present in TeV sources like Westerlund 1.
In Fig. 8 the volumetric filling factors for the magnetic fields in range  = ( min ;  max ) were calculated as follows: where  stands for the total domain volume, and |    ()| is the magnetic field magnitude.

Ejecta material transport
In this paper we modelled the ejecta relaxation in the compact cluster with multiple powerful SWs at a few kyr time-scale.This process is relevant to the YMSCs evolution stage dominated by the most massive stars.While the radiative losses of the hot gas in the cluster do not affect seriously the plasma dynamics in any case, the account of the metal enhancement from the SN ejecta material could be important for the X-ray diffuse emission spectra.
During the first ∼ 1000 yr when the shock wave from a supernova propagates through the cluster, the ejecta replaces interstellar matter in about 50-80 per cent of the cluster core volume, depending on the position of the SN event, see Figs. 7 and 10.A significant amount of this ejected material remains within the cluster core for ∼ 3000 years after the SN event, when 10 M ⊙ of material were expelled.It is interesting to note that ≈ 2000 yr after the SN event there is a total of 17-18 M ⊙ of diffuse gas within the cluster core, of which only ∼ 0.7-1 M ⊙ is the ejected material.Even though it is heavily diluted, this material still occupies about 25-30 per cent of the volume, see Fig. 7.The rate at which the cluster loses mass remains relatively stable at around ∼ 2.5×10 −3 M ⊙ yr −1 during the considered stage of evolution.However, when the SNR shock front reaches the domain boundaries and the dense shell created by the SN starts to leave the volume, the mass loss rate increases by an order of magnitude.As a result, the total mass of gas within the cluster decreases rapidly, but then begins to grow linearly due to the effective CSG-wind mass loss, see Fig. 9.After  snr ∼ 4000 yr there is almost no ejecta left inside the cluster core.Additionally, we conclude that the SN shock wave has a weak impact on the structure of neighboring heavy CSG-wind shells within a timescale of 10 kyr, compare Figs. 3 and 5.
In Fig. 7 we calculated the ejecta mass fraction using passive scalar tracer , which was set to be negative for the ejecta material (see §2.2): The ejecta volumetric filling factor was treated following Eq. 12 with the above condition for (, , ).

X-ray diffuse emission
Some of the YMSCs are bright X-ray sources (e.g.Westerlund 1, Westerlund 2).MHD simulations presented in this paper and Paper I allow modelling the thermal diffuse emission of the cluster core.
Together with the analysis of the X-ray spectra of YMSCs it allows revealing the nature of their X-ray radiation, namely distinguishing between thermal and non-thermal origin of the diffuse emission.In Paper I we presented the modelling of thermal X-ray spectra of the cluster core, containing tens of massive stars, but with no SNe.Here we analysed the impact of the latest SN event in a cluster on X-ray emission spectrum.
During the first 1000-2000 years after the SN event the cluster contains several solar masses of the ejecta.It should be taken into account that the ejecta of Type Ib/c SNe has different from the ISM chemical composition, rich in metals.As the SN ejecta was marked in our simulation, we were able to calculate the thermal spectra, considering the different chemical composition in different parts of the cluster.While for the stellar winds we used the standard table of abundances given in Asplund et al. (2009), which is implemented in     APEC model of the XSPEC, for the ejecta material we constructed the new table of abundances of elements up to Zn using the functionality of VVAPEC.We got the values of elemental abundances in SN Type Ib/c ejecta from Limongi & Chieffi (2018), where the detailed modelling of stellar nucleosynthesis is presented and SN ejecta yields are found.
Even when SN shock has left the cluster core, the ejected material can still substantially enrich with the metals the hot X-ray emitting plasma in the core as it is apparent from Fig. 7 and §4.3.Most of the ejecta filling the cluster volume is cold (see lower left panels of Fig. 2, 3), but in the areas where ejecta collides with fast O-WR winds, it is heated up to a few keV, which is shown in Fig. 10.These areas can affect significantly the normalization and features of the X-ray diffuse thermal spectrum of the cluster.Integrating the emission (found with APEC and VVAPEC models) from each cell of the computational domain with its own ,  and chemical composition ('standard' or 'ejecta'), we managed to obtain the thermal spectrum from YMSC as a whole.
The MHD model we used here is single-fluid and therefore we must set an appropriate electron temperature prescription for the spectra calculation.We assumed that the hot plasma in the fast SN shock downstream has an initial  e ≃ 0.1 p as it was argued by Raymond et al. (2023); Vink et al. (2015).The ions heated there to ∼ 10 keV have not enough time to reach the thermal equilibrium with electrons in the first ∼ 1000 years after the SN event.Thus the temperature of electrons in the hot plasma in the SN shock downstream  e ≃ 0.1 (where  is plasma temperature in one-fluid model) is about keV.On the other hand, the older material, heated by the stellar wind flows before the SN event, had longer time to relax to  e ≃  p and the electron temperature there is  e ≃ /2 which is in keV range as well.
We present the modelled diffuse spectra in the energy range 0.5-12 keV for times up to 900 yr after the SN event in Fig. 11 where ejecta's  e ≃ 0.1 p .The diffuse spectrum before the SN event (same as in Paper I) is also presented.Despite the spectra being given in arbitrary units, the relative fluxes are saved.One can see that after hundreds of years from the SN event the intensity of the X-ray emission is increased due to compression of the keV temperature plasma by SN shock.The metal-rich composition of the ejecta also leads to the increase of the X-ray flux.The decrease of the intensity with time after 300 years of the simulation is because we present the flux from a fixed domain and the compressed SN shell just moves out from the domain.The decrease of all energy profiles in Fig. 6 after ∼ 1000 years is due to the same reason.The absence of thermal equilibrium between electrons and ions, as well as the fact that in ∼ 1000 years most of the heated matter is swept out of the cluster, leads to the softening of the spectrum, and at high energies it can even be lower than the spectrum, calculated before the SN event.
The thermal spectra, modelled in this paper, correspond to time periods, close to the SN event and can be used for the analysis of X-ray observations of YMSCs.Such analysis for Westerlund 2 was performed in Bykov et al. (2023) and it was found that in addition to thermal component, which temperature was obtained with MHD modelling, the presence of a non-thermal component is necessary to fit the Chandra and ART-XC observations.Similar discussion on the origin of the diffuse emission from other YMSCs, e.g.Westerlund 1, is needed and may be the subject of the future work.

Nonthermal particle acceleration and radiation
The knowledge of magnetic field structure is particularly important for modeling of very high energy particle acceleration, their confinement and the synchrotron radiation of the accelerated leptons.About 30% of the cluster volume is filled with the low magnetic fields < 3  and more than a half of the volume has the field magnitude below 10 .The intracluster plasma is perturbed by multiple shocks and the energy-containing MHD motions of a parsec scale, which are supported by the stellar winds and supernovae.Particle acceleration timescale in kinetic models of CR acceleration in compact clusters and superbubbles (which differ by scale sizes) can be estimated as  a ( ) ≈ 9 ( )/ 2  (see e.g.Bykov & Kalyashova 2022), where  ( ) is the diffusion coefficient of a CR particle of momentum , and the velocity dispersion of the bulk compressible plasma motions in the system   is ≳ 1000 km s −1 .The maximal energies of CR accelerated by the Fermi mechanism in the system are limited by the CR loss time   ( ) which is the shortest among the CR escape time or the synchrotron-Compton losses time (for leptons).One may get an estimation for the maximal momentum  max of accelerated CRs from the equation  a (  max ) =   (  max ).Using the approximation of the electron energy loss rate with account taken of the Klein-Nishina suppression of the Compton losses on the intense optical radiation field following Moderski et al. (2005) one can get   max ≳ 10 TeV for a cluster like Westerlund 1 or Westerlund 2. The accelerated leptons of these energies can radiate keV regime X-ray synchrotron radiation in the filaments of high magnetic field magnitudes ≳ 100  which fill about 5% of the cluster volume and dominate the magnetic energy of the system.In Bykov et al. (2023) we performed such calculation for Westerlund 2, which allowed us to explain the presence of a nonthermal component in its diffuse X-ray emission.Also, when SN shocks are passing through the close vicinity of a hot luminous star, the inverse Compton radiation of MeV electrons accelerated at the shock can contribute to the X-ray emission of the cluster.Thus, the clusters are expected to be the sources of non-thermal X-ray radiation.The hot keV plasma in the cluster also emits in X-ray energy range, which was observed in a number of clusters (see e.g.Muno et al. 2006;Townsley et al. 2019;Kavanagh 2020;Bykov et al. 2023) and is discussed in §4.4.

CONCLUSIONS
We present the 3D MHD simulations of the structure and evolution of the plasma flows, temperature and magnetic fields in the core of a young massive star cluster which is disturbed by an SN event.Both central and peripheral positions of a core-collapsed SN of energy 10 51 ergs and 10 M ⊙ of mass ejected in a wind driven circumstellar medium are considered.The model concerns young massive clusters of a few million years age.
The SNR forward shock crushes through the YMSC core, changes its elemental abundances, heats the gas, and creates multiple bow shocks.SWs form magnetized cometary structures which penetrate into the SNR shell leaving heated up dents on its surface.This highly violent plasma state is relaxed to the quasi-stationary state over a period of ∼ 4000-5000 yr, while for the magnetic field structure it takes substantially longer.No significant differences in terms of densities, temperatures, and magnetic field strength are detected for the case of peripheral SN over the central one.
The structure of the cluster's perturbed magnetic fields is highly intermittent and shows a bunch of evolving elongated magnetic filaments of high magnitudes ≳ 100 G (up to ∼ 750 G) and the extended regions with magnetic fields of a few 10 G.During the propagation of the SNR shock, the range of magnetic field magnitudes || ∈ (3; 30) G becomes spatially dominant, yet the magnetic energy is mainly contributed by the fields of amplitudes ≳ 100 G with a few per cent volumetric filling factor.The presence of highly amplified magnetic fields allows the production of non-thermal X-ray synchrotron radiation by multi-TeV particles accelerated by shocks within the cluster.
The metal-rich material ejected by core-collapse supernova reaches the maximal volumetric factor of 50-80% (depending on the SN position) inside the cluster in about a 1000 yr after the SN event, which then decreases to < 5% in the next 3000 yr.Some relatively minor fraction of the SN ejecta material is heated to temperatures above a keV within the cluster core region affecting the thermal X-ray spectrum.Right after the SN event the thermal radiation (1) is increased due to the growth in thermal energy released in a cluster (2) has spectral features, connected with the mixing of a cluster material with the metal-rich SN ejecta.

Figure 1 .
Figure 1.The obtained 1D SNR solution profiles of density, pressure, and velocity at the moment of mapping  map .Here the forward shock is located at  = 0.38 pc and the contact discontinuity between the SN ejecta and progenitor's wind material is right at the density profile spike.

Figure 2 .
Figure 2. The evolution of the SNR placed near the cluster center.Top: density maps and flow streamlines (dark translucent lines) at  snr = 100 and  snr = 500 yr, left to right.Bottom: temperature and magnetic field magnitude maps at  snr = 500 yr.The white lines 1, 2 and 3 mark the positions of the reversed shock, contact discontinuity and forward shock, respectively.In the top left corner the dense CSG-wind shell has formed a wide bow shock with a magnetized tail.The maps are captured in  -plane at  = −2.7 in units of [0.1 pc].

Figure 3 .
Figure 3.The evolution of the SNR placed at the cluster periphery.Top: density maps and flow streamlines (dark translucent lines) at  snr = 100 and  snr = 700 yr, left to right.Bottom: temperature and magnetic field magnitude maps at  snr = 700 yr.The white lines 1, 2 and 3 mark the positions of the reversed shock, RT-unstable contact discontinuity and forward shock, respectively.At the bottom right corner there is a cold and dense CSG-wind shell.The maps are captured in  -plane at  = −10.7 in units of [0.1 pc].

Figure 4 .
Figure 4.The density (top), temperature (middle) and magnetic field (bottom) distributions at the surface of the cluster-centered SNR shell shown as 3D volumetric renders,  snr = 500 yr.The red bulb in the top figure is a CSG-wind shell.The orange 'dents' in the middle figure are the SWs impacting the SNR.

Figure 5 .
Figure 5. Relaxation dynamics after the peripheral SN event.At the time of 3500 yr the overall flow geometry is recovered, and after the next 2500 yr the small-scale density perturbations are smoothed.At the lower right corner of each frame one can see a dynamically stable dense shell which is formed by a CSG-wind.The cross-section here is the same as the one given in Fig. 3 ( -plane), so the axis tick marks are also given in units of [0.1 pc].

Figure 6 .
Figure 6.Total kinetic, thermal, and magnetic energy pools as functions of time  snr .The left picture stands for the central and the right one stands for the peripheral SN event.Due to the resolution down-scaling the magnetic energy curves have a small drop at  snr = 1500 yr.

Figure 7 .
Figure 7. Magenta line shows the part of the cluster core (domain) volume filled with the SN ejecta.Cyan dashed line is a ratio of the ejected mass to the total diffuse mass contained inside the cluster core volume.The left picture stands for the central SN and the right one stands for the peripheral SN event.

Figure 8 .
Figure 8. Distribution of the magnetic field in different bands of absolute value.The left picture stands for the case of centrally located SN and the right one stands for the peripheral SN event.It is clearly visible that during the propagation of the SNR shock through the cluster core volume the range of | | ∈ (3; 30) G becomes dominant in terms of volume.

Figure 9 .
Figure 9.Total mass of the gas in the cluster and its total mass loss rate through the cubic domain borders as functions of time  snr .The left picture stands for the centrally located SN and the right one is for the peripheral SN event.

Figure 10 .
Figure 10.The distribution of the hot ejecta of the central SN in the clusterthe material, consisting of >∼ 50 % ejecta and with  > 1 keV, is marked with the hatching.The colors represent the fraction of the SN ejecta in the material -pure ejecta is marked with red colour, while pure stellar wind material is marked with blue colour.Here,  snr = 900 yr.

Figure 11 .
Figure 11.Cluster thermal X-ray spectrum before the SN event (orange curve) and 300, 600 and 900 years after the event.The modelling is obtained for the central SN and the ejecta   = 0.1  .