Wind of a young massive star colliding with a supernova remnant shell

The fast stellar winds of massive stars, along with supernovae, determine the dynamics within the star-forming regions. Within a compact star cluster, counterpropagating supersonic MHD shock flows associated with winds and supernova remnants can provide favorable conditions for efficient Fermi I particle acceleration up to energies > 10 PeV over a short timescale of several hundred years. To model the nonthermal spectra of such systems it is necessary to know the complex structure of colliding supersonic flows. In this paper using the PLUTO code we study on a subparsec scale a 2D MHD model of the collision of a core-collapse supernova remnant with a magnetized wind of a hot rotating O-star. As a result the detailed high resolution (~ 10−4 pc) maps of density, magnetic field, and temperature during the the wind - supernova shell interaction are presented.


Introduction
The star-forming regions in galaxies are the birthplaces of massive stars. Observations show that the vast majority of massive stars do form in clustered manner [1,2]. More than 80% of O-stars in the Milky Way are believed to be a part of young massive star clusters (YMSCs) or loose associations and about 50% of the remaining population of field stars are identified as possible runaways [3,4,5]. YMSC usually contain several dozens of observable hot stars at different evolutionary stages [6]. Given the typical sizes of especially young (≤ 5 Myrs old) massive clusters being ∼ 1 pc and the typical number of hot stars within them, the average distances between the stars are ∼ 0.1 pc [7].
The dynamics of plasma flows inside YMSCs is completely determined by the fast and dense radiatively driven winds of hot massive stars, as well as by the explosions of Type II/Ibc supernovae (SNe) [8]. Stellar winds of isolated hot massive stars produce wind-blown bubbles which contain an extended region of highly supersonic wind that is far greater than the average distances between the stars in the considered compact YMSCs [9]. Therefore, the structure of interacting flows in the interstellar medium (ISM) of compact clusters turns out to be rather complicated. Numerous interacting supersonic magnetohydrodynamic (MHD) shock flows provide favorable conditions for the effective Fermi I acceleration of cosmic rays. Particularly, in compact YMSCs, such systems arise in a situation where a supernova remnant (SNR) shell collides with the wind's bow shock of a nearby massive star [10,11].
The increased frequency of SNe events in YMSCs is expected due to high concentrarion of massive stars with short lifespans. Core-collapse supernovae during the explosion releases a huge amount of kinetic energy which is typicaly ∼ 10 50 − 10 51 ergs [12]. Young massive stars  [13]. A star of M * > 25 M can shed more than a half of its mass before becoming a supernova, i.e. just over a few Myrs. The mechanical luminosity of these powerful winds varies widely from ∼ 10 35 and up to ∼ 10 38 erg s −1 strongly depending on the evolutionary stage of the central star [14]. During their lifetime, massive stars continuously evolve through a sequence of evolutionary stages (generally: main sequence → red or yellow supergiant → Wolf-Rayet phase) that are very different in their spectral properties. This sequence is mainly determined by the star's initial mass, metallicity, and rotation rate [15].
The complex MHD processes in close colliding stellar wind binaries which are sources of nonthermal emission were studied recently in [16,17,18]. The first large scale 2D hydrodynamic model of the collision between an SNR evolved into Sedov-Taylor phase and a wind-blown bubble was discussed by [19]. In this paper, using the PLUTO code [20,21], we consider on a subparsec scale a 2D MHD model of the collision of a core-collapse SNR shell with a magnetized wind of a hot rotating O-star. Compared to our previous paper [22] this model considers much shorter length scales ∼ 0.1 pc that are typical of interstellar distances in compact YMSCs. The rotating stellar wind parameters are now calculated based on the data of the Geneva stellar evolution code provided in S. Ekstrom et al. (2012). To initialize the SNR, a two-stage algorithm was applied [23], which makes it possible to track the earliest evolution of the SNR from a time of ∼ 0.1 yrs after the core collapse in accordance with self-similar solution of R. Chevalier [24]. In the energy balance equation, in addition to radiative losses, heating from photoionization is included [25]. Possible effects and efficiency of thermal conduction are briefly discussed.

Numerical simulations
We perform a 2D MHD simulation with the PLUTO code [20,21]. The dynamics of a magnetized plasma flow is described by the set of partial differential equations of ideal MHD including losses and heating by optically thin radiation. We use a second-order unsplit cell-centered finite volume method together with the HLLC approximate Riemann solver controlled by the standard Courant-Friedrichs-Levy (CFL) factor. Generally, the numerical schemes in PLUTO are stable up to CF L ∼ 1/N dim , where N dim = 2 is a number of dimensions considered. Here we adopt CF L = 0.2 in order to get lower timesteps for better stability against high Mach number flows. Radiative cooling and heating are taken into account using the net cooling function for a fully ionized medium (cf. [25]). The wind and ISM material are assumed to be at solar metallicity.
The simulations are performed in a cylindrical coordinates (R, z) that cover R ∈ [0, 0.4] pc and z ∈ [−0.1, 0.6] pc with a uniform grid of 1600×2800 cells, i.e. the effective resolution is 2.5×10 −4 pc. At the R = 0 border of the computational domain we use an axisymmetric boundary condition, while at the other three domain borders a free outflow condition is used. In the entire domain we initialize a magnetized ISM of constant temperature and density: T = 10 4 K, n = 0.5 cm −3 , B = −Be R , where B = 3.5 µG. At first, we simulate the ISM formed by the interaction of the winds of two massive stars: pre-supernova Wolf-Rayet (WR) star and O-type main sequence (MS) star. We obtain the wind properties from the analysis of the stellar evolutionary data provided by the Geneva group [15] for  Table 1). The supersonic stellar winds are injected into the domain via spheres (where wind profiles of density, velocity, and etc. are fixed in time) of 40 cells in radius placed at the coordinate origins for the Wolf-Rayet and at (0, 0.3) pc for the O-star. We use a wind compressed disk model of Bjorkman & Cassinelli [26] discussed in detail in [22]. The magnetic field strength at the stellar surface is set to 100 G for both WR and MS star [27], and the angular velocities of stars are as follows: ω M S = 3.85 × 10 −5 and ω W R = 2.00 × 10 −7 rad s −1 . On the considered scale, we assume the winds to be isothermal and their magnetic fields to be purely azimuthal [28]. Once the interaction of winds comes to a steady state we map the precalculated 1D self-similar solution of the SNR at the WR star position. The SNR expansion is launched using the method described in [29]. The mass of the ejecta and the energy release are M ej = 7.7M and E 0 = 10 51 erg.

Results
The general morphology, temperature, and magnetic field features of the flows during the collision are shown in Figures 1 − 3 at a time scale of ∼ 10 yrs. Compared to our previous result [22] here we have higher resolution and a length scale that is consistent with a typical raduis of YMSCs. In Figure 1   speed v snr ≈ 10 4 km s −1 . Some fraction of the kinetic energy released by the stellar wind and SN explosion thermalizes, heating the gas to X-ray-emitting temperatures. At the beginning of the collision, a wing-shaped region of compressed wind material is forming, this structure has temperatures of a few ×10 8 K and densities ∼ 1 cm −3 . The highest temperatures ∼ 10 9 K are observed at the very front of the shock wave having densities ∼ 0.01 cm −3 and propagating behind the star through the region of supersonic wind. At the regions of high compression This counterpropagating supersonic MHD shock flows inside compact young massive clusters can provide favorable conditions for efficient Fermi I particle acceleration up to energies > 10 PeV over a short timescale of a few hundred years [30]. To model the non-thermal spectra of such systems it is necessary to know this complex structure of colliding supersonic flows.
We do not include thermal conduction (TC) in this simulation, as the problem of anisotropic heat conduction efficiency in such complex collisionless astrophysical systems is poorly understood. The most recent measurements of the solar wind plasma (e.g., Parker Solar Probe) have shown that the heat fluxes reach saturation (for Knudsen numbers 0.1) well below the standard "free-streaming" flux value [31,32,33], q f s = 3/2n e k B T e v e , but there are still uncertainties. In terms of MHD anisotropic TC is not expected to change the general shape or morphology of bow shock, yet there are results indicating that it is capable of suppressing the development of hydrodynamical instabilities [34].