{'fraction_non_alphanumeric': 0.031007039535353907, 'fraction_numerical': 0.01290275148418402, 'mean_word_length': 4.432651023464802, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulating stellar winds in AMUSE

We present stellar_wind.py, a module that provides multiple methods of simulating stellar winds using smoothed particle hydrodynamics codes (SPH) within the astrophysical multipurpose software environment (AMUSE) framework. With the simple wind mode, we create SPH wind particles in a spherically symmetric shell. 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. 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.


Introduction
Stars lose mass through stellar winds during various stages of their evolution (e.g. Meyer-Vernet 2007;Owocki 2013;Puls et al. 2015). These winds can affect the gas near the star, creating lower density bubbles (Castor et al. 1975) and regulating star formation (Oey & Clarke 2009). If a binary companion is present, accretion of the stellar wind material can also affect the evolution of that companion (Boffin 2015).
The two most important parameters of the stellar wind are the mass-loss rate,Ṁ, and the terminal wind velocity, v ∞ , which determine the effect of the wind on the environment. Based on these parameters, stellar winds can be broadly divided into three categories (Owocki 2013): 1) Cool main-sequence stars like the Sun have winds with very low mass-loss rates (Ṁ∼10 −14 M /yr) that are thermally or gas pressure driven. 2) Cool giants and super giants have slow (v ∞ ∼5 -30 km/s) high mass-loss rate (Ṁ∼10 −7 -10 −5 M /yr) winds driven mainly by radiation pressure on dust (Höfner 2015). 3) Hot luminous stars have fast (v ∞ ∼200 -3000 km/s) high mass-loss rate (Ṁ∼10 −7 -10 −5 M /yr) line driven winds (Puls et al. 2009). 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 (e.g. Owocki et al. 1988;Blondin et al. 1990; Kudritzki & Puls 2000;Boffin 2015). To investigate the net effect of the stellar wind on the environment, it is often sufficient to simulate the stellar wind using only hydrodynamics (Theuns & Jorissen 1993;Cuadra et al. 2006;Mohamed et al. 2012). 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 (e.g. Agertz et al. 2013;Muratov et al. 2015).
For all these simulations, the astrophysical multipurpose software environment (amuse 1 ; Portegies Zwart et al. 2013;Pelupessy et al. 2013;van Elteren et al. 2014;Portegies Zwart & McMillan 2018;Portegies Zwart et al. 2018) 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 amuse (e.g. , a consistent and properly tested module was still missing. The stellar wind.py code presented in this paper can be combined with the SPH (Smoothed particle hydrodynamics), Nbody, stellar evolution and (with some additional work) radiative transfer codes that are already available. We describe the stellar wind.py code and explain the different modes in which it can be used (Section 2). In Section 3 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 1 amusecode.org 1 arXiv:1903.11098v1 [astro-ph.SR] 26 Mar 2019 Edwin van der Helm et al.: Simulating stellar winds in AMUSE Section 4 with an exposition of some ongoing research projects using this code and ideas for further use.

Methods
The goal of stellar wind.py is to create gas particles that represent the stellar wind from one or more stars. The code requires a number of stars, represented by amuse particles 2 , 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 amuse which simulates the hydrodynamics of the wind.
Creating the SPH particles is only one step in the simulations for which stellar wind.py is used. Following the goal of the amuse framework, the other parts of the simulations are handled by specialized interchangeable codes. For the hydrodynamics, SPH codes such as fi (Pelupessy 2005) and gadget2 (Springel 2005) 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, bridge (Fujii et al. 2007) is available. The stellar properties on which the wind is based will generally be calculated using a stellar evolution code. Both parametrized (e.g. Hurley et al. 2000) and Henyey type (e.g. Paxton et al. 2011) stellar evolution codes are available in amuse. Any or all of these codes can be combined with stellar wind.py to set up a wide variety of simulations (see Section 3). For more information about the codes available within amuse and examples of how to couple them, we refer the reader to Portegies Zwart & McMillan (2018).

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 amuse.
The stellar wind.py 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 amuse. However, the terminal wind velocity, v ∞ , is not calculated by any code currently in amuse. Determining v ∞ 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 stellar wind.py a formula that has been fitted to observations of hot stars (Kudritzki & Puls 2000), and which, they claim, is valid for these stars within 20%. The 2 A particle set is the fundamental data structure in amuse. 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. terminal velocity of the wind is given by: T * ≤ 10 000 K, 1.4 10 000 K < T * < 21 000 K, 2.65 T * ≥ 21 000 K, where v phesc is the photospheric escape velocity (similar to the escape velocity v esc 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, v ∞ ≈ v phesc and this formula is still applicable (Kudritzki, private communication).

Simple wind
Within stellar wind.py, there are currently three wind modes available. The simplest mode creates a spherical shell of particles around the star with radial velocity, v(r) = v ∞ 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 (e.g. Mohamed et al. 2012) and it serves as a starting point for the two other modes described in Sections 2.3 and 2.4. 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) = v ∞ 2 + v esc (r) 2 where v esc (r) = √ 2GM * /r is the local escape velocity at the initial particle radius, r. We calculate this new velocity for each particle because v esc (r) can vary within the thin shell in which we create the particles.
We set the outer radius (r max ) 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 A.1). We scale the particle positions within the shell to follow the density profile matching the velocity profile as described in Appendix A.2.

SPH and initial distributions
SPH is a method to solve the dynamics of a fluid by approximating it with a set of discrete particles (Monaghan 1992). 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 3 and a particle smoothing length (h) are used. In all modes of stellar wind.py, 1. An example of the initial positions of newly created particles using a random (top) and grid (bottom) distribution. A shell of particles was created between 1 and 3 stellar radii (R * ) and the x and y positions of a thin slice (|z| < 0.05 R * ) are shown. The positions are scaled to match the given density profile (see Appendix A). Note that this is merely an illustration of the difference between random and grid initial distributions. In real simulations, the shells would generally be much thinner.
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 (e.g. see Pelupessy 2005).
Creating an initial distribution of SPH particle positions is not trivial (e.g. Diehl et al. 2015). 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 direc-tions (for example a 'glass' initial condition, White 1996;Wang & White 2007). 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 stellar wind.py 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 1. 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 1.
There are many more ways to create initial particle distributions. A good overview of the different methods and their advantages can be found in Diehl et al. (2015). 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 3.

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 < 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 bridge (Fujii et al. 2007). Originally, bridge was designed to couple multiple gravitational codes. In this method, each code is evolved separately for a short, predefined timestep 4 . The mutual gravitational effect is included by bridge using a kick-drift-kick scheme (see e.g. Portegies Zwart & McMillan 2018). 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 stellar wind.py, we use bridge by including an artificial potential near the star, and then use the same kickdrift-kick scheme to ensure a smooth acceleration of the wind particles.  Table 1. For the beta law we show two curves that are good fits for hot massive stars (β = 0.8) and cool giants (β = 2.0) (Lamers & Cassinelli 1999). To illustrate the shape of the acceleration curves, we have chosen v 0 = 0.2 v ∞ , r acc start = 2 R * (for the delayed rsquared function), and r mid = 3 R * and s = α = 10 (for the logistic and agb function).
In Figure 2 and Table 1 we present the acceleration functions (sometimes referred to as acceleration laws) currently implemented in stellar wind.py. For a given velocity (or acceleration) profile, all required quantities are calculated following Table 1. An overview of the acceleration functions currently available in stellar wind.py. Either the acceleration (a) or the velocity (v) is given depending on which is simpler. The corresponding v or a function can be derived using a(r) = v(r) dv dr and the known boundary conditions. Some functions allow user defined parameters to affect the functions (e.g. r acc start , r mid , β, etc.) The first three functions are rough approximations to the last three functions. Their advantage is that they are computationally faster.

Name
Equation Use Hot/cool stars the equations in Appendix A. 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 Nowotny et al. (2010). These winds exhibit a specific acceleration zone, the location of which can be chosen with the parameters r mid 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 r mid = 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 Lamers & Cassinelli (1999) and Maciel (2014). 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 r acc start (with a default value of 2) sets the lower boundary of the acceleration zone. The initial velocity v 0 , 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 v 0 result in high densities which lead to slow simulations, so in many cases a higher value of v 0 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 A.3) 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 stellar wind.py 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 2.2.1. 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 A.9 for each particle. For most acceleration functions, this equation has to be solved numerically, which can severely slow down the simulation (see Section 3.1). We therefore include the option to define a critical timestep, t c . When new wind particles are created, δt, which determines the thickness of the shell of new particles, is compared to t c . If δt < t c 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 A.9 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 2).

Heating wind
The third wind mode is based on the method used in  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 (E mech ) of the accumulated wind, defined as, where L mech 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 stellar wind.py during the simulation. The internal energy of the new particles is set to 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 3.3 and 3.4. As discussed in , 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

Application
To ensure that stellar wind.py 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 3.1 and 3.2 are simulations of processes that happen close to the star. Therefore only the simple and accelerating wind modes are applicable. The tests in Section 3.3 and 3.4 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 3.5, where we couple stellar dynamics, hydrodynamics and stellar winds, shows the power of stellar wind.py within amuse 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 stellar wind.py module is designed to couple different parts of a simulation in amuse. 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 3 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 fi with an adiabatic equation of state and artificial viscosity parameters α = 0.5 and β = 1.0 (following Lombardi et al. 1999). 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 bridge, which is also used in the accelerating wind mode (Section 2.3). The bridge code requires an additional code for calculating the gravitational force. When we simulate only a single star that does not evolve dynamically, we use fastkick 5 . When we simulate multiple stars that evolve dynamically (Section 3.5), we use the Nbody code huayno . For each simulation, we also need to define a number of integration timescales, such as the bridge timestep (t br ), the (maximum) internal timestep (t N−body ) 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 ≥ t N−body 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 = 2t br = 4t N−body . 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 2. 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 stellar wind.py. The initial wind velocity, v 0 = 100 km/s, is based on numerical considerations. As we shall see in Section 3.2, slow (subsonic 6 ) winds are more complex to simulate and for many applications using a higher v 0 is sufficient. In Figure 4 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 fol-5 The fastkick code, developed by N. de Vries, is an unpublished gpu-enabled code in amuse 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 bridge. 6 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. Fig. 4. The analytical and simulated velocity (v), density (ρ) and temperature (T ) as a function of radius (r) for the fast wind test. The results are from simulations using the simple wind mode (top) and the accelerating wind mode (bottom). The stellar and wind parameters can be found it Table 2. To calculate the analytical temperature profile, we assume adiabatic expansion. low 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 5 that for a higher resolution (smaller M gas ), 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 6 we present the time spent in different parts of the simulation code (t code ) divided by the total cpu  Figure 4, but only for a simulation using the accelerating wind mode with the logistic acceleration function. We have varied the resolution by changing the SPH particle mass (M SPH ) and through that the number of particles in the simulation. We have added the smoothing length, h, as a function of radius for each simulation, which is a measure of the local spatial resolution.
(or wall-clock) time (t tot ) as a function of resolution 7 . In the top 7 These simulations where all performed on the same desktop computer using a 4-core Intel Xeon E5507 CPU. The part of the simulation time used for the stellar wind.py code while creating new particles (circles) and accelerating them (stars) compared to the time used by the SPH code (diagonal lines). The three panels show results for simulations with the simple wind mode (top), the accelerating wind mode with the logistic acceleration function without a critical timestep (middle) and the accelerating wind mode with a critical time step (bottom). Marks along each line denote separate simulation runs. At the top axis we give an estimate of the number of SPH particles (N) actually used in the simulation with the corresponding particle mass (M SPH ). The remaining simulation time (white space) was mostly spent on unoptimized administrative tasks like saving snapshots and removing escaping particles.
panel we see that when using the simple wind mode, the time spent in the stellar wind.py 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 A.9 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 (t c , see Section 2.3). In the Table 3. Stellar and wind parameters used in the slow wind test. Derived parameters are indicated with an arrow (→). Since the smoothing length is highly variable with extreme outliers, we include the median value of all gas particles shown in Figure 8. bottom panel we see that by using this approximation, the time spent in stellar wind.py 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 3. 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 7 we present the velocity profiles of simulations using the simple wind mode and varying v ∞ 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 8 we present the results of simulations using the accelerating wind mode. In these simulations we have used the staging shell option (Section 2.3) that enforces the correct particle velocity for all SPH particles with radius r < 1.1R * . 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 = 1500R ).
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 (Nowotny et al. 2010). In these profiles, the velocities in the subsonic region oscillate due to stellar pulsations and do not follow the simple  Figure 4 but for the slow wind simulations using the simple wind mode without gravity where we vary the wind velocity. We have added the expected local sound speed (dotted line) for comparison. 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 4 are also present for the slow wind test in Figure 8. In Figure 9 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 4), we take the hot, massive, luminous star from Section 3.1 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   field following Bonnell et al. (2003). 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 .
In Figure 10 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 (Rimoldi et al. 2016) 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 10 we have drawn a dashed line that shows the radius where the mean density is highest (r ρmax , see Appendix A.4 for details). At the start of the simulation, this radius is unde-  Fig. 11. The radius with the highest mean gas density (r ρmax ) as a function of time, t for the embedded star simulations. Different colours correspond to different resolutions resulting from different SPH particle masses (M SPH ). In the top panel, we show r ρmax for all simulations as a function of t. Lines with open circles correspond to simulations where the gas density, ρ gas = 10 M /pc 3 and lines with filled circles to simulations with a gas density, ρ gas = 100 M /pc 3 . In the bottom panel, we only show simulations with ρ gas = 100 M /pc 3 and subtractṀ/M SPH (the time of the first SPH particle creation) from t. The black solid line shows the analytical solution for the shell radius of an energy driven flow in a constant density medium (Dyson 1984).
fined, 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 11. 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 M ), the heating wind method still results in a dispersion of the gas cloud. The expansion is faster for lower gas density, which is   Fig. 12. The mean gas density as a function of radius at time t = 1000 yr for the supernova test at four resolutions (dashed lines). The analytical solution (solid black line), is the Sedov-Taylor solution for a self-similar blast wave in a uniform medium (Taylor 1950;Sedov 1959).
in agreement with analytical solutions for the shell radius of an energy driven flow in a constant density medium (Dyson 1984). 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 M andṀ = 1 M /Myr, this delay is 1 Myr. In the bottom panel of Figure 11, 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 2.4, the heating wind mode can also be used to simulate the effect of a supernova on the surrounding gas. The supernova test (Table 5) is similar to the embedded star test (Section 3.3). We start the simulation by evolving the star using the stellar evolution code SeBa (Portegies Zwart & Verbunt  Fig. 13. The radius with the highest mean gas density (r ρmax ) as a function of time for the supernova test. We include the analytical solution, which is free expansion followed by the Sedov-Taylor solution.
2012) 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 2.1). When this option is set, stellar wind.py detects the supernova and creates the particles with a combined mass of 12.95 M 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 timestepping artefacts with this high velocity (also see: . In Figure 12 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 M ), 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 M 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 L * (L ) 2 · 10 5 1.4 · 10 5 5.5 · 10 4 temperature T * (K) 7.1 · 10 4 4.5 · 10 4 3.9 · 10 4 orbital period p 19.14 days 130 yr eccentricity e 0 0 inclination i 0 • wind timestep t wind 0.2 days end time t end 190 days particle mass M SPH 10 −11 M the originally expelled gas, the expansion can be approximated as a pure adiabatic expansion, which is described by the selfsimilar Sedov-Taylor solution. Only this last phase can be simulated with the heating wind mode of stellar wind.py, because the particles are given a high internal energy instead of an initial velocity.
In Figure 13 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 (e.g. Hubber et al. 2013). 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 (e.g. Rogers & Pittard 2013) or including magnetic fields (e.g. Körtgen et al. 2016). 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 bridge 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 6), we simulate a triple star system where all three stars have a strong stellar wind. The system we simulate is loosely based 9 on WR48 (θ Muscae), which is a triple system (Sugawara et al. 2008) consisting of a WC5/WC6 + O6/O7V binary with a short period (∼19 days, Hill et al. 2002) and an O9.5/B0Iab star in a longer orbit (> 130 yrs, Dougherty & Williams 2000) around that binary. For this simulation we have used similar numerical parameters to the fast wind test in Section 3.1.
In Figure 14 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 (Hill et al. 2002). 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 stellar wind.py module and not an in-depth investigation into wind interactions in a triple star systems. For example, in the middle panels of Figure 14 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 stellar wind.py module, which can be used to simulate stellar winds within the amuse framework by creating (and accelerating) SPH particles. The code includes three different modes: the simple mode, the accelerating mode and the heating mode (Table 7). 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 stellar wind.py 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 * (Lützgendorf et al. 2016). The accelerating wind mode is used to simulate the accretion of the wind from a red giant onto a close binary companion (Saladino et al. 2018). The heating wind mode is part of a larger simulation to investigate the evolution of the Arches cluster (van der Helm et al. in prep.). In Table 7 we give an overview of the application of the different modes. The code is publicly available in the amuse framework.
investigation is to demonstrate the use of the stellar wind.py code, not to explain the observed system. The gas density (top), temperature (center) and velocity (bottom) in the orbital plane of the colliding wind triple simulation for the inner (left) and outer (right) binary. The sizes of the stars (yellow circles) in the plots on the right hand panels were multiplied by 10 to make them visible. In the left hand panels, the Wolf-Rayet star (star 1, see Table 6) can be seen on the right and star 2 on the left. In the right hand panels, the short period binary (star 1 and 2), can be seen on the left and the O9.5 supergiant (star 3) on the right. In the bottom plots, the arrows indicate the wind direction and larger arrows correspond to higher wind velocities, however, the colors provide a more precise indication of the velocities. Note that the two density plots have separate color bars, while the temperature and velocity plots each share a single color bar. There are many other types of simulations involving stellar winds that could be done with the amuse framework and corresponding modes could be added to stellar wind.py. 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 amuse, 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 amuse already (Wall et al. 2017, N. Clementel, private communication In general this equation has to be solved numerically 10 for r(t), although for some velocity functions we can solve it analytically, for example if v(r) = v ∞ then A.2. 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), (A.5) 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) . (A.6) 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 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 For the desired density profile based on a given velocity profile, we rewrite equation A.7 in terms of v using equation A.6 x v = 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 A.9 to get the new particle radius r p . In general 10 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. this equation has to be solved numerically, although for some velocity functions we can solve it analytically, for example if v(r) = v ∞ then r p − R * r(t) − R * , r p = R * + x(r(t) − R * ). (A.10)

A.3. 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) γ , (A.11) 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, (A.12) where u is the internal energy of the gas particles defined by 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 A.11 and A.12 we get K = u(γ − 1)ρ(R * ) 1−γ . (A.14) The acceleration caused by the gradient of the gas pressure is a P (r) = − 1 ρ(r) ∂P(r) ∂r , (A.15) which we can rewrite using equations A.11 and A.6 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 . (A.16)

A.4. Density as a function of radius
In Sections 3.3 and 3.4 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.