Issue 
A&A
Volume 625, May 2019



Article Number  A85  
Number of page(s)  15  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201732020  
Published online  17 May 2019 
Simulating stellar winds in AMUSE
^{1}
Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
email: edwin.vanderhelm@gmail.com
^{2}
Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands
Received:
30
September
2017
Accepted:
24
March
2019
Aims. 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.
Methods. 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.
Results. 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.
Key words: stars: winds, outflows / methods: numerical / hydrodynamics
© ESO 2019
1. Introduction
Stars lose mass through stellar winds during various stages of their evolution (e.g. MeyerVernet 2007; Owocki et al. 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 2014).
The two most important parameters of the stellar wind are the massloss 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 (Owocki et al. 2013): (1) Cool mainsequence stars like the Sun have winds with very low massloss rates (Ṁ∼10^{−14} M_{⊙} yr^{−1}) that are thermally or gas pressure driven. (2) Cool giants and super giants have slow (υ_{∞}∼5−30 km s^{−1}) high massloss rate (Ṁ∼10^{−7}−10^{−5} M_{⊙} yr^{−1}) winds driven mainly by radiation pressure on dust (Höfner 2015). (3) Hot luminous stars have fast (υ_{∞}∼200−3000 km s^{−1}) high massloss rate (Ṁ∼10^{−7}−10^{−5} M_{⊙} yr^{−1}) 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 2014). 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 subgrid 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, 2018; Pelupessy et al. 2013; van Elteren et al. 2014; Portegies Zwart & McMillan 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. Pelupessy & Portegies Zwart 2012), a consistent and properly tested module was still missing.
The STELLAR_WIND.PY code presented in this paper can be combined with the smoothed particle hydrodynamics (SPH), 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 (Sect. 2). In Sect. 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 Sect. 4 with an exposition of some ongoing research projects using this code and ideas for further use.
2. 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 massloss 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 Nbody 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 Sect. 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).
2.1. Calculating stellar wind properties
To simulate the stellar winds, the stellar parameters (mass, radius, temperature and position) and wind parameters (massloss 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 userfriendly routines to derive some of the stellar parameters such as stellar mass, massloss rate, stellar radius and effective temperature from one of the stellar evolution codes within AMUSE. However, the terminal wind velocity, υ_{∞}, is not calculated by any code currently in AMUSE. 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 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 terminal velocity of the wind is given by:
$${\upsilon}_{\infty}=C({T}_{\ast}){\upsilon}_{\mathrm{phesc}},$$
where,
$$\begin{array}{cc}\hfill & C({T}_{\ast})=\{\begin{array}{cc}1\hfill & {T}_{\ast}\le 10\phantom{\rule{0.166667em}{0ex}}000\phantom{\rule{3.33333pt}{0ex}}\mathrm{K},\hfill \\ 1.4\hfill & 10\phantom{\rule{0.166667em}{0ex}}000\phantom{\rule{3.33333pt}{0ex}}\mathrm{K}<{T}_{\ast}<21\phantom{\rule{0.166667em}{0ex}}000\phantom{\rule{3.33333pt}{0ex}}\mathrm{K},\hfill \\ 2.65\hfill & {T}_{\ast}\ge 21\phantom{\rule{0.166667em}{0ex}}000\phantom{\rule{3.33333pt}{0ex}}\mathrm{K},\hfill \end{array}\hfill \\ \hfill & {\upsilon}_{\mathrm{phesc}}=\sqrt{2{g}_{\ast}{R}_{\ast}(1\mathrm{\Gamma})},\hfill \\ \hfill & {g}_{\ast}=\frac{G{M}_{\ast}}{{R}_{\ast}^{2}},\hfill \\ \hfill & \mathrm{\Gamma}=7.66\times {10}^{5}{\sigma}_{\mathrm{e}}\frac{{L}_{\ast}/{L}_{\odot}}{{M}_{\ast}/{M}_{\odot}},\hfill \\ \hfill & {\sigma}_{\mathrm{e}}=0.398\frac{1+{I}_{\mathrm{He}}\mathrm{Y}}{1+4\mathrm{Y}},\hfill \end{array}$$(1)
where υ_{phesc} is the photospheric escape velocity (similar to the escape velocity υ_{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, υ_{∞} ≈ υ_{phesc} and this formula is still applicable (Kudritzki, priv. comm.).
2.2. 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, υ(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 (e.g. Mohamed et al. 2012) and it serves as a starting point for the two other modes described in Sects. 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 $\upsilon (r)=\sqrt{{{\upsilon}_{\infty}}^{2}+{\upsilon}_{\mathrm{esc}}{(r)}^{2}}$ where ${\upsilon}_{\mathrm{esc}}(r)=\sqrt{2G{M}_{\ast}/r}$ is the local escape velocity at the initial particle radius, r. We calculate this new velocity for each particle because υ_{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.
2.2.1. 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, 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 directions (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 Fig. 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 Fig. 1.
Fig. 1. 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. 
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 Sect. 3.
2.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, υ < υ_{∞}. 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 kickdriftkick scheme (see e.g. Portegies Zwart & McMillan 2018). This method can also be used to gravitationally couple a pure Nbody 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 kickdriftkick scheme to ensure a smooth acceleration of the wind particles.
In Fig. 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 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 timeaveraged 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 betalaw, 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 υ_{0}, which is used in all functions except for constant_velocity, is of the order of a few km s^{−1} due to microturbulence in the stellar atmosphere where the material is launched. We note that low values of υ_{0} result in high densities which lead to slow simulations, so in many cases a higher value of υ_{0} can be used as an approximation. In addition to these predefined functions, new user defined velocity functions can easily be incorporated.
Fig. 2. Radial velocity profiles for the wind acceleration functions currently available in STELLAR_WIND.PY. The formulae for these functions can be found in 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 υ_{0} = 0.2 υ_{∞}, r_{acc_start} = 2 R_{*} (for the delayed_rsquared function), and r_{mid} = 3 R_{*} and s = α = 10 (for the logistic and agb function). 
Overview of the acceleration functions currently available in STELLAR_WIND.PY.
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 Sect. 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 Eq. (A.9) for each particle. For most acceleration functions, this equation has to be solved numerically, which can severely slow down the simulation (see Sect. 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 Eq. (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 Fig. 2).
2.4. Heating wind
The third wind mode is based on the method used in Pelupessy & Portegies Zwart (2012) 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,
$$\begin{array}{cc}& {E}_{\mathrm{mech}}={\displaystyle {\int}_{{t}_{0}}^{{t}_{1}}{L}_{\mathrm{mech}}(t)\phantom{\rule{0.166667em}{0ex}}\mathrm{d}t},\hfill \\ \hfill & {L}_{\mathrm{mech}}(t)=\frac{1}{2}\dot{M}(t){\upsilon}_{\infty}{(t)}^{2},\hfill \end{array}$$(2)
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
$$\begin{array}{c}\hfill u={f}_{\mathrm{fb}}\frac{{E}_{\mathrm{mech}}}{\mathrm{\Delta}{M}_{\ast}},\end{array}$$(3)
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 Sects. 3.3 and 3.4.
As discussed in Pelupessy & Portegies Zwart (2012), 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 timestepping artefacts (e.g. Pelupessy & Portegies Zwart 2012). One way to prevent this is to use a very small timestep (preliminary tests suggest ∼10 yr) shortly after a supernova.
3. 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 Sects. 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 Sects. 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 Sect. 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 Fig. 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). Selfgravity 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 (Sect. 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 (Sect. 3.5), we use the Nbody code HUAYNO (Pelupessy et al. 2012). 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_{Nbody}) of the SPH and Nbody 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 leapfrog algorithm to work, we should set t_{br} ≥ t_{Nbody} 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_{Nbody}. For the hydro simulation we also need to set the SPH particle mass (M_{SPH}).
Fig. 3. Diagram of the AMUSE codes (boxes) and their relations (labelled arrows). Dotted lines indicate optional codes that can be added depending on the astrophysical application. 
3.1. 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, υ_{0} = 100 km s^{−1}, is based on numerical considerations. As we shall see in Sect. 3.2, slow (subsonic^{6}) winds are more complex to simulate and for many applications using a higher υ_{0} is sufficient.
Stellar and wind parameters used in the fast wind test.
In Fig. 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 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 Fig. 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.
Fig. 4. Analytical and simulated velocity (υ), 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. 
Fig. 5. Same as Fig. 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. 
In addition to being accurate, it is also important that a simulation code is fast. In Fig. 6 we present the time spent in different parts of the simulation code (t_{code}) divided by the total cpu (or wallclock) time (t_{tot}) as a function of resolution^{7}. In the top 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 Eq. (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 Sect. 2.3). In the bottom panel we see that by using this approximation, the time spent in STELLAR_WIND.PY reduces to <1% for >10^{4} particles.
Fig. 6. 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. 
3.2. 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.
Stellar and wind parameters used in the slow wind test.
In Fig. 7 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.
Fig. 7. Same as Fig. 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. 
In Fig. 8 we present the results of simulations using the accelerating wind mode. In these simulations we have used the staging shell option (Sect. 2.3) 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 R_{⊙}).
Fig. 8. Same as Fig. 7 but for the accelerating wind mode with four different acceleration functions. 
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 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 Fig. 4 are also present for the slow wind test in Fig. 8. In Fig. 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.
3.3. Embedded star
For the embedded star test (Table 4), we take the hot, massive, luminous star from Sect. 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 star clusters and it is what the heating wind mode is designed for. The initial gas is distributed along a grid^{8} to ensure a constant density and a divergencefree random Gaussian velocity field following Bonnell et al. (2003). To ensure that the medium is stable, we use periodic boundary conditions and stop the simulation when the windblown 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 Pelupessy & Portegies Zwart (2012).
Parameters used in the embedded star test.
In Fig. 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.
Fig. 10. Gas density in the x − y plane after 0.2 Myr (top) and after 0.6 Myr (bottom) for the embedded star simulation with M_{SPH} = 0.1 M_{⊙} and ρ_{gas} = 100 M_{⊙} pc^{−3}. The embedded star is positioned at the origin (yellow dot) and the red dashed circle shows the radius with the largest mean density. 
In Fig. 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 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 Fig. 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 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^{−1}, this delay is 1 Myr. In the bottom panel of Fig. 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 massloss rate and the delay and expansion profile that would be acceptable in the desired simulations.
Fig. 11. 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}). Top panel: 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}. Bottom panel: 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). 
3.4. Supernova
As discussed in Sect. 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 (Sect. 3.3). We start the simulation by evolving the star using the stellar evolution code SEBA (Portegies Zwart & Verbunt 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 Sect. 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: Pelupessy & Portegies Zwart 2012).
Parameters used in the supernova test.
In Fig. 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.
Fig. 12. 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 selfsimilar blast wave in a uniform medium (Taylor 1950; Sedov 1959). 
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 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 Fig. 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.
Fig. 13. 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. 
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.
3.5. 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 yr, Dougherty & Williams 2000) around that binary. For this simulation we have used similar numerical parameters to the fast wind test in Sect. 3.1.
Stellar, wind and orbital parameters of the colliding wind triple simulation.
In Fig. 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.
Fig. 14. 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. Left hand panels: WolfRayet star (star 1, see Table 6) can be seen on the right and star 2 on the left. 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. 
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 indepth investigation into wind interactions in a triple star systems. For example, in the middle panels of Fig. 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.
4. 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.
Overview of the modes in STELLAR_WIND.PY and their suggested application domains.
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 Sstars onto the supermassive 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.
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, Clementel, priv. comm.).
As mentioned in Sect. 2.2.1, this can introduce preferential directions and a glass or other relaxed system should be considered for most applications.
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.
Acknowledgments
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.
References
 Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Blondin, J. M., Kallman, T. R., Fryxell, B. A., & Taam, R. E. 1990, ApJ, 356, 591 [NASA ADS] [CrossRef] [Google Scholar]
 Boffin, H. M. J. 2014, ASSL, 413, 153 [Google Scholar]
 Bonnell, I. A., Bate, M. R., & Vine, S. G. 2003, MNRAS, 343, 413 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J., McCray, R., & Weaver, R. 1975, ApJ, 200, L107 [NASA ADS] [CrossRef] [Google Scholar]
 Cuadra, J., Nayakshin, S., Springel, V., & Di Matteo, T. 2006, MNRAS, 366, 358 [NASA ADS] [CrossRef] [Google Scholar]
 Diehl, S., Rockefeller, G., Fryer, C. L., Riethmiller, D., & Statler, T. S. 2015, PASA, 32, e048 [NASA ADS] [CrossRef] [Google Scholar]
 Dougherty, S. M., & Williams, P. M. 2000, MNRAS, 319, 1005 [NASA ADS] [CrossRef] [Google Scholar]
 Dyson, J. E. 1984, Ap&SS, 106, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Fujii, M., Iwasawa, M., Funato, Y., & Makino, J. 2007, PASJ, 59, 1095 [NASA ADS] [Google Scholar]
 Hill, G. M., Moffat, A. F. J., & StLouis, N. 2002, MNRAS, 335, 1069 [NASA ADS] [CrossRef] [Google Scholar]
 Höfner, S. 2015, Why Galaxies Care about AGB Stars III: A Closer Look in Space and Time (Vienna, Austria), 497, 333 [NASA ADS] [Google Scholar]
 Hubber, D. A., Falle, S. A. E. G., & Goodwin, S. P. 2013, MNRAS, 432, 711 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Körtgen, B., Seifried, D., Banerjee, R., VázquezSemadeni, E., & ZamoraAvilés, M. 2016, MNRAS, 459, 3460 [NASA ADS] [CrossRef] [Google Scholar]
 Kudritzki, R.P., & Puls, J. 2000, ARA&A, 38, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999, Introduction to Stellar Winds (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Lombardi, J. C., Sills, A., Rasio, F. A., & Shapiro, S. L. 1999, J. Comput. Phys., 152, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Lützgendorf, N., van der Helm, E., Pelupessy, F. I., & Portegies Zwart, S. 2016, MNRAS, 456, 3645 [NASA ADS] [CrossRef] [Google Scholar]
 Maciel, W. J. 2014, Hydrodynamics and Stellar Winds: An Introduction (Cham: Springer) [CrossRef] [Google Scholar]
 MeyerVernet, N. 2007, Basics of the Solar Wind (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Mohamed, S., Mackey, J., & Langer, N. 2012, A&A, 541, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Monaghan, J. J. 1992, ARA&A, 30, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Muratov, A. L., Kereš, D., FaucherGiguère, C.A., et al. 2015, MNRAS, 454, 2691 [NASA ADS] [CrossRef] [Google Scholar]
 Nowotny, W., Höfner, S., & Aringer, B. 2010, A&A, 514, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oey, M. S., & Clarke, C. J. 2009, Massive Stars: Feedback Effects in the Local Universe (Cambridge: Cambridge University Press), 74 [Google Scholar]
 Owocki, S. 2013, in Stellar Winds, eds. T. D. Oswalt, & M. A. Barstow, 735 [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Pelupessy, F. I. 2005, PhD Thesis, Leiden Observatory, Leiden University, Leiden, The Netherlands [Google Scholar]
 Pelupessy, F. I., & Portegies Zwart, S. 2012, MNRAS, 420, 1503 [NASA ADS] [CrossRef] [Google Scholar]
 Pelupessy, F. I., Jänes, J., & Portegies Zwart, S. 2012, New Ast., 17, 711 [Google Scholar]
 Pelupessy, F. I., van Elteren, A., de Vries, N., et al. 2013, A&A, 557, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Portegies Zwart, S., & McMillan, S. 2018, Astrophysical Recipes: The Art of AMUSE (Bristol: IOP Publishing) [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., & Verbunt, F. 2012, Astrophysics Source Code Library [record ascl:1201.003] [Google Scholar]
 Portegies Zwart, S., McMillan, S. L. W., van Elteren, E., Pelupessy, I., & de Vries, N. 2013, Comput. Phys. Commun., 183, 456 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S., Arjen, V., Inti, P., et al. 2018, AMUSE: the Astrophysical Multipurpose Software Environment, DOI: 10.5281/zenodo.1443252 [CrossRef] [Google Scholar]
 Puls, J., Sundqvist, J. O., Najarro, F., & Hanson, M. M. 2009, AIP Conf. Proc., 1171, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Puls, J., Sundqvist, J. O., & Markova, N. 2015, New Windows Massive Stars, 307, 25 [NASA ADS] [Google Scholar]
 Rimoldi, A., Portegies Zwart, S., & Rossi, E. M. 2016, Comput. Astrophys. Cosmol., 3, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, H., & Pittard, J. M. 2013, MNRAS, 431, 1337 [NASA ADS] [CrossRef] [Google Scholar]
 Saladino, M. I., Pols, O. R., van der Helm, E., Pelupessy, I., & Portegies Zwart, S. 2018, A&A, 618, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (New York: Academic) [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [NASA ADS] [CrossRef] [Google Scholar]
 Sugawara, Y., Tsuboi, Y., & Maeda, Y. 2008, A&A, 490, 259 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Taylor, G. 1950, Proc. R. Soc. London Ser. A, 201, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Theuns, T., & Jorissen, A. 1993, MNRAS, 265, 946 [NASA ADS] [CrossRef] [Google Scholar]
 van Elteren, A., Pelupessy, I., & Portegies Zwart, S. 2014, R. Soc. London Philos. Trans. Ser. A, 372, 30385 [NASA ADS] [CrossRef] [Google Scholar]
 Wall, J., McMillan, S. L. W., Mac Low, M. M., et al. 2017, AAS Meeting #229, 229, 153.02 [Google Scholar]
 Wang, J., & White, S. D. M. 2007, MNRAS, 380, 93 [NASA ADS] [CrossRef] [Google Scholar]
 White, S. D. M. 1996, in Cosmology and Large Scale Structure, eds. R. Schaeffer, J. Silk, M. Spiro, & J. ZinnJustin, 349 [Google Scholar]
Appendix A: Equations
In this appendix, we calculate the analytical predictions for a stationary, spherically symmetric wind which are used in STELLAR_WIND.PY. For these calculations, we assume that the massloss rate (Ṁ) and the velocity as a function of radius (υ(r)) are known and we define the acceleration
$$\begin{array}{c}\hfill a(r)=\frac{\mathrm{d}\upsilon}{\mathrm{d}t}=\frac{\mathrm{d}\upsilon}{\mathrm{d}r}\frac{\mathrm{d}r}{\mathrm{d}t}=\upsilon (r)\frac{\mathrm{d}\upsilon}{\mathrm{d}r}\xb7\end{array}$$(A.1)
A.1. 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 υ(r) is known, we can write
$$\begin{array}{cc}& \upsilon (r(t))=\frac{\mathrm{d}r(t)}{\mathrm{d}t},\hfill \\ \hfill & \mathrm{d}t=\frac{\mathrm{d}r}{\upsilon (r)},\hfill \end{array}$$(A.2)
which is solved by,
$$\begin{array}{c}\hfill t={\displaystyle {\int}_{{R}_{\ast}}^{r(t)}\frac{1}{\upsilon (r)}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}r}.\end{array}$$(A.3)
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 υ(r) = υ_{∞} then
$$\begin{array}{ccc}& t=\frac{1}{{\upsilon}_{\infty}}{\displaystyle {\int}_{{R}_{\ast}}^{r(t)}\mathrm{d}r}\hfill & \hfill =\frac{r(t){R}_{\ast}}{{\upsilon}_{\infty}},\\ \hfill & r(t)={R}_{\ast}+t\ast {\upsilon}_{\infty}.\hfill \end{array}$$(A.4)
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 massloss rate, Ṁ is related to the density and the velocity at any point of the wind via the equation of mass continuity,
$$\begin{array}{c}\hfill \dot{M}=4\pi {r}^{2}\rho (r)\upsilon (r),\end{array}$$(A.5)
where ρ is the density of the wind. Because we assume that Ṁ and υ(r) are known, we can rewrite this as
$$\begin{array}{c}\hfill \rho (r)=\frac{\dot{M}}{4\pi {r}^{2}\upsilon (r)}\xb7\end{array}$$(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
$$\begin{array}{c}\hfill x=\frac{{\int}_{{R}_{\ast}}^{{r}_{p}}\pi {r}^{2}\rho (r)\mathrm{d}r}{{\int}_{{R}_{\ast}}^{r(t)}\pi {r}^{2}\rho (r)\mathrm{d}r},\end{array}$$(A.7)
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
$$\begin{array}{c}\hfill {x}_{u}=\frac{{\int}_{{R}_{\ast}}^{{r}_{p}}{r}^{2}\mathrm{d}r}{{\int}_{{R}_{\ast}}^{r(t)}{r}^{2}\mathrm{d}r}=\frac{{r}_{p}^{3}{R}_{\ast}^{3}}{r{(t)}^{3}{R}_{\ast}^{3}}\xb7\end{array}$$(A.8)
For the desired density profile based on a given velocity profile, we rewrite Eq. (A.7) in terms of υ using Eq. (A.6)
$$\begin{array}{c}\hfill {x}_{\upsilon}=\frac{{\int}_{{R}_{\ast}}^{{r}_{p}}\frac{\dot{M}}{\upsilon (r)}\mathrm{d}r}{{\int}_{{R}_{\ast}}^{r(t)}\frac{\dot{M}}{\upsilon (r)}\mathrm{d}r}=\frac{{\int}_{{R}_{\ast}}^{{r}_{p}}\frac{1}{\upsilon (r)}\mathrm{d}r}{{\int}_{{R}_{\ast}}^{r(t)}\frac{1}{\upsilon (r)}\mathrm{d}r}\xb7\end{array}$$(A.9)
We then set x_{u} = x_{υ} where x_{u} is calculated with the old particle radius (of the generated uniform density shell). The last step is to solve Eq. (A.9) 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 υ(r) = υ_{∞} then
$$\begin{array}{cc}& x=\frac{{\int}_{{R}_{\ast}}^{{r}_{p}}\frac{1}{{\upsilon}_{\infty}}\mathrm{d}r}{{\int}_{{R}_{\ast}}^{r(t)}\frac{1}{{\upsilon}_{\infty}}\mathrm{d}r}=\frac{{\int}_{{R}_{\ast}}^{{r}_{p}}\mathrm{d}r}{{\int}_{{R}_{\ast}}^{r(t)}\mathrm{d}r}=\frac{{r}_{p}{R}_{\ast}}{r(t){R}_{\ast}},\hfill \\ \hfill & {r}_{p}={R}_{\ast}+x(r(t){R}_{\ast}).\hfill \end{array}$$(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,
$$\begin{array}{c}\hfill P=K\rho {(r)}^{\gamma},\end{array}$$(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
$$\begin{array}{c}\hfill P(r)=(\gamma 1)\rho (r)u,\end{array}$$(A.12)
where u is the internal energy of the gas particles defined by
$$\begin{array}{c}\hfill u=\frac{3}{2}\frac{{k}_{B}{T}_{\ast}}{\mu},\end{array}$$(A.13)
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 Eqs. (A.11) and (A.12) we get
$$\begin{array}{c}\hfill K=u(\gamma 1)\rho {({R}_{\ast})}^{1\gamma}.\end{array}$$(A.14)
The acceleration caused by the gradient of the gas pressure is
$$\begin{array}{c}\hfill {a}_{P}(r)=\frac{1}{\rho (r)}\frac{\partial P(r)}{\partial r},\end{array}$$(A.15)
which we can rewrite using Eqs. (A.11) and (A.6)
$$\begin{array}{cc}\hfill {a}_{P}(r)& =\frac{K}{\rho (r)}\frac{\partial {\rho}^{\gamma}}{\partial r}\hfill \\ \hfill & =\frac{K}{\rho (r)}\gamma \rho {(r)}^{\gamma 1}\frac{\partial}{\partial r}\frac{\dot{M}}{4\pi {r}^{2}\upsilon (r)}\hfill \\ \hfill & =K\gamma \rho {(r)}^{\gamma 1}(\frac{2}{r}+\frac{1}{\upsilon (r)}\frac{\mathrm{d}\upsilon (r)}{\mathrm{d}r})\xb7\hfill \end{array}$$(A.16)
A.4. Density as a function of radius
In Sects. 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.
All Tables
Overview of the modes in STELLAR_WIND.PY and their suggested application domains.
All Figures
Fig. 1. 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. 

In the text 
Fig. 2. Radial velocity profiles for the wind acceleration functions currently available in STELLAR_WIND.PY. The formulae for these functions can be found in 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 υ_{0} = 0.2 υ_{∞}, 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 the text 
Fig. 3. Diagram of the AMUSE codes (boxes) and their relations (labelled arrows). Dotted lines indicate optional codes that can be added depending on the astrophysical application. 

In the text 
Fig. 4. Analytical and simulated velocity (υ), 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. 

In the text 
Fig. 5. Same as Fig. 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. 

In the text 
Fig. 6. 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. 

In the text 
Fig. 7. Same as Fig. 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. 

In the text 
Fig. 8. Same as Fig. 7 but for the accelerating wind mode with four different acceleration functions. 

In the text 
Fig. 9. Same as Fig. 5 but for the slow wind test. 

In the text 
Fig. 10. Gas density in the x − y plane after 0.2 Myr (top) and after 0.6 Myr (bottom) for the embedded star simulation with M_{SPH} = 0.1 M_{⊙} and ρ_{gas} = 100 M_{⊙} pc^{−3}. The embedded star is positioned at the origin (yellow dot) and the red dashed circle shows the radius with the largest mean density. 

In the text 
Fig. 11. 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}). Top panel: 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}. Bottom panel: 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). 

In the text 
Fig. 12. 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 selfsimilar blast wave in a uniform medium (Taylor 1950; Sedov 1959). 

In the text 
Fig. 13. 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. 

In the text 
Fig. 14. 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. Left hand panels: WolfRayet star (star 1, see Table 6) can be seen on the right and star 2 on the left. 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. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.