Numerical Study of Colliding Winds in Massive Stars

We run a numerical experiment ejecting stellar winds in a very massive binary system measuring the properties of the resulting colliding wind structure and accreted mass onto the companion under different conditions. Colliding massive binaries interact and create a colliding wind structure with a shape that depends on the momentum ratio, orbital motion, distance between the stars, and other factors. We run simulations of a static LBV-WR binary and in each simulation abruptly varying the mass loss rate of the LBV from the fiducial value. The modified wind front propagates and interacts with the previous colliding wind structure, and modifies its shape. We calculate the emitted X-ray from the interaction and investigate the proprieties of the new shape. We derive the mass accretion rate onto the secondary, and find that it depends on the momentum ratio of the winds. We then add orbital velocity that reduces the mass accretion rate, a similar behaviour as the analytical estimates based on modified Bondi–Hoyle–Lyttleton. Creating a large set of simulations like those presented here can allow constraining parameters for specific colliding wind binaries and derive their stellar parameters and orbital solution.


Introduction
Very massive stars are known to blow strong winds that intensify as these stars evolve off the main-sequence (MS) (e.g., [1,2]). When a binary system has two stars that both eject winds, the two winds collide in a process which may produce observational signatures across the spectrum. Colliding wind systems appear in O and B stars and WR binaries, O-O stars, and LBV binaries [3][4][5][6][7]. They are also present in symbiotic systems, where ejecta from the WD collides with the AGB wind and produces soft X-ray emission (e.g., [8][9][10]), and in low-mass binaries (e.g., [11]). Some colliding wind binaries present challenges to theory and common modelling, such as the Wolf-Rayet (WR) binary (or perhaps multiple) system Apep [12] that showed discrepancies in the observed wind velocities that suggest strongly anisotropic wind, never previously observed in colliding wind systems and challenging the understanding of WR winds. Presently colliding wind binary systems are being observed in the Magellanic clouds and in neighbouring galaxies (e.g., [13][14][15]). We will surely see more observations of such systems in the near future, and the field is ready for advanced theoretical interpretation with the aid of 3D hydrodynamic simulations.
The colliding wind structure (CWS) is classically defined as an infinitely long and thin axis-symmetric conical-like shape with thickness that varies with distance from the stars [5,6]. In the absence of orbital motion, its symmetry line is a line connecting the two stars. This structure is split by a contact discontinuity with two shock waves on both its sides. The shocked gas flows away along the sides of this conical-like structure. The shocked gas mainly resides between each of the two shocks and the contact discontinuity. The main parameter that determines the shape of the contact discontinuity is the momentum ratio of the two winds: whereṀ 2 and v 2 are the mass loss rate and the velocity of the wind of the secondary, andṀ 1 and v 1 are the same for the primary. The primary is considered to be the star with stronger wind, i.e., largerṀv. Thus, in WR-O binary systems the WR is the primary, and in LBV-WR systems the LBV is the primary. From Equation (1), essentially η < 1, and in some cases (that will be studied in this paper) η 1. The ability of the shocked gas to dissipate its thermal energy affects the shape of the CWS. Radiative shocks manage to radiate away the thermal energy of the compressed gas and therefore cool and become thin and dense. It is also more prone to instabilities, such as the Rayleigh-Taylor instability and the non-linear thin shell instability (NTSI, [16,17]). If energy cannot be lost to radiation its only way to cool is by adiabatic expansion, thus the CWS with such adiabatic shocks is wider and has higher temperatures.
The asymptotic semi-opening angle of the CWS is obtained from a non-linear differential equation with η being the only free parameter. It was calculated numerically by [18], and for intermediate values of η there is an approximation derived by [19]: which the authors state to be accurate to about 1%. It can be approximated, as done by [20] who approximated the contact discontinuity between the shocks as an hyperboloid, and gave analytic expressions for the width of the shocks (the distance from the contact discontinuity to the shock wave), on both sides. The expression for the shape of the CWS was later revisited a number of times (e.g., [21,22]). Indeed, the actual shape of the CWS is more complicated, taking other factors into account. A partial list of these factors includes: • The value of η is location dependent, due to the acceleration of the wind. • Clumps in the winds (small-scale inhomogeneities common in winds of massive stars; (e.g., [23][24][25]) make the CWS inhomogeneous. • Instabilities arise when the winds collide, as a result of cooling that will form clumps in the colliding wind region. • The radiation field of the secondary (the star with the weaker wind) may decelerate the primary incoming wind (radiative inhibition; [26,27]). • Mixing can result in a range of flow speeds instead of a single uniform one [21]. • Orbital motion completely alters the CWS shape, making it rotate and wind around itself. In eccentric orbits the distance between the stars and the orbital velocity is time dependent, causing periodic modulation in the CWS.
Simulations of colliding wind binaries usually focused on one system with its particular conditions. To date, there is no set of general simulations that ran over a range of each of the parameters and isolated the influence of each physical effect. While covering a large parameter space is challenging and requires vast computational resources, isolating specific effects while restricting to a fiducial set of parameters is possible.
When the two winds are uneven in their momentum, accretion of primary gas to the secondary can take place. Accretion onto a companion star is a process that is only beginning to be studied in numerical simulations. Akashi et al. [28] were the first to derive accretion in a colliding wind binary simulation in which the secondary has a strong wind. They found that for the massive binary system η Car mass is being accreted onto the companion close to periastron passage. A similar accretion process, though for a much longer period of time, was obtained for the O-O binary system HD 16674 [29].
Since the parameter space is very large it will have to be covered by large number of simulations. In this paper, as a first attempt, we run an experiment that quantifies accretion onto the secondary star, and the opening angle of the CWS.

Simulating Colliding Winds
We run the simulations using the hydrodynamic code FLASH, originally described by [30]. Our 3D Cartesian grid extends over x = (−8, 8) AU, y = (−8, 8) AU and z = (−4, 4) AU. The grid is centred around the secondary star, and the orbital plane is on the z = 0 plane. The grid has 5 levels of refinement. The largest cell has a side of 18 R and the smallest has a side of 1.1 R . The apex of the colliding winds has the second finest resolution, with resolution of 2.2 R . To solve the hydrodynamic equations, we use the FLASH version of the dimensionally split piecewise parabolic method (PPM) solver [31]. The internal parts of the stars are not modelled and do not participate in the simulations. The gravity is treated as two-point gravity since self-gravity is not important. We include radiative cooling from [32]. The simulation includes radiation transfer with multigroup diffusion as implemented by the FLASH code. The simulations and post-processing described here were performed on the Ariel University GALAXY HPC cluster, for ∼10 6 core-hours.
For most of the experiments presented here, the binary system remains stationary, and we do not include the orbital motion as the experiment aims to isolate the effects of the mass loss rate. We place the two stars at a distance of a = 6 AU, and allow each star to accelerate their own wind. The mass loss rates of the primary and the secondary arė respectively. The winds are ejected from a spherical layer above the star. The winds are radiatively accelerated according to the CAK model [33] force multiplier, with an acceleration that corresponds to β = 1 for the primary and β = 0.8 for the secondary. These values are motivated by a number of studies (e.g., [34,35]) and were chosen for simplicity. The acceleration also occurs outside the ejection zone. The terminal velocities are v 1,∞ = 500 km s −1 and v 2,∞ = 3000 km s −1 . This gives a momentum ratio of η = 0.2. We name this set of parameters as Run 1, and treat it as our fiducial run. Figure 1 shows the density of the wind for the simulations we run. The stronger the primary wind, the smaller the opening angle of the colliding wind structure. The winds blow for 68.5 days until the colliding wind structure extends to the edges of the grid, and then the experiment continues to the next step in which we turn to run the other simulations (runs listed below as 2-5). In Runs 2-5, we make a sudden transition of the mass loss rate of the primary, and by doing so lower the value of the momentum ratio η. As it takes time for the wind to propagate, the change is not immediate, but rather takes time to establish. Table 1 lists the parameters of the simulations we run, as well as a summary of the results. Figure 2 shows density slices for some of our simulations.   With the increase in the mass loss rate of the primary wind, the secondary wind and radiation cannot overcome the inflowing primary wind. Therefore, the secondary accretes mass from the primary wind. Let us describe the accretion process and the way we measure the accreted mass.
As the simulations run, the colliding wind experiences instabilities, especially in the side facing the secondary, that cools faster, and high-density clumps are formed. These clumps approach the injection zone of the secondary wind and even reach the cells of the secondary itself. A number of approaches are examined for the best numerical implementation for accretion [36], and we here briefly describe the one we use. To calculate the accreted mass we check whenever the actual density ρ a,cell in a cell in the injection zone supersedes the expected undisturbed value of ρ u,cell (expected as a result of the ejection of radiation-driven wind), and count the extra mass as accreted. Namely, the accreted mass is where V cell is the volume of the cell. The contributions from all cells in the injection zone are then summed up to obtain the total mass accreted for that time step. As the simulation time steps are much shorter than the duration of the accretion, this procedure brings accurate results. This process was tested in [29,36], where more details can be found. Figure 3 shows the mass accretion rate (lower panel) and the accreted mass as a function of time for Runs 2-5 (for Run 1 there is no accretion).

Results
Our runs continue for another 130 days, which is more than enough time for the new outflow to cover the entire grid and to obtain enough data for calculating average accretion rate, and X-ray luminosity, as we describe below. In one of the runs, Run 4, we include orbital motion in a circular orbit. This run can be compared with Run 3 that does not include orbital motion, as for both runs η = 0.005. A 3D density contour plot of this run is shown in Figure 4. The curving of the CWS is clearly seen. The orbital velocity in Run 4 is about a quarter of the wind velocity of the primary, namely, it has a moderate effect. The volume inside the cavity of the secondary wind does not include an undisturbed secondary wind but rather clumps from the colliding winds.
The X-ray emission is computed by interpolating from emission tables according to the implementation of the yt-project [37]. Figure 5 shows the 2-10 KeV X-ray emission on a slice of the simulation grid. For Run 1, the emission arrives from the colliding wind structure. For stronger mass loss rates in the rest of the runs, the emission arrives from shocks close to the secondary, and highlights the attempts of the secondary wind to blow against the collapsing clumps of the CWS around the secondary. As the figure shows, the smaller the value of η, the more concentrated the emission on regions closer to the secondary. Figure 6 shows the X-ray luminosity L x emitted from the entire simulation grid, as a function of time. We do not include the absorption, which can be large and can greatly attenuate the observed flux by a large factor, and is also line-of-sight dependent. The enhancement in the mass loss rate of the primary creates a dense front that propagates and hits the colliding wind structure. As a result of this collision, the X-ray emission increases and a 10 days peak appears. After this time, which is the time that it takes for the intense wind to replace the quiescence wind (as in Run 1), the X-ray emission reaches a new quasi-stable level. The fluctuations in the level are the result of the interaction between the secondary wind that struggles to push against the primary wind.   The total emitted X-ray from the entire simulation grid (absorption is not included). Note the different scaling for some of the runs. The initial peak is the result of interaction of the enhanced wind with the previous CWS from the fiducial run (Run 1). The average X-ray luminosity (excluding the initial peak) is given in Table 1 and Figure 7. In Run 4 that included the orbital motion, the luminosity is stronger than in Run 3, which had the same parameters without orbital motion.
As explained above, we measure the mass accreted onto the secondary. In Run 1, no mass is accreted, but in all other runs, accretion onto the secondary takes place, either continuously or intermittently. Figure 7 shows the average mass accretion rate onto the secondary star (blue circles), and average X-ray emission luminosity (black squares). The grey symbols are for Run 4 that includes orbital motion. We find that the amount of mass accreted onto the secondary increases with decreasing η, but not in a linear way. The accretion dependency on η will be investigated in a future paper that will include more simulations. In this accretion domain, the X-ray emission is generally stronger for smaller η, however, this trend breaks for Run 6 as the primary wind is so strong to the extent that the secondary wind almost cannot blow and therefore create shocks that emit X-rays.  Figure 7. Average mass accretion rate onto the secondary star (blue circles), and average X-ray emission luminosity (black squares). The average X-ray emission excludes the initial X-ray peak which is the result of interaction of the enhanced wind with the previous CWS from the fiducial run (Run 1). The orange symbols are for Run 4 that includes orbital motion.
Previous scaling relations of η with L x have been proposed before based on theoretical arguments (e.g., [5]), and simulations (e.g., [22]). These were mostly obtained for adiabatic wind collisions. However, previous scaling relations are not applicable for very low η as we examine here, especially when accretion takes place. Scaling relations as above use the total X-ray luminosity, which is a volume integrated quantity. They do not address the question of where the X-ray emission comes from, which is a different place than for higher η values, as explained above and shown in Figure 5.
We use image processing methods to fit the shape line of the wind, an imaginary line dividing the turbulent flow into two domains, and the semi-opening angle, an asymptotic angle that describes the CWS. We use Canny edge detector [38], a multi-stage operation filter based on the derivative of a Gaussian, in order to compute the intensity of the gradients and by that determine a line that traces the shape of the wind. To reduce noise in the resulting image, we use morphological operations (dilation and erosion). Then, based on the shape of the wind and the location of the primary and secondary, we can derive the semi-opening angle of the wind. To obtain the optimal rays that form the opening angle, we calculate the Euclidean distance of each point in the line to the shape and sum over all the distances. We then minimise the function using the quasi-Newton BFGS method to obtain the best solution. Figure 8 shows the resulting opening angle for our six runs.
We plot the wind semi-opening angle as a function of η in Figure 9. For most of our simulations, except Run 6, the simulations obtain larger opening angles than Equation (2). As we have an accelerating wind, η is different at different locations, and this might affect the calculation. However, the more significant effect is the formation of instabilities, and requires the treatment we described above. For Run 6, we obtain a smaller semiopening angle than the result from Equation (2). The equation is actually not applicable for such a small value of η. In this run, the strong mass loss from the primary chokes the secondary wind almost completely. Therefore, the semi-opening angle does not really describe the same problem of colliding winds of two stars, but rather the collision of the primary wind with itself after being focused by the secondary and its gravitational field.

Summary and Discussion
We explore the possibility of studying physical effects in colliding winds by simulating them in a systematic way. We start by quantifying a parameter that is not often discussedaccretion onto the companion star. We find that the mass accretion rate depends on the momentum ratio of the winds, and generally increases with decreasing η. The instabilities have the velocity component in the direction of the flow but also have the acceleration of the secondary. Depending on the parameters, some of the clumps can be pulled towards the secondary and be accreted. In very low eta, the accretion is directly onto the secondary as there is no colliding wind structure. Our results also hint at saturation in the accreted mass for very low η. The actual quantification of this dependency requires more simulations, which will be carried out in the near future.
We chose certain values for the mass loss rates and other parameters with the aim of describing a general system. These parameters do not represent the typical LBV-WR system. In future works, it will be necessary to check whether changing them makes a significant change in the results.
The accreted wind can have more effects, especially at high accretion rates. One of them is the reduction in the effective temperature of the secondary, as shown in [39]. At very high rates, the accreted gas might interact with the stellar envelope and modify the evolutionary path of the star. We do not model this in the present paper.
Our simulations also show that orbital motion reduces the mass accretion rate. This is in agreement with analytical estimates based on modified Bondi-Hoyle-Lyttleton. We note that the accretion is not Bondi-Hoyle-Lyttleton like and most of the accreted gas, it actually arrives from the side facing the primary and not from behind.
We also calculate the emitted 2-10 KeV X-ray luminosity from the colliding wind systems. The emission shows a first large peak at the time the enhanced wind of the primary reaches the pre-existing colliding wind structure. After that, there is a quiescence level of emission, associated with the new meta-stable colliding wind structure. The actual observed X-ray flux from such a system will depend on the absorption in the line of sight, and is expected to look different from different directions.
We used an advanced method to model the the colliding wind (asymptotic) semiopening angle. While this practice is relevant for stationary stars, for orbiting stars it is not very well defined as the orbital motion curves the CWS, and the asymptotic semi-opening angle, are not well defined. We found that for most cases the angle is larger than the theoretical expression predicts, except for very small η in which the assumptions of two winds that collide break and the primary wind in fact collides with itself.
The parameter space we chose here does not match particular systems. The aim of this work is to eventually be able to minimise the number of attempts required when performing simulations of real systems. For that, a better cover of the parameter space and especially testing simulations with different values of eccentricity is still needed. This present range of large mass accretion rates and very small values of η can be helpful in simulating giant eruptions.
Running more simulations and modifying other parameters would allow a better understanding of the changes in the observables resulting from each of them. We expect such a process will be useful for studying specific massive binary wind systems, by reducing the number of attempts needed to simulate each system in order to match observations. This, in turn, will provide better constraints on massive colliding wind binaries' physical properties.