Simulating stellar winds in AMUSE
Edwin van der Helm Martha I. Saladino Simon Portegies Zwart Onno Pols

    

    
    

    
        Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, The Netherlands Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands
        


    

    
        We present , a module that provides multiple methods of simulating stellar winds using smoothed particle hydrodynamics codes (SPH) within the astrophysical multipurpose software environment () framework.
        
        The module currently includes three ways of simulating stellar winds:
        With the simple wind mode, we create SPH wind particles in a spherically symmetric shell for which the inner boundary is located at the radius of the star.
        We inject the wind particles with a velocity equal to their terminal velocity.
        The accelerating wind mode is similar, but with this method particles can be injected with a lower initial velocity than the terminal velocity and they are accelerated away from the star according to an acceleration function.
        With the heating wind mode, SPH particles are created with zero initial velocity with respect to the star, but instead wind particles are given an internal energy based on the integrated mechanical luminosity of the star.
        This mode is designed to be used on longer timescales and larger spatial scales compared to the other two modes and assumes that the star is embedded in a gas cloud.
        
        We present a number of tests and compare the results and performance of the different methods.
        For fast winds, we find that both the simple and accelerating mode can reproduce the desired velocity, density and temperature profiles.
        For slow winds, the simple wind mode is insufficient due to dominant hydrodynamical effects that change the wind velocities.
        The accelerating mode, with additional options to account for these hydrodynamical effects, can still reproduce the desired wind profiles.
        We test the heating mode by simulating both a normal wind and a supernova explosion of a single star in a uniform density medium.
        The stellar wind simulation results matches the analytical solution for an expanding wind bubble.
        The supernova simulation gives qualitatively correct results, but the simulated bubble expands faster than the analytical solution predicts.
        We conclude with an example of a triple star system which includes the colliding winds of all three stars.
        

    
    Simulating stellar winds in AMUSE
    
        Edwin van der Helm 1 Martha I. Saladino 2,1 Simon Portegies Zwart 1 Onno Pols 2
        
    Received October 2, 2017; Accepted March 23, 2019
=====================================================================================================





§ INTRODUCTION

        
        Stars lose mass through stellar winds during various stages of their evolution <cit.>.
        These winds can affect the gas near the star, creating lower density bubbles <cit.> and regulating star formation <cit.>.
        If a binary companion is present, accretion of the stellar wind material can also affect the evolution of that companion <cit.>.

        
        The two most important parameters of the stellar wind are the mass-loss rate, , and the terminal wind velocity, , which determine the effect of the wind on the environment.
        Based on these parameters, stellar winds can be broadly divided into three categories <cit.>:
        1) Cool main-sequence stars like the Sun have winds with very low mass-loss rates (10^-14 /yr) that are thermally or gas pressure driven.
        2) Cool giants and super giants have slow (5  30 km/s) high mass-loss rate (10^-7 10^-5 /yr) winds driven mainly by radiation pressure on dust <cit.>.
        3) Hot luminous stars have fast (200  3000 km/s) high mass-loss rate (10^-7 10^-5 /yr) line driven winds <cit.>.
        The second and third category have the highest kinetic output and therefore have the most pronounced effect on the stellar environment (not including cumulative effects).

        
        To simulate stellar winds in detail, a combination of hydrodynamics, radiative transfer, dust formation and chemical abundances is required.
        Such simulations have been done for many years although they are extremely computationally expensive.
        In most cases simulations are limited to 1D or 2D models <cit.>.
        To investigate the net effect of the stellar wind on the environment, it is often sufficient to simulate the stellar wind using only hydrodynamics    <cit.>.
        For larger scale simulations, stellar wind feedback is often included using a sub-grid model as it can influence star formation and launch galactic winds <cit.>.

        
        For all these simulations, the astrophysical multipurpose software environment <cit.> can be useful.
        It provides a uniform interface for many types of simulations with a large and growing set of simulation codes.
        The consistent python interface makes it possible to quickly set up a scientific simulation and easily exchange different parts of these simulations.
        While stellar winds have been simulated before using  <cit.>, a consistent and properly tested module was still missing.

        
        The  code presented in this paper can be combined with the SPH (Smoothed particle hydrodynamics), N-body, stellar evolution and (with some additional work) radiative transfer codes that are already available.
        We describe the  code and explain the different modes in which it can be used (Section <ref>).
        In Section <ref> we present a series of tests in which we compare the results from the different modes with theoretical expectations and previous work.
        We conclude in Section <ref> with an exposition of some ongoing research projects using this code and ideas for further use.



§ METHODS

        

        
        The goal of  is to create gas particles that represent the stellar wind from one or more stars.
        The code requires a number of stars, represented by  particles[A particle set is the fundamental data structure in . It is an array of particles (stars, SPH particles etc) which contain information to control the data. Each element (particle) of the particle set has certain attributes, such as mass, position, velocity, etc.], with stellar properties that can be derived from observations or stellar evolution simulations.
        Using this, SPH particles are created with the appropriate wind properties in an initially spherically symmetric shell with inner boundary at the radius of the star.
        The number of SPH particles is computed according to the mass-loss rate associated with the star undergoing mass loss and the predefined SPH particle mass, M_SPH.
        These particles can be added to any SPH code in  which simulates the hydrodynamics of the wind.

        
        Creating the SPH particles is only one step in the simulations for which  is used.
        Following the goal of the  framework, the other parts of the simulations are handled by specialized interchangeable codes.
        For the hydrodynamics, SPH codes such as  <cit.> and  <cit.> can be used.
        In many applications, the stars move, for which a large number of N-body codes are available.
        To couple the stellar dynamics to the hydrodynamics gravitationally,  <cit.> is available.
        The stellar properties on which the wind is based will generally be calculated using a stellar evolution code.
        Both parametrized <cit.> and Henyey type <cit.> stellar evolution codes are available in .
        Any or all of these codes can be combined with  to set up a wide variety of simulations (see Section <ref>). For more information about the codes available within  and examples of how to couple them, we refer the reader to <cit.>.

    

 §.§ Calculating stellar wind properties
 
        
        To simulate the stellar winds, the stellar parameters (mass, radius, temperature and position) and wind parameters (mass-loss rate, initial and terminal wind velocity) are required.
        All these parameters can simply be set directly, however some of them can be derived directly from stellar evolution codes available in .

        The  module includes user-friendly routines to derive some of the stellar parameters such as stellar mass, mass-loss rate, stellar radius and effective temperature from one of the stellar evolution codes within .
        However, the terminal wind velocity, , is not calculated by any code currently in .
        Determining  requires detailed and computationally expensive stellar wind simulations which include radiative transfer.
        For this reason, in order to compute the terminal velocities of hot stars, we provide within  a formula that has been fitted to observations of hot stars <cit.>, and which, they claim, is valid for these stars within 20%.
        The terminal velocity of the wind is given by:
        
        
    = C(T_*) , 
    where,
    
                C(T_*)     =
                
                    1     T_* ≤ 10 000 K, 
    
                    1.4     10 000 K < T_* < 21 000 K, 
    
                    2.65     T_* ≥ 21 000 K, 
    
        = √(2 g_* R_*  ( 1 - Γ )), 
    
                g_*     = G M_*/R_*^2, 
    Γ    = 7.66 · 10^-5σ_e L_* / /M_* / , 
    σ_e     = 0.398  1 + I_HeY/1 + 4 Y,

        
        where  is the photospheric escape velocity (similar to the escape velocity  with a correction term for Thomson scattering), G is the gravitational constant, M_*, R_*, L_* and T_* are the mass, radius, luminosity and effective temperature of the star respectively, Γ is the ratio of radiative Thomson acceleration to gravitational acceleration, σ_e is the Thomson absorption coefficient, Y is the Helium fraction and I_He is the number of electrons per Helium nucleus (in this paper we use default values of I_He = 2 and Y = 0.25).
        For cooler stars, ≈ and this formula is still applicable (Kudritzki, private communication).

    

 §.§ Simple wind

        
        
        Within , there are currently three wind modes available.
        The simplest mode creates a spherical shell of particles around the star with radial velocity, v(r) = and initial temperature equal to the effective temperature of the star.
        While this may sound simplistic, a similar setup has been used effectively for a number of scientific problems <cit.> and it serves as a starting point for the two other modes described in Sections <ref>  and <ref>.
        When the gravitational attraction of the star on the wind is included in the simulation, however, this will not result in the desired terminal wind velocity.
        We therefore release the wind with a larger velocity v(r) = √(^2 + (r)^2) where (r) = √(2 G M_* / r) is the local escape velocity at the initial particle radius, r.
        We calculate this new velocity for each particle because (r) can vary within the thin shell in which we create the particles.

        
        We set the outer radius () of the shell of new particles at the radius that the innermost part of the previously released shell (R_*) should have reached in the elapsed simulation time δ t (see Appendix <ref>).
        We scale the particle positions within the shell to follow the density profile matching the velocity profile as described in Appendix <ref>.

    

  §.§.§ SPH and initial distributions

        
        

        
        SPH is a method to solve the dynamics of a fluid by approximating it with a set of discrete particles <cit.>.
        Each particle has both a mass and a density, where the density is calculated using the distance to, and mass of, other particles that are nearby.
        To determine which other particles are taken into account (how nearby they have to be) a kernel function[In the SPH code used in this paper we use the spline kernel.] and a particle smoothing length (h) are used.
        In all modes of , the SPH particle mass is required to be fixed and the same for each particle.
        The smoothing length, h, is set adaptively by fixing the number of neighbouring particles that fall within one h <cit.>.

        
        Creating an initial distribution of SPH particle positions is not trivial <cit.>.
        Randomly distributed positions are clumpy which can introduce shot noise that can affect the entire simulation.
        A better alternative is to have more regular spaces between particles, for instance a distribution that follows a grid.
        However, a regular grid tends to introduce preferred directions in the simulation that can affect the overall results.
        To solve this, it is common to start with either a random or grid distribution and let the system evolve (relax) to a steady state where the positions are regularly spaced without preferred directions <cit.>.
        While some form of relaxation is preferred for simulations where all particles are created at once, for continuous particle creation like we describe here, this is not generally required.

        
        In  we implement two methods in which wind particles can be initially distributed.
        One is a random distribution and the second follows a uniform grid.
        We present an example of both in Figure <ref>.
        The random initial distribution (top panel) is available so that users can investigate if it has advantages for their simulations.
        In this case, a shell with uniform density is created and then the radii are scaled to ensure the correct density profile.
        In the other option we have included (bottom panel), each new shell is created by cutting it out from a cube with positions following a uniform grid.
        The number of particles in this shell is generally not exactly the desired number of particles, N_desired = δ t ·Ṁ/M_SPH.
        We therefore remove a number of randomly selected particles from the grid (typically 30%) to ensure the correct number of SPH particles.
        The grid can be randomly rotated each time a new shell is generated to avoid introducing preferred directions into the resulting wind.
        The positions of the particles are also radially scaled to ensure that the desired density profile is achieved.
        This is the cause of the curved appearance in the grid in Figure <ref>.

        
        There are many more ways to create initial particle distributions.
        A good overview of the different methods and their advantages can be found in <cit.>.
        Our method is a mix between the `stretched lattice' and the `shell' methods described there.
        The reason we do not use the more advanced methods described there is that they would require some form of computationally expensive relaxation for every new shell.
        This is a common issue with continuous particle creation methods.
        If the current methods are found to be unsatisfactory for a specific simulation, the code is set up in a modular way so adding a new particle distribution method is relatively easy.
        The uniform grid with random rotation is the default option used throughout this paper.
        However, due to the small number of particles in a single shell, the difference between this option and a random distribution is negligible for all the tests in Section <ref>.

    

 §.§ Accelerating wind

        

        

        

        
        Near the surface of the star, usually within a few stellar radii, the wind is accelerated to the terminal wind velocity.
        In the accelerating wind mode, the wind particles are created in the same way as in the simple wind mode, but with a lower velocity, v <.
        All particles near the star are artificially accelerated in such a way that the wind follows a predefined velocity profile.

        
        The artificial acceleration is implemented using  <cit.>.
        Originally,  was designed to couple multiple gravitational codes.
        In this method, each code is evolved separately for a short, predefined timestep[
            This should not be confused with the internal timesteps for each code which may be variable, meaning particles have different timesteps depending on conditions such as local particle density.].
        The mutual gravitational effect is included by  using a kick-drift-kick scheme <cit.>.
        This method can also be used to gravitationally couple a pure N-body code with an SPH code, or to apply a gravitational potential to the particles in one or more codes.
        In , we use  by including an artificial potential near the star, and then use the same kick-drift-kick scheme to ensure a smooth acceleration of the wind particles.

        
        In Figure <ref> and Table <ref> we present the acceleration functions (sometimes referred to as acceleration laws) currently implemented in .
        For a given velocity (or acceleration) profile, all required quantities are calculated following the equations in Appendix <ref>.
        The constant velocity function is similar to the simple wind mode in that when the wind particles are created, they already have the terminal velocity.
        However, as noted below, when used in the accelerating mode we can add extra terms to compensate for the gravity of the star, as well as the pressure of the gas on the wind.
        These extra accelerating terms are added after the particles have been created, which is not possible with the simple wind mode.
        In this way, we guarantee the desired constant velocity profile.
        The logistic and agb functions provide a fit to the time-averaged behaviour of dynamical models of asymptotic giant branch (AGB) winds from <cit.>.
        These winds exhibit a specific acceleration zone, the location of which can be chosen with the parameters  and either s or α (for the logistic and agb function, respectively).
        These parameters determine the center and the width of the acceleration zone.
        The default values = 3 and  s = α = 10 are chosen to fit the dynamic models.
        The beta_law function, which was derived using a combination of observations and theoretical wind models, was taken from <cit.> and <cit.>.
        The β parameter indicates the steepness of the acceleration curve and is often derived from observations.
        The example values β = 0.8 and β = 2 are typical for hot and cool stars respectively.
        The rsquared and delayed_rsquared functions can be used as rough approximations to the wind profiles of hot stars and cool giants, respectively.
        They have the advantage of being computationally faster than the beta-law, agb and logistic functions.
        In the delayed_rsquared model the parameter  (with a default value of 2) sets the lower boundary of the acceleration zone.
        The initial velocity , which is used in all functions except for constant_velocity, is of the order of a few km/s due to microturbulence in the stellar atmosphere where the material is launched.
        We note that low values of  result in high densities which lead to slow simulations, so in many cases a higher value of  can be used as an approximation.
        In addition to these predefined functions, new user defined velocity functions can easily be incorporated.

        
        When the gravity of the star is included, an additional acceleration term can be added to compensate for it and ensure that the wind particles follow the chosen velocity profile.
        The gas pressure can also exert an acceleration on the wind.
        We therefore provide the option to subtract the expected gas pressure acceleration (see Appendix <ref>) from the applied artificial acceleration.
        If we do not include a hydrodynamical simulation of the stellar interior, unphysical boundary effects near the surface of the star can influence the wind evolution or even prevent the wind from being launched.
        We therefore provide the option to create a “staging shell” near the star, generally at least twice the thickness of the newly created shells.
        Within this shell, the accelerations are chosen in such a way that the desired velocities are enforced regardless of the gas dynamics.
        This shell then provides a better boundary condition for the rest of the simulation.

        
        For many simulations using the simple and accelerating wind we can start the simulation with a vacuum around the star into which the wind particles are released.
        However, this can lead to an extra acceleration at the outer radius as the vacuum does not exert any pressure on the outermost particles.
        This can be problematic, especially for slow winds, where this spurious acceleration can significantly increase the velocities.
        We therefore include a function in  that creates an initial set of SPH particles following the desired temperature, density and velocity profiles up to a given radius.
        This function uses the same initial grid distribution described in Section <ref>.
        Since the whole grid is created at once, without random rotations between different shells, this can introduce preferential directions.
        We advise that any scientific measurements are started after all these particles have left the area of interest.

        
        When particles are created, we ensure that they follow the desired density profile by solving equation <ref> for each particle.
        For most acceleration functions, this equation has to be solved numerically, which can severely slow down the simulation (see Section <ref>).
        We therefore include the option to define a critical timestep, .
        When new wind particles are created, δ t, which determines the thickness of the shell of new particles, is compared to .
        If δ t <  it means that the new wind particles have not reached the accelerating region yet.
        For this reason, the acceleration function is approximated by the constant velocity function, for which equation <ref> is solved analytically.
        This approximation is only valid for acceleration functions for which the velocity near the star is close to constant, like the logistic function, not for acceleration functions with a large acceleration near the stellar surface, like the beta_law function (see Figure <ref>).

    

 §.§ Heating wind

        
        
        The third wind mode is based on the method used in <cit.> and is designed for use in large scale simulations, e.g. embedded star clusters.
        For these simulations, the main effect of the wind is that it adds mass and energy to the surrounding gas, therefore this mode cannot be used for a star in a vacuum.
        Studying the detailed kinematics of the wind near the star is not the goal of these simulations and therefore a simpler approximation of the wind interaction is used.
        The advantage of this approximate approach is that it can be used at far lower resolution (longer timesteps and higher SPH particle mass) which greatly speeds up the simulations.
        If particles were created with a high velocity, small timesteps would be required to completely sample the particle trajectory and interactions with other particles.

        
        The basic idea of the heating wind mode is that new wind particles do not have an initial velocity relative to the star.
        Instead they have an internal energy, u, which corresponds to the mechanical energy () of the accumulated wind, defined as,
        
    = ∫_t_0^t_1(t)   dt, 
    (t)     = 1/2(t) (t)^2,

        where  is the instantaneous mechanical luminosity and t_0 and t_1 are the previous and current wind release time respectively.
        The integral is numerically approximated in  during the simulation.
        The internal energy of the new particles is set to
        
    u = f_fb/Δ M_*,

        where Δ M_* is the mass lost and f_fb is the feedback efficiency parameter that accounts for radiative losses.
        Typical values for these parameters can be found in the examples shown in Sections <ref> and <ref>.

        
        As discussed in <cit.>, this method of creating particles with appropriate internal energy can also be used to simulate a supernova.
        If a star goes supernova, the calculated mechanical energy is ignored, and instead 10^51 erg of energy is divided over the newly created particles.
        It should be noted that the injection of so much energy in the surrounding gas will cause the gas to evolve very rapidly, which can lead to time-stepping artefacts <cit.>.
        One way to prevent this is to use a very small timestep (preliminary tests suggest 10 yr) shortly after a supernova.



§ APPLICATION

        
        


        
        To ensure that  performs as expected, we run test simulations with different initial conditions and wind modes.
        In this section we present the results of these tests.
        The tests in Sections <ref> and <ref> are simulations of processes that happen close to the star.
        Therefore only the simple and accelerating wind modes are applicable.
        The tests in Section <ref> and <ref> are large scale simulations of the interaction between the stellar wind and the gas in which the star is embedded.
        For this type of simulation the heating wind mode is applicable.
        The final test in Section <ref>, where we couple stellar dynamics, hydrodynamics and stellar winds, shows the power of  within  by simulating the colliding winds from a stable triple star system.
        The initial distribution of the wind particles is the one based on a grid for all the models in these tests.
        The particular parameters used are described accordingly in the following subsections.

        
        The  module is designed to couple different parts of a simulation in .
        Testing the code requires the use of other simulation codes, which have their own parameters to be set.
        Here we describe the general method and general parameters used in our test models.
        In Figure <ref> we present a flow diagram to illustrate the codes and the relationships between them.
        Note that the codes and coupling strategies used here are merely an example and should be modified in order to be suitable for any specific application.
        To simulate the gas we use the SPH code  with an adiabatic equation of state and artificial viscosity parameters α = 0.5 and β=1.0 <cit.>.
        Self-gravity of the gas is only included in the tests which do not require periodic boundary conditions, i.e. those tests where the simple or accelerating mode are used.
        For simulating the gravitational attraction of the star on the gas we need , which is also used in the accelerating wind mode (Section <ref>).
        The  code requires an additional code for calculating the gravitational force.
        When we simulate only a single star that does not evolve dynamically, we use [The  code, developed by N. de Vries, is an unpublished gpu-enabled code in  that can calculate the gravitational force of one set of particles on another set of particles.
                      It is ideal for the gravitational coupling between particles in different codes via .].
        When we simulate multiple stars that evolve dynamically (Section <ref>), we use the N-body code  <cit.>.
        For each simulation, we also need to define a number of integration timescales, such as the  timestep (t_br), the (maximum) internal timestep () of the SPH and N-body codes and the wind release timestep (t_wind), as well as the end time (t_end) of the simulation.
        The choice of these timesteps depends on the problem and the type of simulation.
        For the bridge leap-frog algorithm to work, we should set t_br≥ and the wind code requires t_wind≥ t_br.
        It is also a good idea to ensure that larger timescales are integer multiples of smaller timescales.
        For the simulations in this paper, we only define t_wind and choose t_wind = 2 t_br = 4.
        For the hydro simulation we also need to set the SPH particle mass (M_SPH).

    

 §.§ Fast winds

        

        

        

        

        

        
        In this section we present the results of a set of simulations using the simple and accelerating wind modes.
        We simulate the wind from a single, hot, massive, luminous star, for which we present the parameters in Table <ref>.
        Note that these values were not chosen using a specific stellar model and this test should only be considered as an example of the use of .
        The initial wind velocity, = 100 km/s, is based on numerical considerations.
        As we shall see in Section <ref>, slow (subsonic[If the wind speed near the star is lower than the local sound speed, the wind is called subsonic. On the other hand if the wind speed is higher than the local sound speed, i.e. past the 'sonic point', the wind is called supersonic.]) winds are more complex to simulate and for many applications using a higher  is sufficient.

        

        In Figure <ref> we show the velocity, density and temperature profiles for the simple wind mode and two accelerating functions: the beta_law and logistic function.
        In all cases the velocity profiles follow the desired analytical velocity curve.
        For the simple wind mode, the density and temperature profiles also follow the desired curve, but with more scatter.
        For the accelerating wind curves, we see that in regions with high acceleration the densities and temperatures in the simulation are too high.
        This is a result of the low resolution in combination with the way densities are calculated in SPH, using a kernel function that 'smears out' these variables.
        We show in Figure <ref> that for a higher resolution (smaller M), the desired density and temperature curve are recovered.
        Note that the logistic acceleration function is not a good representation for the velocity profile of a hot massive luminous star.
        This example was chosen to illustrate the discrepancies that can potentially occur.
        For any scientific application of this code, a detailed convergence test for the selected setup will still be required.

        
        In addition to being accurate, it is also important that a simulation code is fast.
        In Figure <ref> we present the time spent in different parts of the simulation code () divided by the total cpu (or wall-clock) time () as a function of resolution[These simulations where all performed on the same desktop computer using a 4-core Intel Xeon E5507 CPU.].
        In the top panel we see that when using the simple wind mode, the time spent in the  code is less than 1% when using more than 10^4 particles.
        Most of the simulation time is therefore spent in the SPH code itself, which is what we want.
        When we use the accelerating wind mode however (middle panel), the particle creation becomes a major bottleneck because numerically solving equation <ref> is slow.
        To speed up the simulation, we have included the option to approximate the acceleration function with a constant velocity when particles are created near the star by defining a critical timestep (, see Section <ref>).
        In the bottom panel we see that by using this approximation, the time spent in  reduces to < 1% for > 10^4 particles.

    

 §.§ Slow wind

        

        

        

        

        

        

        
        For the slow wind test, we simulate the wind from a single, cool, giant star, for which we present the parameters in Table <ref>.
        The values of these parameters are not computed with a stellar evolution code, but they correspond to typical values for AGB stars.
        Part of the stellar wind is subsonic and therefore hydrodynamical effects are no longer negligible, unlike in the fast wind test.

        
        In Figure <ref> we present the velocity profiles of simulations using the simple wind mode and varying  at t = 5 days.
        We have not included the stellar gravity in these simulations, therefore the escape speed is not a relevant factor.
        However, when the wind speed is near or below the sound speed, the gas pressure gradient dominates and affects the terminal wind velocity.
        For very low wind speeds, this can cause wind particles to move inside the star, which is unphysical.
        We therefore conclude that the simple wind mode is not reliable for slow winds.

        
        In Figure <ref> we present the results of simulations using the accelerating wind mode.
        In these simulations we have used the staging shell option (Section <ref>) that enforces the correct particle velocity for all SPH particles with radius r < 1.1 R_*.
        We have also used the option to subtract the expected gas pressure acceleration from the acceleration to ensure that the particles follow the desired velocity profile.
        To avoid spurious acceleration of the outer particles, we start the simulation after initially creating a sphere of particles following the desired velocity and density profiles throughout the simulation area (up to r = 1500).

        
        For the constant, rsquared and beta_law velocity profiles, the particles follow the desired velocity profiles.
        For the logistic velocity profile, where particles are subsonic for a longer time, the simulated velocity profile deviates from the desired velocity profile in the subsonic region.
        However, the corrections described above ensure that the resulting velocities after particles pass the sonic point follow the desired velocity profile.
        To see how this deviation will affect the results of simulations where the subsonic region is of interest, we have compared this velocity profile with detailed velocity profiles of AGB stars <cit.>.
        In these profiles, the velocities in the subsonic region oscillate due to stellar pulsations and do not follow the simple acceleration functions we have used here.
        We therefore advise caution when interpreting results of simulations in the subsonic region.

        
        Similar discrepancies in density and temperature as seen for the fast wind test in Figure <ref> are also present for the slow wind test in Figure <ref>.
        In Figure <ref> we present the results of a resolution test for the slow wind test.
        We see that the discrepancies in density and temperature decrease with higher resolution, as expected.
        The deviations in the velocity profile in the subsonic region also decrease, although some differences are still present even at high resolution.

    

 §.§ Embedded star

        

        

        

        

        
        For the embedded star test (Table <ref>), we take the hot, massive, luminous star from Section <ref> and embed it in a constant density medium.
        The stellar wind will heat the gas and create a cavity around the star.
        This situation is quite common in embedded star clusters and it is what the heating wind mode is designed for.
        The initial gas is distributed along a grid[As mentioned in Section <ref>, this can introduce preferential directions and a glass or other relaxed system should be considered for most applications.] to ensure a constant density and a divergence-free random Gaussian velocity field following <cit.>.
        To ensure that the medium is stable, we use periodic boundary conditions and stop the simulation when the wind-blown bubble covers more than half the simulation box.
        For the heating wind mode, the outer radius for new wind particles (r_wind) is set manually to r_wind = 0.01 pc and the feedback efficiency is set to f_fb = 0.01 following <cit.>.

        
        In Figure <ref> we show the gas density when the bubble has just started to form (t=0.2 Myr) and when it has had time to grow (t=0.6 Myr).
        The stellar wind creates an approximately spherical bubble of lower density as the gas is swept up in a high density shell around it.
        To understand why the bubble is not perfectly spherical, we note that the finite number of gas particles cause small numerical fluctuations in the initial gas density.
        When we then introduce a small number of wind particles with higher energy than the surrounding gas, these small fluctuations grow into a larger asymmetry in the wind bubble.
        This growth of small initial asymmetries was observed in SPH simulations of supernovae explosions <cit.> where they found that if the injected energy is spread out over more particles, the asymmetric effects diminish.
        If the asymmetry in the wind bubbles would become a problem for specific simulations, then the wind energy could be spread out in a similar way.

        
        In Figure <ref> we have drawn a dashed line that shows the radius where the mean density is highest (r_ρmax, see Appendix <ref> for details).
        At the start of the simulation, this radius is undefined, because the gas has a constant density.
        As the bubble grows and gas is swept up in an approximately spherical shell, the radius of maximum density matches the shell radius, which is what we plot as a function of time in Figure <ref>.
        Note that r_ρmax is slightly larger than the shell radius because of the asymmetry of the wind bubble.
        We present this expansion for different values of M_SPH (different resolutions) and two different gas densities.
        Even when the resolution is very low (M_SPH = 1.0), the heating wind method still results in a dispersion of the gas cloud.
        The expansion is faster for lower gas density, which is in agreement with analytical solutions for the shell radius of an energy driven flow in a constant density medium <cit.>.
        However, the bubble expansion starts later for simulations with a lower resolution.
        This delay corresponds to the time it takes for the star to lose enough mass to create the first wind particle.
        For example, if M_SPH = 1.0 and = 1 /Myr, this delay is 1 Myr.
        In the bottom panel of Figure <ref>, we show the shell expansion starting at the moment of the first wind injection.
        We see that lower resolution results in a faster expansion, caused by the larger energy injected in a single SPH particle.
        The expansion profile approaches the analytical solution for high resolution (small M_SPH).
        We therefore advise that the choice of M_SPH be based on the stellar mass-loss rate and the delay and expansion profile that would be acceptable in the desired simulations.

    

 §.§ Supernova

        

        

        

        

        
        As discussed in Section <ref>, the heating wind mode can also be used to simulate the effect of a supernova on the surrounding gas.
        The supernova test (Table <ref>) is similar to the embedded star test (Section <ref>).
        We start the simulation by evolving the star using the stellar evolution code  <cit.> to a few timesteps (40 yr) before the star goes supernova (9.78 Myr).
        We place the star inside the uniform density gas medium and use the option to derive the stellar wind parameters from the output of a stellar evolution calculation (see Section <ref>).
        When this option is set,  detects the supernova and creates the particles with a combined mass of 12.95  and an energy of 10^51 erg.
        After the supernova feedback is generated we trace the resulting blast wave.
        Similar to the embedded star test, we get a sphere of high density material moving away from the star, however, due to the higher energy input, the radial velocity is higher.
        Using test simulations, we have found that a timestep of t_wind = 20 yr is required to avoid time-stepping artefacts with this high velocity <cit.>.

        
        In Figure <ref> we present the radial mean density profile for simulations with four different particle masses.
        We find that for the low resolution simulation (M_SPH = 1 ), the gas has moved away from the star, but the shape of the shockfront, where the density is highest, is only loosely defined.
        For the higher resolution simulations, we do see the shape of the main shockfront clearly, and the simulations agree on the radius of highest density at t = 1000 yr.
        However, the radius of the shockfront does not converge to the analytical solution.
        There is an increased density at r = 0 for the M_SPH = 0.01  simulation.
        This feature is present at some point for all high resolution simulations and is the result of a reverse density wave within the outgoing shockwave.
        These waves are an artefact of the hydrodynamical simulation method used, but they are not an accurate representation of the true physical process.
        They should not be confused with the reverse shock that takes place in real supernova remnants.
        While these density waves do subside after a few thousand years, the density inside the supernova bubble shortly after the explosion should not be considered correct.

        
        The time evolution of the expansion of gas from a supernova explosion is usually modelled in separate phases.
        The first phase is a free expansion, where the ejected gas moves at an approximately constant velocity, sweeping up the gas in the interstellar medium.
        After the mass of swept up gas is equal to the mass of the originally expelled gas, the expansion can be approximated as a pure adiabatic expansion, which is described by the self-similar Sedov-Taylor solution.
        Only this last phase can be simulated with the heating wind mode of , because the particles are given a high internal energy instead of an initial velocity.

        
        In Figure <ref> we present the time evolution of r_ρmax, which we calculate in the same way as for the embedded star test.
        We now compare it to the analytical solution for the two phases of a supernova blast wave, the free expansion and the Sedov-Taylor solution.
        We find that the simulations do approach the analytical solution and roughly follow the same shape, but even the highest resolution simulation expands faster than the analytical solution.
        We do not model the initial free expansion phase and shockwaves are in general difficult to simulate using SPH <cit.>.
        Differences with the analytical solution are therefore to be expected and this type of simulation should be interpreted with care.

        
        Given these caveats, both the use of SPH and the chosen approximations may not seem to be the ideal choice for simulating a supernova explosion in a gaseous medium.
        Indeed, depending on the goal of the simulations, other available methods could be more suitable, for example using a grid based simulation code <cit.> or including magnetic fields <cit.>.
        However, the method presented here has two main advantages:
        1) It is simple and scales well to very low SPH resolution, making it computationally faster than more detailed simulation techniques.
        2) The use of SPH combined with  allows easy gravitational coupling between the gas and the stars.
        We can therefore use this code to run large scale simulations of multiple supernova explosions in a gaseous medium also containing many dynamic stars.
        These advantages allow us to model a very turbulent stage in the evolution of embedded star clusters.

    

 §.§ Colliding wind triple

        

        

        

        
        The previous tests were for single stars and therefore the geometry of the outflow was not modified by the environment.
        In this test (Table <ref>), we simulate a triple star system where all three stars have a strong stellar wind.
        The system we simulate is loosely based[
            We are aware that the chosen values do not match the most up-to-date observations of the WR48 system.
            However, the goal of this investigation is to demonstrate the use of the  code, not to explain the observed system.]
        on WR48 (θ Muscae), which is a triple system <cit.> consisting of a WC5/WC6 + O6/O7V binary with a short period <cit.> and an O9.5/B0Iab star in a longer orbit <cit.> around that binary.
        For this simulation we have used similar numerical parameters to the fast wind test in Section <ref>.

        
        In Figure <ref> we show the gas density, temperature and velocity at the end of the simulation.
        Due to the large difference between the inner and outer orbital periods, the system appears similar to a normal colliding wind binary, which was assumed in previous models of WR48 <cit.>.
        However, the orbital motion of the inner binary creates a spiral pattern in the density and temperature distribution, which is very different from the wind from a single star.
        This spiral pattern creates high and low temperature regions in the shockfront where the wind from the inner binary collides with the wind from the third star.
        The observations of this shockfront can therefore be quite different from observations of a normal colliding wind binary.

        
        It is important to note that this simulation is just an example of what is possible with the  module and not an in-depth investigation into wind interactions in a triple star systems.
        For example, in the middle panels of Figure <ref> we can see that the temperature of the wind from the inner binary is extremely high (>10^8 K).
        These high temperatures are unrealistic because in reality the gas would cool, which is not taken into account in this simulation.
        When using this code for simulations that are to be compared with observations, gas cooling and a convergence test of the shockfront regions should be performed.



§ DISCUSSION AND CONCLUSION

        

        

        
        We have presented and tested the  module, which can be used to simulate stellar winds within the  framework by creating (and accelerating) SPH particles.
        The code includes three different modes: the simple mode, the accelerating mode  and the heating mode (Table <ref>).
        We have tested the code for single stars with fast and slow winds, as well as an embedded star with both wind and a supernova explosion.
        For both fast and slow winds, the simple and accelerating wind modes perform well, although subsonic winds must be simulated with the latter.
        For the embedded star, the heating wind mode creates a wind bubble, even at low resolution; with higher resolution the expansion profile approaches the analytical solution.
        After a supernova, the heating wind mode creates an expanding shell with velocities similar to the analytical solution if small enough timesteps are used.
        Finally we have shown an example of how this module can be used to tackle new problems, by simulating a colliding wind triple system.

        
        The  module can be used for many different simulations that involve stellar winds and several projects are already in progress.
        The simple wind mode has been used to simulate the accretion of gas from the winds of the S-stars onto the super-massive black hole Sgr A^* <cit.>.
        The accelerating wind mode is used to simulate the accretion of the wind from a red giant onto a close binary companion <cit.>.
        The heating wind mode is part of a larger simulation to investigate the evolution of the Arches cluster van der Helm et al..
        In Table <ref> we give an overview of the application of the different modes.
        The code is publicly available in the  framework.

        
        There are many other types of simulations involving stellar winds that could be done with the  framework and corresponding modes could be added to .
        It would be possible to add the mass and corresponding energy lost by stars to existing SPH particles.
        This would make it possible to run simulations of embedded stars with even lower resolution (higher SPH particle mass).
        However, this would result in unequal mass particles, which requires advanced treatment in the SPH codes.
        In the other extreme, since radiative transfer codes are available in , it would be possible to add a mode that solves the radiative hydrodynamics of the wind and this would make detailed stellar wind simulations possible.
        In fact, such coupled simulations have been performed with  already <cit.>.


    We thank N. Lützgendorf for testing and improving the simple wind mode, R. P. Kudritzki for his advice on the  scaling law and F. I. Pelupessi for providing his code for the heating wind mode.
    This work was supported by the Netherlands Research Council NWO.
    


    aa
    



§ EQUATIONS

        

        In this appendix, we calculate the analytical predictions for a stationary, spherically symmetric wind which are used in .
        For these calculations, we assume that the mass-loss rate () and the velocity as a function of radius (v(r)) are known and we define the acceleration
        
    a(r) = dv/dt = dv/drdr/dt = v(r) dv/dr.


    

 §.§ Radius as a function of time

        
        To calculate the outer radius of a new wind shell, we need to know the radius as a function of time (r(t)) where the wind starts at the stellar surface, so r(0) = R_*.
        Since v(r) is known, we can write
        
    v(r(t))     = d r(t)/dt, 
    
                dt     = dr/v(r),

        which is solved by,
        
    t = ∫_R_*^r(t)1/v(r)  dr.

        In general this equation has to be solved numerically[
            When solving the equations mentioned here numerically, we use the python library scipy (scipy.org).
            For integrals we use scipy.integrate.quad and for finding a root we use scipy.optimize.brentq.
            See docs.scipy.org for the details of these methods.
        ] for r(t), although for some velocity functions we can solve it analytically, for example if v(r) = then
        
    t     = 1/∫_R_*^r(t) dr = r(t) - R_*/, 
    
                r(t)     = R_* + t * .


    

 §.§ New particle radii

        
        When we create a new shell of particles, we want the density profile in the shell to match the density profile corresponding to the chosen velocity profile.
        To calculate that density profile, we first note that the mass-loss rate, Ṁ is related to the density and the velocity at any point of the wind via the equation of mass continuity,
        
    = 4π r^2 ρ(r) v(r),

        where ρ is the density of the wind.
        Because we assume that  and v(r) are known, we can rewrite this as
        
    ρ(r) = /4π r^2 v(r).


        To generate the positions of new particles, we start with a cube filled with particle positions with a uniform density.
        In our code, this can be a simple grid or randomly generated positions.
        From that cube, we remove all particles that are not inside the desired shell to get a shell of particles with uniform density.
        After that, we shift the particle positions in the radial direction to get the desired density profile.

        To find the new particle radius, we define the relative enclosed mass, x as
        
    x = ∫_R_*^r_pπ r^2 ρ(r) dr/∫_R_*^r(t)π r^2 ρ(r) dr,

        where r_p is the radius of the particle and R_* and r(t) are the inner and outer radius of the shell respectively.
        For the uniform density shell that was generated, this reduces to
        
    x_u = ∫_R_*^r_p r^2 dr/∫_R_*^r(t) r^2 dr = r_p^3 - R_*^3/r(t)^3 - R_*^3.

        For the desired density profile based on a given velocity profile, we rewrite equation <ref> in terms of v using equation <ref>
        
    x_v = ∫_R_*^r_p/v(r) dr/∫_R_*^r(t)/v(r) dr = ∫_R_*^r_p1/v(r) dr/∫_R_*^r(t)1/v(r) dr.

        We then set x_u = x_v where x_u is calculated with the old particle radius (of the generated uniform density shell).
        The last step is to solve equation <ref> to get the new particle radius r_p.
        In general this equation has to be solved numerically, although for some velocity functions we can solve it analytically, for example if v(r) = then
        
    x     = ∫_R_*^r_p1/ dr/∫_R_*^r(t)1/ dr = ∫_R_*^r_p dr/∫_R_*^r(t) dr = r_p - R_*/r(t) - R_*, 
    
                r_p     = R_* + x (r(t) - R_*).


    

 §.§ Gas pressure

        
        To calculate the expected acceleration, a_P(r) caused by the gradient of the gas pressure, P(r) we assume a polytropic equation of state,
        
    P = K ρ(r)^γ,

        where K is the polytropic constant and γ = 5/3 is the adiabatic index for a monoatomic ideal gas.
        Because K is constant we can calculate it at the surface of the star and use that value for the entire wind.
        To calculate P(R_*) we use
        
    P(r) = (γ - 1)ρ(r) u,

        where u is the internal energy of the gas particles defined by
        
    u = 3/2k_B T_*/μ,

        where k_B is the Boltzmann constant, T_* is the temperature at the photosphere of the star and μ is the mean molecular weight of the gas particles.
        Combining equations <ref> and <ref> we get
        
    K = u(γ-1)ρ(R_*)^1-γ.


        The acceleration caused by the gradient of the gas pressure is
        
    a_P(r) = - 1/ρ(r)∂ P(r)/∂ r,

        which we can rewrite using equations <ref> and <ref>
        
    a_P(r)     = - K/ρ(r)∂ρ^γ/∂ r
        = - K/ρ(r)γρ(r)^γ - 1∂/∂ r/4π r^2 v(r)
        = K γρ(r)^γ - 1( 2/r + 1/v(r)dv(r)/dr ).


    

 §.§ Density as a function of radius

        
        In Sections <ref> and <ref> we calculate the density as a function of radius.
        For each radius r, we take six points in six directions (+r and -r along each axis x, y and z) and calculate the SPH density at those points.
        Note that there does not need to be an SPH particle at that point to calculate the density.
        We then take the mean of these 6 densities to be the density at that radius.
        To calculate the radius with maximum density r_ρmax, we calculate this for a grid of radii and select the radius with the largest density.

