Simulating stellar winds in AMUSE

TL;DR

This paper introduces a versatile module for simulating stellar winds within the AMUSE framework, offering multiple methods that accurately reproduce wind profiles and bubble expansions, applicable to various astrophysical scenarios.

Contribution

The paper presents a new stellar_wind.py module with three simulation modes for stellar winds in SPH, improving modeling flexibility and accuracy in the AMUSE environment.

Findings

Simple and accelerating modes reproduce desired wind profiles for fast winds.

Heating mode successfully simulates wind bubbles and supernova explosions.

The module can handle complex systems like colliding stellar winds in multiple star systems.

Abstract

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.

Figures

Tables

Table 1. Table 2: Stellar and wind parameters used in the fast wind test. Derived parameters are indicated with an arrow ( → → \rightarrow ). Since the smoothing length is highly variable with extreme outliers, we include the median value of all gas particles shown in Figure  4 .

name	parameter	value
mass-loss rate	$\dot{M}$	${10}^{-6}$ $M_{\odot }$/yr
terminal wind velocity	$v_{\infty }$	700 km/s
initial wind velocity	$v_0$	100 km/s
stellar mass	$M_{\ast }$	20 $M_{\odot }$
stellar radius	$R_{\ast }$	30 $R_{\odot }$
stellar luminosity	$L_{\ast }$	100 000 $L_{\odot }$
stellar surface temperature	$T_{\ast }$	20 000 K
escape velocity at $R_{\ast }$	$v_{\mathrm{esc}}\left(R_{\ast }\right)$	$\to 504.5$ km/s
wind timestep	$t_{\mathrm{wind}}$	0.02 days
end time	$t_{\mathrm{end}}$	5 days
SPH particle mass	$M_{\mathrm{SPH}}$	${10}^{-11}$ $M_{\odot }$
particles per shell	$N_{\mathrm{shell}}$	$\to \sim 5$
particles in simulation	$N_{\mathrm{tot}}$	$\to \sim 1378$
median smoothing length	$h$	$\to \sim 32$ $R_{\odot }$

Table 2. Table 3: Stellar and wind parameters used in the slow wind test. Derived parameters are indicated with an arrow ( → → \rightarrow ). Since the smoothing length is highly variable with extreme outliers, we include the median value of all gas particles shown in Figure  8 .

name	parameter	value
mass-loss rate	$\dot{M}$	$5\cdot {10}^{-7}$ $M_{\odot }$/yr
terminal wind velocity	$v_{\infty }$	25 km/s
initial wind velocity	$v_0$	2 km/s
stellar mass	$M_{\ast }$	2 $M_{\odot }$
stellar radius	$R_{\ast }$	300 $R_{\odot }$
stellar luminosity	$L_{\ast }$	8 000 $L_{\odot }$
stellar surface temperature	$T_{\ast }$	3 000 K
escape velocity at $R_{\ast }$	$v_{\mathrm{esc}}$($R_{\ast }$)	$\to$ 50.45 km/s
wind timestep	$t_{\mathrm{wind}}$	2 days
end time	$t_{\mathrm{end}}$	2 000 days
SPH particle mass	$M_{\mathrm{SPH}}$	${10}^{-9}$ $M_{\odot }$
particles per shell	$N_{\mathrm{shell}}$	$\to \sim 5$
particles in simulation	$N_{\mathrm{tot}}$	$\to \sim 8476$
median smoothing length	$h$	$\to \sim 159$ $R_{\odot }$

Table 3. Table 4: The parameters used in the embedded star test. The stellar and wind parameters are the same as in Table  2 .

name	parameter	value
gas density	${\rho }_{\mathrm{gas}}$	10 and 100 $M_{\odot }$/pc3
wind timestep	$t_{\mathrm{wind}}$	$2\cdot {10}^4$ yr
end time	$t_{\mathrm{end}}$	${10}^6$ yr
SPH particle mass	$M_{\mathrm{SPH}}$	0.05 to 1 $M_{\odot }$
wind release radius	$r_{\mathrm{wind}}$	0.01 pc
feedback efficiency	$f_{\mathrm{fb}}$	0.01
particles per shell	$N_{\mathrm{shell}}$	$\to$ 0 or 1
new particles in simulation	$N_{\mathrm{new}}$	$\to$ 5 to 100
median smoothing length	$h$	$\to \sim 0.2$ pc

Table 4. Table 5: The parameters used in the supernova test.

name	parameter	value
initial stellar mass	$M_{\mathrm{ZAMS}}$	20 $M_{\odot }$
stellar age	$T_{\ast }$	9.78 Myr
gas density	${\rho }_{\mathrm{gas}}$	10 $M_{\odot }$/pc3
end time	$t_{\mathrm{end}}$	2000 yr
SPH particle mass	$M_{\mathrm{SPH}}$	0.01  – 1 $M_{\odot }$
wind timestep	$t_{\mathrm{wind}}$	20 yr
wind release radius	$r_{\mathrm{wind}}$	0.01 pc
supernova energy	$E_{\mathrm{SN}}$	${10}^{51}$ erg
mass-loss	$\Delta M_{\ast }$	$12.95$ $M_{\odot }$
feedback efficiency	$f_{\mathrm{fb}}$	0.01

Table 5. Table 7: An overview of the modes in stellar_wind.py  and their suggested application domains.

mode	section	description	application
simple	2.2	Creates particles with a radial velocity given by the desired terminal wind velocity.	Detailed wind interaction simulations well outside the acceleration zone and past the sonic point.
accelerating	2.3	Similar to simple wind, but also accelerate particles near the star.	Detailed wind simulation near or inside the acceleration zone and near the sonic point.
heating	2.4	Does not give new particles a radial velocity, but instead adds internal energy to the particles.	Large scale, low resolution simulations of wind from embedded stars, including the effect of a supernova.

Equations

$$\begin{array}{rl}{v}_{\infty }& =C(T_{\ast })v_{\mathrm{phesc}},\\ \mathrm{where},& \\ C(T_{\ast })& =\{\begin{array}{ll}1& T_{\ast }\le 10000\text{leavevmode}\mathrm{K},\\ 1.4& 10000\text{leavevmode}\mathrm{K}<T_{\ast }<21000\text{leavevmode}\mathrm{K},\\ 2.65& T_{\ast }\ge 21000\text{leavevmode}\mathrm{K},\end{array}\\ v_{\mathrm{phesc}}& =\sqrt{2g_{\ast }{R}_{\ast }\left(1-\Gamma \right)},\\ g_{\ast }& =\frac{G{M}_{\ast }}{R_{\ast }^2},\\ \Gamma & =7.66\cdot 10^{-5}{\sigma }_e\frac{L_{\ast }/{\mathrm{L}}_{\odot }}{M_{\ast }/{\mathrm{M}}_{\odot }},\\ {\sigma }_e& =0.398\frac{1+I_{\mathrm{He}}\mathrm{Y}}{1+4\mathrm{Y}},\end{array}$$

$$\begin{array}{rl}{E}_{\mathrm{mech}}& ={\int }_{t_0}^{t_1}{L}_{\mathrm{mech}}(t)dt,\\ L_{\mathrm{mech}}(t)& =\frac{1}{2}\dot{M}(t)v_{\infty }(t{)}^2,\end{array}$$

$$u=f_{\mathrm{fb}}\frac{E_{\mathrm{mech}}}{\Delta M_{\ast }},$$

$$a(r)=\frac{dv}{dt}=\frac{dv}{dr}\frac{dr}{dt}=v(r)\frac{dv}{dr}.$$

$$\begin{array}{rl}v(r(t))& =\frac{dr(t)}{dt},\\ dt& =\frac{dr}{v(r)},\end{array}$$

$$t={\int }_{R_{\ast }}^{r(t)}\frac{1}{v(r)}dr.$$

$$\begin{array}{rl}t& =\frac{1}{v_{\infty }}{\int }_{R_{\ast }}^{r(t)}dr=\frac{r(t)-R_{\ast }}{v_{\infty }},\\ r(t)& =R_{\ast }+t\ast v_{\infty }.\end{array}$$

$$\dot{M}=4\pi r^2\rho (r)v(r),$$

$$\rho (r)=\frac{\dot{M}}{4\pi r^{2}v(r)}.$$

$$x=\frac{{\int }_{R_{\ast }}^{r_p}\pi r^2\rho (r)dr}{{\int }_{R_{\ast }}^{r(t)}\pi r^2\rho (r)dr},$$

$$x_u=\frac{{\int }_{R_{\ast }}^{r_p}{r}^{2}dr}{{\int }_{R_{\ast }}^{r(t)}{r}^{2}dr}=\frac{r_p^3-R_{\ast }^3}{r(t{)}^3-R_{\ast }^3}.$$

$$x_v=\frac{{\int }_{R_{\ast }}^{r_p}\frac{\dot{M}}{v(r)}dr}{{\int }_{R_{\ast }}^{r(t)}\frac{\dot{M}}{v(r)}dr}=\frac{{\int }_{R_{\ast }}^{r_p}\frac{1}{v(r)}dr}{{\int }_{R_{\ast }}^{r(t)}\frac{1}{v(r)}dr}.$$

$$\begin{array}{rl}x& =\frac{{\int }_{R_{\ast }}^{r_p}\frac{1}{v_{\infty }}dr}{{\int }_{R_{\ast }}^{r(t)}\frac{1}{v_{\infty }}dr}=\frac{{\int }_{R_{\ast }}^{r_p}dr}{{\int }_{R_{\ast }}^{r(t)}dr}=\frac{r_p-R_{\ast }}{r(t)-R_{\ast }},\\ r_p& =R_{\ast }+x(r(t)-R_{\ast }).\end{array}$$

$$P=K\rho (r{)}^{\gamma },$$

$$P(r)=(\gamma -1)\rho (r)u,$$

$$u=\frac{3}{2}\frac{k_B{T}_{\ast }}{\mu },$$

$$K=u(\gamma -1)\rho (R_{\ast }{)}^{1-\gamma }.$$

$$a_P(r)=-\frac{1}{\rho (r)}\frac{\partial P(r)}{\partial r},$$

$$\begin{array}{rl}{a}_P(r)& =-\frac{K}{\rho (r)}\frac{\partial {\rho }^{\gamma }}{\partial r}\\ & =-\frac{K}{\rho (r)}\gamma \rho (r{)}^{\gamma -1}\frac{\partial }{\partial r}\frac{\dot{M}}{4\pi r^{2}v(r)}\\ & =K\gamma \rho (r{)}^{\gamma -1}\left(\frac{2}{r}+\frac{1}{v(r)}\frac{dv(r)}{dr}\right).\end{array}$$

Full text

*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.

1 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, $\dot{M}$, and the terminal wind velocity, $v_{\infty}$, 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 ($\dot{M}\mathord{\sim}10^{-14}$ $\mathrm{M}_{\odot}$/yr) that are thermally or gas pressure driven.
2. Cool giants and super giants have slow ($v_{\infty}\mathord{\sim}5\,\text{--}\,30$ km/s) high mass-loss rate ($\dot{M}\mathord{\sim}10^{-7}\,\text{--}\,10^{-5}$ $\mathrm{M}_{\odot}$/yr) winds driven mainly by radiation pressure on dust (Höfner 2015).
3. Hot luminous stars have fast ($v_{\infty}\mathord{\sim}200\,\text{--}\,3000$ km/s) high mass-loss rate ($\dot{M}\mathord{\sim}10^{-7}\,\text{--}\,10^{-5}$ $\mathrm{M}_{\odot}$/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 (amuse111amusecode.org; 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. 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 SPH (Smoothed particle hydrodynamics), N-body, 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 Section 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 particles222A 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., 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_{\mathrm{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).

2.1 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_{\infty}$, is not calculated by any code currently in amuse. Determining $v_{\infty}$ 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:

[TABLE]

where $v_{\mathrm{phesc}}$ is the photospheric escape velocity (similar to the escape velocity $v_{\mathrm{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, $\Gamma$ is the ratio of radiative Thomson acceleration to gravitational acceleration, $\sigma_{e}$ is the Thomson absorption coefficient, $Y$ is the Helium fraction and $I_{\mathrm{He}}$ is the number of electrons per Helium nucleus (in this paper we use default values of $I_{\mathrm{He}}=2$ and Y = 0.25). For cooler stars, $v_{\infty}\approx v_{\mathrm{phesc}}$ and this formula is still applicable (Kudritzki, private communication).

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, $v(r)=v_{\infty}$ 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)=\sqrt{{v_{\infty}}^{2}+v_{\mathrm{esc}}(r)^{2}}$ where $v_{\mathrm{esc}}(r)=\sqrt{2GM_{*}/r}$ is the local escape velocity at the initial particle radius, $r$. We calculate this new velocity for each particle because $v_{\mathrm{esc}}(r)$ can vary within the thin shell in which we create the particles.

We set the outer radius ($r_{\mathrm{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 $\delta 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 function333In the SPH code used in this paper we use the spline kernel. 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 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_{\mathrm{desired}}=\delta t\cdot\dot{M}/M_{\mathrm{SPH}}$. We therefore remove a number of randomly selected particles from the grid (typically $\mathord{\sim}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.

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, $v<v_{\infty}$. 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 timestep444 This should not be confused with the internal timesteps for each code which may be variable, meaning particles have different timesteps depending on conditions such as local particle density.. The mutual gravitational effect is included by 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 kick-drift-kick scheme to ensure a smooth acceleration of the wind particles.

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 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_{\mathrm{mid}}$ and either s or $\alpha$ (for the logistic and agb function, respectively). These parameters determine the center and the width of the acceleration zone. The default values $r_{\mathrm{mid}}=3$ and $s=\alpha=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 $\beta$ parameter indicates the steepness of the acceleration curve and is often derived from observations. The example values $\beta=0.8$ and $\beta=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_{\mathrm{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 12 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_{\mathrm{c}}$. When new wind particles are created, $\delta t$, which determines the thickness of the shell of new particles, is compared to $t_{\mathrm{c}}$. If $\delta t<$ $t_{\mathrm{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 12 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).

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_{\mathrm{mech}}$) of the accumulated wind, defined as,

[TABLE]

where $L_{\mathrm{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

[TABLE]

where $\Delta M_{*}$ is the mass lost and $f_{\mathrm{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 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 time-stepping artefacts (e.g. Pelupessy & Portegies Zwart 2012). One way to prevent this is to use a very small timestep (preliminary tests suggest $\mathord{\sim}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 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 $\alpha=0.5$ and $\beta=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 fastkick555The 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.. When we simulate multiple stars that evolve dynamically (Section 3.5), we use the N-body 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_{\mathrm{br}}$), the (maximum) internal timestep ($t_{\mathrm{N-body}}$) of the SPH and N-body codes and the wind release timestep ($t_{\mathrm{wind}}$), as well as the end time ($t_{\mathrm{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_{\mathrm{br}}\geq t_{\mathrm{N-body}}$ and the wind code requires $t_{\mathrm{wind}}\geq t_{\mathrm{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_{\mathrm{wind}}$ and choose $t_{\mathrm{wind}}=2t_{\mathrm{br}}=4t_{\mathrm{N-body}}$. For the hydro simulation we also need to set the SPH particle mass ($M_{\mathrm{SPH}}$).

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, $v_{0}=100$ km/s, is based on numerical considerations. As we shall see in Section 3.2, slow (subsonic666If the wind speed near the star is lower than the local sound speed, the wind is called subsonic. On the other hand if the wind speed is higher than the local sound speed, i.e. past the ’sonic point’, the wind is called supersonic.) winds are more complex to simulate and for many applications using a higher $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 follow the desired curve, but with more scatter. For the accelerating wind curves, we see that in regions with high acceleration the densities and temperatures in the simulation are too high. This is a result of the low resolution in combination with the way densities are calculated in SPH, using a kernel function that ’smears out’ these variables. We show in Figure 5 that for a higher resolution (smaller $M_{\mathrm{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_{\mathrm{code}}$) divided by the total cpu (or wall-clock) time ($t_{\mathrm{tot}}$) as a function of resolution777These simulations where all performed on the same desktop computer using a 4-core Intel Xeon E5507 CPU.. In the top panel we see that when using the simple wind mode, the time spent in the stellar_wind.py code is less than 1% when using more than $\mathord{\sim}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 12 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_{\mathrm{c}}$, see Section 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.

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.

In Figure 7 we present the velocity profiles of simulations using the simple wind mode and varying $v_{\infty}$ 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=1500\mathrm{R}_{\odot}$).

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 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.

3.3 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 star clusters and it is what the heating wind mode is designed for. The initial gas is distributed along a grid888As mentioned in Section 2.2.1, this can introduce preferential directions and a glass or other relaxed system should be considered for most applications. to ensure a constant density and a divergence-free random Gaussian velocity field following 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_{\mathrm{wind}}$) is set manually to $r_{\mathrm{wind}}=0.01$ pc and the feedback efficiency is set to $f_{\mathrm{fb}}=0.01$ following Pelupessy & Portegies Zwart (2012).

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_{\rho\mathrm{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 Figure 11. Note that $r_{\rho\mathrm{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_{\mathrm{SPH}}$ (different resolutions) and two different gas densities. Even when the resolution is very low ($M_{\mathrm{SPH}}=1.0\leavevmode\nobreak\ \mathrm{M}_{\odot}$), 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_{\mathrm{SPH}}=1.0\leavevmode\nobreak\ \mathrm{M}_{\odot}$ and $\dot{M}=1$ $\mathrm{M}_{\odot}$/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_{\mathrm{SPH}}$). We therefore advise that the choice of $M_{\mathrm{SPH}}$ be based on the stellar mass-loss rate and the delay and expansion profile that would be acceptable in the desired simulations.

3.4 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 2012) to a few timesteps ($\mathord{\sim}40$ yr) before the star goes supernova ($\mathord{\sim}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 $\mathrm{M}_{\odot}$ 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_{\mathrm{wind}}=20$ yr is required to avoid time-stepping artefacts with this high velocity (also see: Pelupessy & Portegies Zwart 2012).

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_{\mathrm{SPH}}=1$ $\mathrm{M}_{\odot}$), 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_{\mathrm{SPH}}=0.01$ $\mathrm{M}_{\odot}$ simulation. This feature is present at some point for all high resolution simulations and is the result of a reverse density wave within the outgoing shockwave. These waves are an artefact of the hydrodynamical simulation method used, but they are not an accurate representation of the true physical process. They should not be confused with the reverse shock that takes place in real supernova remnants. While these density waves do subside after a few thousand years, the density inside the supernova bubble shortly after the explosion should not be considered correct.

The time evolution of the expansion of gas from a supernova explosion is usually modelled in separate phases. The first phase is a free expansion, where the ejected gas moves at an approximately constant velocity, sweeping up the gas in the interstellar medium. After the mass of swept up gas is equal to the mass of the originally expelled gas, the expansion can be approximated as a pure adiabatic expansion, which is described by the self-similar Sedov-Taylor solution. Only this last phase can be simulated with the heating wind mode of 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_{\rho\mathrm{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.

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 based999 We are aware that the chosen values do not match the most up-to-date observations of the WR48 system. However, the goal of this investigation is to demonstrate the use of the stellar_wind.py code, not to explain the observed system. on WR48 ($\theta$ Muscae), which is a triple system (Sugawara et al. 2008) consisting of a WC5/WC6 + O6/O7V binary with a short period ($\mathord{\sim}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.

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.

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.

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).

Acknowledgements.

We thank N. Lützgendorf for testing and improving the simple wind mode, R. P. Kudritzki for his advice on the $v_{\infty}$ 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.

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 mass-loss rate ($\dot{M}$) and the velocity as a function of radius ($v(r)$) are known and we define the acceleration

[TABLE]

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 $v(r)$ is known, we can write

[TABLE]

which is solved by,

[TABLE]

In general this equation has to be solved numerically101010 When solving the equations mentioned here numerically, we use the python library scipy (scipy.org). For integrals we use scipy.integrate.quad and for finding a root we use scipy.optimize.brentq. See docs.scipy.org for the details of these methods. for $r(t)$, although for some velocity functions we can solve it analytically, for example if $v(r)=v_{\infty}$ then

[TABLE]

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, $\dot{M}$ is related to the density and the velocity at any point of the wind via the equation of mass continuity,

[TABLE]

where $\rho$ is the density of the wind. Because we assume that $\dot{M}$ and $v(r)$ are known, we can rewrite this as

[TABLE]

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

[TABLE]

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

[TABLE]

For the desired density profile based on a given velocity profile, we rewrite equation 10 in terms of $v$ using equation 9

[TABLE]

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 12 to get the new particle radius $r_{p}$. In general this equation has to be solved numerically, although for some velocity functions we can solve it analytically, for example if $v(r)=v_{\infty}$ then

[TABLE]

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,

[TABLE]

where $K$ is the polytropic constant and $\gamma=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

[TABLE]

where $u$ is the internal energy of the gas particles defined by

[TABLE]

where $k_{B}$ is the Boltzmann constant, $T_{*}$ is the temperature at the photosphere of the star and $\mu$ is the mean molecular weight of the gas particles. Combining equations 14 and 15 we get

[TABLE]

The acceleration caused by the gradient of the gas pressure is

[TABLE]

which we can rewrite using equations 14 and 9

[TABLE]

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_{\rho\mathrm{max}}$, we calculate this for a grid of radii and select the radius with the largest density.