# JKAS Archive

Journal of the Korean Astronomical Society - Vol. 51 , No. 2

 [ Article ] Journal of the Korean Astronomical Society - Vol. 51, No. 2, pp.37-48 Abbreviation: JKAS ISSN: 1225-4614 (Print) 2288-890X (Online) Print publication date 30 Apr 2018 Received 12 Mar 2018 Accepted 10 Apr 2018 DOI: https://doi.org/10.5303/JKAS.2018.51.2.37 THE CONTRIBUTION OF STELLAR WINDS TO COSMIC RAY PRODUCTION JEONGBHIN SEO1 ; HYESUNG KANG1 ; DONGSU RYU2 1Department of Earth Sciences, Pusan National University, 2 Busandaehak-ro, Geumjeong-gu, Busan 46241, Korea (hskang@pusan.ac.kr) 2Department of Physics, School of Natural Sciences, UNIST, 50 UNIST-gil, Ulsan 44919, Korea (ryu@sirius.unist.ac.kr) Correspondence to : H. Kang JKAS is published under Creative Commons license CC BY-SA 4.0. Funding Information ▼National Research Foundation of Korea2017R1D1A1A09000567 / 2016R1A5A1013277 / 2017R1A2A1A05071429

Abstract

Massive stars blow powerful stellar winds throughout their evolutionary stages from the main sequence to Wolf-Rayet phases. The amount of mechanical energy deposited in the interstellar medium by the wind from a massive star can be comparable to the explosion energy of a core-collapse supernova that detonates at the end of its life. In this study, we estimate the kinetic energy deposition by massive stars in our Galaxy by considering the integrated Galactic initial mass function and modeling the stellar wind luminosity. The mass loss rate and terminal velocity of stellar winds during the main sequence, red supergiant, and Wolf-Rayet stages are estimated by adopting theoretical calculations and observational data published in the literature. We find that the total stellar wind luminosity due to all massive stars in the Galaxy is about Lw ≈ 1.1 × 1041 erg s−1, which is about 1/4 of the power of supernova explosions, LSN ≈ 4.8 × 1041 erg s−1. If we assume that ~ 1 − 10 % of the wind luminosity could be converted to Galactic cosmic rays (GCRs) through collisonless shocks such as termination shocks in stellar bubbles and superbubbles, colliding-wind shocks in binaries, and bow-shocks of massive runaway stars, stellar winds might be expected to make a significant contribution to GCR production, though lower than that of supernova remnants.

 Keywords: stars: massive, stars: mass-loss, stars: winds, outflows, acceleration of particles, shock waves, cosmic rays

1. INTRODUCTION

Massive stars with MZAMS ≳ 10M deposit a significant amount of mechanical energy into the interstellar medium (ISM) through stellar winds during main sequence (MS), red supergiant (RSG), and Wolf-Rayet (WR) stages, as well as supernova (SN) explosions at the end of their lives (e.g., Ekström et al. 2012; Georgy et al. 2012; Smith 2014; Yoon 2015). Hereafter, MZAMS refers to the stellar mass at zero-age MS (ZAMS). At each stage of their lives, the mechanical power of stellar winds due to a massive star can be characterized by the mass loss rate, M˙ , and the terminal velocity, ν. For a star with 35M, for instance, (1) M˙ ~ 10−6 M yr−1 and ν ~ 3, 000 km s−1 during the MS stage, (2) M˙ ~ 10−4 M yr−1 and ν ~ 25 km s−1 during the RSG stage, (3) M˙ ~ 3 × 10−5 M yr−1 and ν ~ 5, 000 km s−1 during the WR stage (Garcia-Segura et al. 1996b; Smith 2014). The so-called wind luminosity, Lw ≡ (1/2) M˙v2, ranges ~ 1034 − 1038 erg s−1 during the MS stage for MZAMS ≈ 15 − 120 M (Georgy et al. 2013). The uncertainties on mass loss, along with convective overshooting, rotation, and magnetic fields, greatly hinder full understanding of the evolution of massive stars (e.g., Yoon et al. 2010; Ekström et al. 2012; Georgy et al. 2012). For example, depending on their rotation, stars with ≳ 20−25 M are expected to experience the WR phase, but those with ≳ 32−40 M enter the WR phase without going through the RSG phase (Georgy et al. 2012).

The interaction of stellar winds with the circumstellar medium has been studied extensively by hydrodynamic simulations that are equipped with the mass loss parameters (i.e., M˙ and ν) from stellar evolution calculations (e.g., Garcia-Segura et al. 1996a,b; Freyer et al. 2003; Georgy et al. 2013). The idealized spherical structure of the so-called stellar wind bubble can be found in many previous studies such as Weaver et al. (1977) and Freyer et al. (2003). Typically, the expanding wind is terminated by a (reverse) termination shock that propagates into the wind flow, while a forward shock expands into the photo-ionized circumstellar medium. The termination shock decelerates and heats the wind gas to 106−108 K, creating a hot bubble around the central star. If we assume that an average of <Lw> ~ 1036 erg s−1 per star is deposited to the ISM through stellar winds and that there are 105 OB stars in the Galaxy, then the total wind power from massive stars is roughly in the order of Lw ~ 1041 erg s−1. On the other hand, with the SN explosion rate (SNER) of a few per century (Reed 2005) and the explosion energy of ESN ~ 1051 ergs, the energy deposition rate of SNe in the Galaxy is estimated as LSN ~ 3 × 1041 erg s−1. So the contribution of the mechanical energy from stellar winds could be comparable to that from SN explosions in the Milky Way.

Galactic cosmic rays (GCRs) with energies lower than ~ 106 GeV/nucleon are thought to be accelerated mainly within the Galactic disk (see Blandford & Eichler 1987; Hillas 2005; Drury 2012, for reviews). They are transported into the Galactic halo and escape from our Galaxy through galactic winds with a roughly constant leakage rate. From the observed ratios of the secondary to primary CRs and the observed anisotropy of CR distribution, the power of CR sources required to maintain the energy density of GCRs is estimated to be 3 − 8 × 1040 erg s−1 (Strong et al. 2010; Drury 2012). Dominant source candidates that can replenish such escape of GCRs are all related to massive stars: supernova remnants (SNRs), stellar winds, and pulsar winds.

Nonthermal particles are known to be accelerated via diffusive shock acceleration (DSA) at collisonless shocks that are ubiquitous in astrophysical environments, from the Earth’s bow-shock to merger shocks in galaxies clusters (Drury 1983). The DSA efficiency at such shocks depends mainly on the shock Mach number, magnetic field obliquity angle, and strength of MHD turbulence responsible for particle scattering (e.g., Treumann 2009). For instance, CR protons are accelerated efficiently at quasi-parallel shocks, while CR electrons are accelerated preferentially at quasiperpendicular shocks (Riquelme & Spitkovsky 2011). Recent plasma simulations have indicated that about 10% of the shock kinetic energy is transferred to the energy of CR protons at strong shocks with the shock Mach number Ms ≳ 10 and ‘quasi-parallel’ magnetic field configuration (Caprioli & Spitkovsky 2014).

It is well established that CR ions with atomic charge Za can be accelerated up to the knee energy of 3×1015Za eV by strong SNR blast shocks and that approximately 10 % of the explosion energy can be transferred to GCRs in the Galaxy (see Blasi 2013; Caprioli 2015, for recent reviews). Similarly, a significant fraction of wind mechanical power, Lw, from massive stars in pre-supernova stages is expect to be converted to GCR energy through DSA at termination shocks inside stellar wind bubbles (Casse & Paul 1980; Völk & Forman 1982; Drury 1983), shocks in particle-accelerating colliding-wind binaries (PACWBs) (De Becker et al. 2017), and bow-shocks of massive runaway stars (del Valle & Romero 2012). CR electron acceleration at strong shocks produced by colliding winds in binary (or multiple) stellar systems have been observed (De Becker & Raucq 2013). In addition, bow-shocks produced by the interaction between strong winds from massive runaway stars and the ISM may provide a minor contribution to GCR production as well (del Valle et al. 2015). On the other hand, pulsar winds are ultra-relativistic plasmas composed of primarily electron-positron pairs, and termination shocks are predominantly perpendicular with toroidal magnetic fields. Thus, mainly CR leptonic components are thought to be accelerated via shock drift acceleration and/or magnetic reconnection at pulsar wind termination shocks (Amato 2014; Sironi & Cerutti 2017).

In stellar wind bubbles around massive stars, typical termination shocks with ν ~ 3 × 103 km s−1 are strong shocks that stop the wind flow with ρwr−2 and Tw ~ 104 K, photo-ionized by the central star (Weaver et al. 1977). Yet, contrary to SNRs, the conversion of wind energy to CRs at termination shocks has not been estimated quantitatively. Unlike SNR blast waves with their relatively well preserved spherical symmetry, the structures inside wind bubbles are complex with the termination shock, contact discontinuity, and multiple shells due to winds at different stages. They are prone to several instabilities such as Rayleigh-Taylor, thin-shell, and thermal instabilities (Garcia-Segura et al. 1996b). Moreover, the wind flow may contain strong MHD turbulence and density clumps. So it is difficult to estimate the key physical parameters for the DSA process such as the radius, lifetime, magnetic field strength and obliquity of the termination shock as well as MHD turbulence of the preshock wind plasma. Moreover, if the star moves relative to the ISM (Meyer et al. 2014), or if the circumstellar medium is not uniform, then the simple spherical geometry of the wind bubble is distorted. Thus, it is very challenging to calculate the overall DSA efficiency at such complex and turbulent structures. As far as we know, an accurate estimation for the CR conversion rate of wind mechanical energy at the termination shock is not available in the literature.

Another complication is the fact that massive stars form not in isolation, but as binary systems, OB associations, or stellar clusters inside dense molecular clouds, which may lead to interacting multiple winds or superbubbles (e.g., De Becker 2007; Zinnecker & Yorke 2007; van Marle et al. 2012). The DSA efficiency at shocks formed in the wind-wind interaction region of PACWBs is likely to be higher than that at termination shocks around individual member stars. However, it depends on the fraction of the shock surface area of the windwind interaction region, in addition to the uncertain DSA parameters such as magnetic field obliquity and preshock turbulence (De Becker & Raucq 2013). Note that PACWBs have been confirmed observationally by nonthermal synchrotron radiation, although direct observational evidence for nonthermal emission from individual wind bubbles has yet to be established.

The bulk of core-collapse SNe are observed to be clustered in superbubbles, which are created by previous episodes of stellar winds and SN explosions inside OB associations (Higdon & Lingenfelter 2013). It has been shown that the observed composition of GCRs can be explained by the CR acceleration at SNRs expanding inside metal-enriched superbubbles (e.g., Higdon et al 1998; Binns et al. 2005; Bykov 2014). In addition, gamma-ray emission due to pion production in pp collisions in the Cygnus cocoon, as detected by the Fermi-LAT telescope, is interpreted as the first direct evidence for the CR acceleration in a superbubble with OB-star complexes (Ackermann et al. 2011). Obviously, it would be very difficult to quantify the CR acceleration by shocks associated with stellar winds, separately from SNRs, inside these complex superbubbles. Hence, these issues are not addressed here.

Considering these issues, here we assume that approximately 1 − 10% of the wind luminosity could be transferred to GCRs at various shocks associated with massive stars including termination shocks inside stellar bubbles and superbubbles, shocks formed by colliding winds in multiple star systems, and bow-shocks around runaway stars.

In this study, we attempt to estimate the relative importance of wind mechanical energy deposition at different stages, i.e., MS, RSG, and WR phases. To that end, we first model the wind luminosity at different stages as a function of stellar mass by adopting theoretical estimates or observational data for M˙ and ν. Adopting the galaxy-wide initial mass function of massive stars, we then estimate the number of massive stars existing in the Galactic disk and their kinetic energy contribution due to stellar winds as a function of stellar mass. Finally, we compare the total wind power deposited from stellar winds to the SNR explosion power in the Milky Way.

In the next section, we describe how we model the integrated Galactic initial mass function, mass loss rate, terminal velocity, and luminosity of stellar winds at different stellar types. In Section 3, we calculate the integrated wind power due to all massive stars in the Galaxy. Section 4 presents a brief summary.

2. MODELS

In this section, we define the mass distribution function of massive stars in the Galaxy. Then we explain how we model the mass loss rate and the terminal velocity of stellar winds at different stages. We use them to estimate wind mechanical energy deposition as a function of stellar mass.

2.1. Integrated Galactic Initial Mass Function

The initial mass function (IMF), ξ(m) = dn/dm, describes the observed distribution of the initial mass of stars in groups such as stellar clusters (see Kroupa et al. 2013, for a review). Here ξ(m)dm is the number of stars in the unit volume whose initial mass, mMZAMS, lies between m and m + dm. It is very difficult to predict the IMF theoretically, since star formation involves the complex interplay of many physical processes including gravity, hydrodynamics, radiative transfer, turbulence, magnetic fields, and external radiation field (see McKee & Ostriker 2007, for a review). However, the observed IMF is found to be remarkably universal in a wide range of environments, and can be represented well by the following canonical power-law form,

 ξ(m)=A⋅m-α, (1)

where the power-law index is α ≈ 2.3−2.7 for stars with 1 Mm ≤ 150 M (Salpeter 1955; Schmidt 1959; Miller 1979; Scalo 1986; Kroupa et al. 2002). Hereafter, m is expressed in units of M, so the normalization factor A is given in units of pc−3.

Massive stars form predominantly inside OB associations or stellar clusters. The observed distribution function of cluster mass, Mcl, formed in our Galaxy can also be described by a similar power-law function,

 ξclMcl=Acl⋅Mcl-β, (2)

for the cluster mass range of Mcl,minMclMcl,max (Kroupa et al. 2013). According to Weidner et al. (2013), the power-law index is β ≈ 2.0 for the star formation rate (SFR) of ~ 1 M yr−1, and the maximum cluster mass also depends on SFR as

 Mcl,maxM⊙=8.5×104⋅SFR(t)M⊙ yr-10.75 (3)

The smallest mass of observed clusters can be taken as Mcl,min ≈ 5 M.

Following Kroupa et al. (2013), we define the Integrated Galactic IMF (IGIMF) as the galaxy-wide IMF at a given time for all stars contained in the entire population of stellar clusters in the Galaxy:

 ξIGIMFm;t=∫Mcl,minMcl,max     tξm

Here, mmax(Mcl) is the maximum mass of the member stars contained in a cluster with Mcl. In general, the upper limit of the integration, Mcl,max (t), depends on time, for instance, as given in Equation (3), since SFR changes with time. So depending on the time variation of SFR and Mcl,max (t) in our Galaxy, IGIMF can have a power-law distribution steeper than the canonical IMF (Weidner et al. 2013).

The determination of the normalization factors, A and Acl, as well as the power-law indices, α and β, based on observed stellar populations is limited due to severe interstellar extinction, since massive stars are born in the Galactic disk and located mainly near spiral arms. However, those parameters can be estimated indirectly by comparing certain theoretical predictions with observed quantities. For example, the chemical composition of the ISM is a product of the sum of all star-formation events and ensuing chemical enrichment throughout the history of our Galaxy (Kroupa et al. 2013).

Here, adopting the results of Weidner et al. (2013), IGIMF is assumed to have the power-law form of ξIGIMF(m) ∝ m−2.6 for a SFR of 1 M yr−1. Thus we assume that the distribution function of all massive stars contained in the Galactic disk at the present time has the following form:

 N(m)=AOB⋅m-2.6. (5)

So N(m)dm represents the total number of stars with the initial mass between m and m+dm in the presentday Galaxy. As defined above, m = MZAMS is expressed in units of solar masses, so the normalization factor, AOB, is dimensionless.

We can estimate AOB approximately, using the fact that SNER in our Galaxy is 1-2 in 100 years. In other words, if each star heavier than 10 M explodes as a core-collapse SN after its MS lifetime, τMS(m), SNER can be calculated approximately by

 SNER≈∫10M⊙150M⊙N(m)τMS(m)dm. (6)

For the MS lifetimes we adopt the results of the stellar evolution calculation due to Schaller (1992), modelled by the following fitting form (Zakhozhay 2013):

 logτMS(yr)≈9.96-3.32logm+0.63(logm)2+0.19(logm)3-0.057(logm)4. (7)

where MS is given in units of years. Inserting Equations (5) and (7) into Equation (6) gives an estimated value of AOB ≈ 4.4−8.8×106 stars. Adopting this normalization, the total number of massive stars in the Galactic disk becomes 10M150MN(m)dm ≈ (0.93 − 1.85) × 105. This is fairly consistent with the results of Reed (2005), who predicted that the number of stars more massive than 10 M inside the solar circle is N(> 10M) ≈ 2×105 and that the Galactic SNER ≈ 1−2 per century. We will use the mass distribution function in Equation (5) with AOB ≈ 6.6 × 106 stars to estimate the relative contribution of wind mechanical energy from stars of different masses.

2.2. Massive Star Evolution

Figure 1 shows the evolution of the stellar mass, M (t), and the mass loss rate, M˙ (t), for massive stars (m = 15 − 120 M), which are adopted from the grid of the stellar evolution computation for nonrotating stars reported by Ekström et al. (2012). As mentioned in the Introduction, stellar rotation greatly affects stellar evolution, leading to the prediction of different evolutionary tracks in the Hertzsprung-Russell diagram, and resulting in different mass losses ∆M and lifetimes τ at different stages. Without rotation, for example, stars with ≳ 25 M experience the WR stage, while stars with ≳ 40 M do not go through the RSG stage. Stellar rotation reduces the former mass limit to ~ 20 M and the latter mass limit to ~ 32 M (Georgy et al. 2012). Thus, the parametrizations for stellar properties such as M˙ and ν adopted in this work should be taken as approximations with uncertainties of at least a factor of a few. Figure 1 demonstrates that the mass loss rate during the MS stage depends strongly on the initial mass with a range of 10−8 − 10−5 M yr−1. Although M˙(t) stays more or less constant during the MS stage, its time variation increases drastically afterwards.

Figure 1.
Time evolution of the stellar mass M (t) (left) and the mass loss rate M˙ (t) (right) taken from the grid of the stellar evolution computation for nonrotating stars presented by Ekström et al. (2012). Here M (t), M˙ , and t are expressed in units of M, M yr−1, and years, respectively.

The left-hand panel of Figure 2 shows the mass loss ∆Mk(m) during the MS, RSG, and WR stages for nonrotating stars, which are taken from Table 1 of Georgy et al. (2013). This was based on the stellar evolution grid of Ekström et al. (2012). Stars with 25 M, for example, lose about 1.4, 13.5, and 0.4 M during the MS, RSG, and WR phases, respectively. But stars with 120 M lose up to 100 M, i.e., about 32.5 and 68.4 M during the MS and WR phases, respectively. The solid lines show our fitting forms for ∆Mk(m):

Figure 2.
Left panel: Mass loss, ∆Mk(M), as a function of the initial mass m, where k stands for the MS, RSG, and WR stages. The solid lines show our fitting forms in Equations (8)-(10). Right panel: Lifetime, τk(yr), as a function of the initial mass m of the MS, RSG, and WR stages. The solid lines show our fitting forms in Equations (11)-(13). Both ∆Mk and τk are calculated with the grid of the stellar evolution computation of Ekström et al. (2012).

 logΔMMS(m)≈2.36logm-3.15, (8)
 logΔMRSG(m)≈0.85logm-0.24, (9)
 logΔMWR(m)≈1.64logm-1.57, (10)

where k stands for the MS, RSG, and WR stages.

Moreover, we estimate the lifetime of each stage by identifying the epoch t that corresponds to the mass coordinate at the end of each stage (e.g., Mend,MS = m−∆MMS) in the evolutionary tracks of Ekström et al. (2012). The right-hand panel of Figure 2 shows τk(m).

The solid lines show our fitting forms for τk(m):

 logτMS(m)≈-0.77logm+7.91, (11)
 logτRSG(m)≈-2.76logm+9.38, (12)
 logτWR(m)≈0.042logm+5.51, (13)

where the lifetimes are expressed in unit of years. For m ≥ 10 M, Equation (11) matches approximately the MS lifetime given in Equation (7), which is based on the stellar evolution calculation by Schaller (1992).

2.3. MS and WR Winds

Stellar winds from hot massive stars are thought to be driven by the transfer of energy and momentum from the radiation field to the atmospheric gas through atomic line transitions (e.g., Puls et al. 2008, for a review). So wind parameters such as the mass loss rate and terminal velocity can be estimated by solving the complex dynamic equations numerically with the line acceleration, which includes radiative transfer calculations (e.g., Vink et al. 2000, 2001; Krtic̆ka & Kub´at 2010; Muijres et al. 2012). Just like the star formation process, the theoretical determination of M˙ and ν is a very challenging problem, because it involves uncertain physics concerning non-LTE processes, opacity, wind clumping, magnetic fields, turbulence, and stellar rotation (Vink 2015).

The mass loss rate, M˙ = 4πr2ρ(r)ν, depends on the density ρ(r) at a radius where the wind has reached its terminal velocity, ν. This terminal velocity can be determined empirically by analyzing the P Cygni profiles of Hα and UV resonance lines, while M˙ can be estimated by adopting a wind density model (e.g., Lamers & Leitherer 1993; Lamers et al. 1995; Puls et al. 2008). Lamers et al. (1995) measured ν of stellar winds from stars of O-F types by analyzing the P Cygni profiles of UV lines. They found that the observed ratio of ν/νesc changes abruptly from 1.3 to 2.6 at Teff ≈ 2.5 × 104K. Here, νesc = (2GM(1 − Γe)/R)1/2 is the photospheric escape velocity corrected for the radiation pressure by electron scattering, parametrized by Γe. This so-called ‘bi-stability’ of winds comes from a shift in the ionization balance of iron (Fe III) that dominates the line acceleration in the lower part of the wind flow (Vink et al. 1999). Although uncertainties in the empirical values of ν/νesc are estimated to be about 30-40%, theoretical predictions indicate much larger variations of this ratio, ranging 1.0 ≲ ν/νesc ≲ 5.5 (e.g., Krtic̆ka 2014; Muijres et al. 2012).

In general, these wind parameters depend on stellar properties such as the stellar mass, M, luminosity, L, effective temperature, Teff, and metallicity, Z, at a given stage. For example, the theoretical prescription for M˙ (L, M, Teff, ν/νesc,Z) presented by Vink et al. (2001) (VKL01) is widely adopted in stellar evolution calculations (e.g., Ekström et al. 2012). Note that M is usually smaller than the initial mass, m, due to mass loss. Moreover, the theoretically predicted values of both M˙ and ν are expected to change significantly with time during the RSG and WR stages, while they are relatively constant during the MS stage (e.g., Freyer et al. 2003; Georgy et al. 2012). As a result, the wind parameters may not be represented accurately by some simple functions of M only, and the distribution of their observed values may have substantial variances for a given mass range. So here we attempt to model ‘timeaveraged’ values of M˙(M) and ν(M) as functions of M in order to obtain the wind mechanical power at a given evolutionary stage.

For MS winds, we extract M˙ at ZAMS from the stellar evolution grid for non-rotating stars with Z = 0.014 of Ekström et al. (2012), which is based on the VKL01 recipe. We also take more recent theoretical estimations for M˙ given in Table 2 of Krtic̆ka (2014) for B-type stars and Table 1 of Muijres et al. (2012) for O-type stars. The left-hand panel of Figure 3 compares M˙ for MS stars taken from these three references, where M˙ is given in units of M yr−1. For the case of Muijres et al. (2012), we choose the results based on the VKL01 recipe, which are larger by a factor of 2-3 than the values based on their different wind models. The predictions from both Krtic̆ka (2014) and Muijres et al. (2012) are larger by a factor of two or so than those of Ekström et al. (2012). This demonstrates the levels of uncertainties in the theoretical predictions for the wind parameters. Below we adopt the fitting form for the results of Ekström et al. (2012) (green solid line).

Figure 3.
Left panel: Mass loss rate M˙ of MS winds as a function of the MS mass MMS. Theoretical predictions are taken from Ekström et al. (2012) (green dots), Krtic̆ka (2014) (blue dots), and Muijres et al. (2012) (red dots). Estimates of Ekström et al. (2012) and Muijres et al. (2012) were calculated based on the VKL01 recipe. Right panel: Mass loss rate M˙ of WR winds as a function of the WR mass MWR. Observational data for WN (yellow dots) and WC (green dots) stars are taken from Nugis & Lamers (2000). The solid lines show our fitting formulas in Equations (14)-(16).

The theoretical predictions for ν given in Krtic̆ka (2014) and Muijres et al. (2012) tend to have somewhat large variations in the ratio of ν/νesc. So we adopt their estimates of νesc and use the empirical relations suggested by Lamers et al. (1995), that is, ν = 2.6νesc for stars earlier than B1 and ν = 1.3νesc for stars later than B1. Here, the spectral type B1 corresponds to Teff ≈ 2.5 × 104 K and M ≈ 12 M. The left-hand panel of Figure 4 shows νesc for MS stars taken from those two references and ν with the aforementioned bi-stability at M ≈ 12 M, where the velocities are given in units of km s−1.

Figure 4.
Left panel: Photospheric escape velocity, νesc, and terminal velocity, ν, for MS stars as a function of the MS mass MMS. Theoretical predictions for νesc are taken from Krtic̆ka (2014) (purple dots) and Muijres et al. (2012) (green dots), while ν = 1.3νesc (blue dots) for M > 12 M and ν = 2.6νesc (red dots) for M > 12 M. Right-hand panel: Terminal velocity ν of WR winds as a function of the WR mass MWR. Observational data for WN (blue dots) and WC (green dots) stars are taken from Nugis & Lamers (2000). The solid lines show our fitting formulas in Equations (17)-(20).

For WR winds, on the other hand, we take the observational data for M˙ and ν in Table 5 for WN stars and Table 6 for WC stars reported by Nugis & Lamers (2000). They are shown in the right-hand panels of Figures 3 and 4. As mentioned above, the observed data points for theWR stage exhibit significant scatter when plotted as a function of M because of the time-variation of the wind parameters.

The solid lines in Figure 3 show our polynomial fitting forms for M˙ :

 logM.MS≈-3.38(logM)2+14.59logM-20.84, (14)
 logM.WN≈2.17logM-7.27, (15)
 logM.WC≈1.32logM-6.15, (16)

where M represents the mass at a given stage, i.e., MMS and MWR at the MS and WR stages, respectively. Similarly, the solid lines in Figure 4 show our linear fitting forms for ν:

 logv∞,B≈0.21logM+2.85, (17)
 logv∞,O≈0.08logM+3.28, (18)
 logv∞,WN≈-0.27logM+3.49, (19)
 logv∞,WC≈0.63logM+2.64. (20)

The bi-stability transition of ν/νesc occurs at B1 stars with ~ 12 M. But, for simplicity’s sake, we refer stars earlier than B1 as “O-type” MS stars and stars later than B1 as “B-type” MS stars hereafter. The mass loss rate varies over a wide range as a function of M: M˙ ~ 10−10−10−5 M yr−1 in the MS stage and ~ 10−5.5 − 10−4.0 M yr−1 in the WR stage. On the other hand, the terminal velocity depends only weakly on M: ν ~ 103.0 −103.4 km s−1 during the MS stage and ν ~ 103.0 − 103.8 km s−1 during the WR stage. Although both M˙ and ν for the WR stage show somewhat significant scatters by a factor of up to 3 − 6 in Figures 3 and 4, we adopt the simple fitting forms.

2.4. RSG Winds

The basic mechanism driving mass loss from cool RSG stars involves the pulsation of outer layers with ensuing dust condensation and the acceleration of dust by radiation pressure, which is yet to be fully understood (Smith 2014). The observed values range as ν ~ 10 − 40 km s−1 (with ν/νesc ~ 0.2 − 0.8) and M˙ ~ 10−6 − 10−4 M yr−1 (de Jager et al. 1988; Jura & Kleinmann 1990; Mauron & Josselin 2011; Smith 2014). Since a quantitative physical model for cool winds is not available yet, here we consider an empirical parametrization due to de Jager et al. (1988), who constructed an interpolation formula that reproduces observed values of M˙ for stars of O-M types. Later Nieuwenhuijzen & de Jager (1990) published the following slightly adjusted form:

 logM.=-7.93+1.64log⁡L+0.16log⁡M-1.61log⁡Teff, (21)

where L and M are given in units of solar values. Adopting the ML relation for RSG stars, M ≈ 0.14L0.41 (Mauron & Josselin 2011), and taking the average value of Teff ≈ 3750 K, Equation (21) can be approximated as a function of M only,

 logM.RSG≈-10.27+4.16log⁡M, (22)

for the RSG mass of 10 − 32 M.

Mauron & Josselin (2011) constructed a fitting formula for ν of the RSG stars observed in the solar neighborhood as ν ≈ 20 · (L/105)0.35 km s−1. With the ML relation for RSG stars adopted above, we obtain the following empirical parametrization:

 v∞,RSG≈1.9⋅M0.85 km s-1 . (23)
2.5. Wind Luminosity

Adopting the parametrization for M˙ and ν described in Sections 2.3 and 2.4, the wind luminosity, Lw, can be approximated as

 log⁡Lw,B≈-3.38log⁡M2+15.02log⁡M+20.36, (24)
 log⁡Lw,O≈-3.38log⁡M2+14.77log⁡M+21.21, (25)
 log⁡Lw,RSG≈5.86log⁡M+25.79, (26)
 log⁡Lw,WN≈1.63log⁡M+35.21, (27)
 log⁡Lw,WC≈2.58log⁡M+34.63, (28)

where Lw is given in units of erg s−1, and again M is the mass at a given stage in units of M.

The left-hand panel of Figure 5 shows Lw(M) at different stages. In the case of the MS phase, Lw,MS increases with increasing M from 1032 erg s−1 for 10 M to 1037.5 erg s−1 for 150 M. The discontinuous change at M ≈ 12 M is due to the aforementioned bi-stability of the ratio ν/νesc. The stellar wind power is highest at the WR phase with Lw,WR ≈ 1036 − 1039 erg s−1. Although the mass of observed WR stars ranges M ≈ 5 − 30 M, they start with the initial mass m ≈ 25 − 150 M. As can be seen in Figure 2 of Georgy et al. (2013), Lw,WR(m) > Lw,MS(m) for a given mass m.

Figure 5.
Left panel: Wind luminosity Lw,k(M) as a function the stellar mass M at different evolutionary stages, where k stands for the MS, RSG, and WR stages. Right panel: Wind luminosity of the Galaxy Nk(m) · <Lw,k(m)> as a function the the initial mass m for different evolutionary stages. Here Nk(m) and <Lw,k(m)> are the galaxy-wide IMF and the time-averaged wind luminosity, respectively, for the stars that are born with the initial mass m and are now in the k phase.

3. RESULTS
3.1. Wind Energy Deposition

Next we estimate how much mechanical energy is deposited by a star with the initial m by integrating the wind luminosity over the lifetime of a given stage:

 Ew,km=∫Lw,km,tdt, (29)

where k represents the MS, RSG, and WR stages. The wind luminosity formula, Lw,k(M), in Equations (24)- (28) is given as a function of the stellar mass M (t) at a given time. Since M (t) decreases significantly throughout the lifetime of a massive star, as shown in Figure 1, it is not straightforward to relate M (t) with its initial mass m, unless we know the mass-loss history of a specific star calculated through stellar-evolution calculations. For example, there is no a priori way to find the initial mass m of an observed WR star with MWR.

During the MS phase, fortunately, the wind luminosity is almost constant in time for m ≲ 32 M, while it increases slightly by a factor of less than two for m ≳ 32 M (see Figure 2 of Georgy et al. 2013). As a result, we can assume that the time-averaged luminosity, <Lw,MS(m)>, is similar to Lw,MS(M) during the MS stage, and so the wind energy deposition can be approximated by

 Ew,MSm≈⟨Lw,MSm⟩τMSm≈Lw,MSmτMSm, (30)

where τMS(m) is the lifetime for the MS phase in Equation (11).

During the RSG and WR stages, on the other hand, M (t) and M˙ (t) change significantly. Moreover, the wind parameters shown in Figures 3 and 4 are estimated using the stellar properties such as L, M, and Teff at a given time, instead of m. It is not possible to estimate accurately the initial mass from the observed stellar properties or to convert Lw,k(M) to Lw,k(m) for the RSG or WR phases. However, using the fact that the terminal velocity is roughly constant during these two stages (Georgy et al. 2013), the wind energy deposition can be approximated as

 Ew,km≈12v∞,k2m∫M.mdt=12v∞,k2mΔMkm, (31)

where ∆Mk(m) for the RSG and WR phases is given in Equations (9) and (10).

Considering that ν(M) depends only weakly on M, we translate the relations for ν(M) to those for ν(m) as follows. The mass at the beginning of the RSG and WR stages are extracted from the stellar evolution grid of Ekström et al. (2012) and fitted by the following forms: logMRSG,i ≈ 0.79 logm + 0.25 and logMWR,i ≈ 1.35 logm − 0.81. Assuming that MMRSG,i for the RSG phase and MMWR,i for the WR phase and inserting these relations into Equations (19), (20), and (23), we can obtain approximate relations for ν(m), which are used to calculate Ew,k(m) in Equation (31).

The left-hand panel of Figure 6 shows Ew,k(m), during the MS, RSG, and WR stages. It shows that the time-integrated energy deposition during the MS phase is greater than that during the WR phase, although Lw,WR(m) is higher than Lw,MS(m) for a given m. In other words, the MS winds are less energetic but last much longer than the WR winds. For O stars with m ≳ 80 M, Ew,MS ≈ 1051 ergs, which is comparable to the explosion energy of typical SNe. On the other hands, the RSG winds are less powerful and contribute the least mechanical energy among the three wind types. Moreover, being cool, dense, and slow winds with a relatively short lifetime, the termination shock of RSG winds is not expected to be important for production of GCRs.

Figure 6.
Left panel: Wind energy deposition, Ew,k(m) = ∫ Lw,kdt (in units of ergs), as a function the initial mass m (M) at different stages, where k stands for the MS, RSG, and WR stages. Right panel: Wind energy deposition of the Galaxy N(m) · <Ew,k(m)> (ergs) as a function the initial mass m (M) for different evolutionary stages. Here N(m) is the galaxy-wide IMF.

We assumed above the time-averaged wind luminosity, <Lw,MS(m)> ≈ Lw,MS(M) for the MS phase. For the RSG and WR phases, on the other hand, we calculate <Lw,k(m)> ≈ Ew,k(m)/τk(m) by using the estimation for Ew,k(m). We will use <Lw,k(m)> to estimate the galaxy-wide wind luminosity below.

3.2. Wind Luminosity of the Galaxy

So far, we have estimated the time-averaged wind luminosity, <Lw,k(m)>, and energy deposition, Ew,k(m), for a star with initial mass m. We now consider the same quantities for all massive stars in the Galaxy with the mass distribution, N(m), in Equation (5). Note that N(m)dm represents the number of stars in the present-day Galaxy that were born with the initial mass in the range of [m,m + dm]. Then we assume that the fraction of stars that are in the k stage is proportional to the lifetime τk(m) of each stage as fk(m) ≈ τk(m)/τ(m), where τk(m) is given in Equations (11)-(13) and τ(m) ≈ τMS + τRSG + τWR. Then the number of stars in the Galaxy at each stage can be approximated as follows:

 Nkm=Nmfkm≈AOB⋅m-2.6⋅τkmτm. (32)

This galaxy-wide mass distribution function will be used to estimate the relative contribution of stellar wind luminosity from stars with the initial mass m.

The right-hand panel of Figure 5 shows the galaxywide wind luminosity, Nk(m) · <Lw,k(m)>, from stars at different phases. For the MS winds, Lw,MS(m) increases with M monotonously, but NMS · <Lw,MS> peaks at 65 M and then decreases at higher mass due to the power-law mass function. As mentioned above, <Lw,RSG(m)> ≈ Ew,RSG/τRSG for the RSG stage, and <Lw,WR(m)> ≈ Ew,WR/τWR for the WR stage. For m ≳ 40 M, the galaxy-wide wind luminosity of the Galaxy due to O-type MS stars is higher than that due to WR stars, although <Lw,WR(m)> is higher than <Lw,MS(m)>. This is because the fraction of MS stars, fMS(m), is much larger than that ofWR stars, fWR(m). The contributions from B-type MS stars and RSG stars to the galaxy-wide wind luminosity are relatively unimportant.

On the other hand, the right-hand panel of Figure 6 shows the galaxy-wide wind energy deposition, N(m)· ·Ew,k(m), which represents the time-integrated wind mechanical energy deposited during different stages from stars that are born with the initial mass m. Due to the power-law mass distribution, both N(mEw,RSG and N(m) · Ew,WR decrease with m. Stars with m ≈ 25 − 60 M during the MS and WR stages contribute the most wind energy to the ISM. Again, the contributions from B-type MS stars and RSG stars are negligible.

The total wind luminosity emitted by all massive star in the present-day Galaxy is calculated as follows:

 Lw=∫10M⊙150M⊙NMS(m)Lw,MS(m)dm+∫10M⊙40M⊙NRSG(m)Lw,RSG(m)dm+∫25M⊙150M⊙NWR(m)Lw,WR(m)dm. (33)

We find Lw ≈ 1.1×1041 erg s−1 with the various phases contributing as: Lw,MSB ≈ 3.2×1036 erg s−1, Lw,MSO ≈ 7.3 × 1040 erg s−1, Lw,RSG ≈ 7.5 × 1036 erg s−1, and Lw,WR ≈ 4.1 × 1040 erg s−1. So O-type MS stars contribute the most wind luminosity to Lw, which is about 1/4 of the SN luminosity of LSN ≈ 4.8 × 1041 erg s−1. The galaxy-wide wind luminosity from WR stars is about 10% of LSN, while the contributions from B-type MS winds and RSG winds are insignificant.

Since Lw is somewhat smaller than LSN, the relative importance between the two processes in generating GCRs could be controlled by the CR acceleration efficiencies at shocks associated with stellar winds and SNRs. As discussed in the Introduction, about 10 % of SN explosion energy is expected to be transferred to CRs at strong SNR shocks. However, ~ 1−10 % of the wind mechanical energy might be transferred to CRs at wind termination shocks, because wind bubbles have complex and turbulent structures that include unstable contact surfaces and multiple shells. Similar DSA efficiencies are expected also for shocks in PACWBs (De Becker & Raucq 2013) and bow-shocks of runaway stars (del Valle & Romero 2012). The results of our study therefore confirm that SNRs are indeed the primary sources of GCRs.

4. SUMMARY

Massive stars are born mainly in the Galactic disk and strongly influence the surrounding ISM through photoionization, stellar winds, and SN explosions. Mass loss through stellar winds is one of several key processes that govern the evolution of massive stars, but remains to be fully elucidated. In particular, wind parameters such as mass loss rate and terminal velocity have not been determined accurately because of complex physics involved in the wind dynamics, such as non-LTE processes, wind clumping, radiative transfer and turbulence among others.

In this study, we attempt to estimate quantitatively the wind mechanical energy deposition from stars more massive than 10 M in the Galaxy by adopting the following models:

• 1. We assume the Integrated Galactic IMF (IGIMF), ξIGIMF(m) ∝ m−2.6. So the number of stars, N(m)dm, formed with the initial mass in the range [m, m+dm] can be approximated by the power-law form in Equation (5) (Kroupa et al. 2013; Weidner et al. 2013).
• 2. The mass loss rate M˙ and the wind terminal velocity ν can be expressed as functions of the stellar mass M as described in Sections 2.3 and 2.4 (Vink et al. 2000, 2001; Krtic̆ka 2014; Muijres et al. 2012; Nugis & Lamers 2000; Mauron & Josselin 2011). Then the wind luminosity Lw = (1/2) M˙v2 can be estimated.
• 3. Assuming that the number of stars in different evolutionary stages is proportional to the lifetimes of each stage, Nk(m) = N(m)τk/τ, we estimate the contribution of galaxy-wide wind mechanical luminosity from stars at the MS, RSG and WR phases.
• 4. We use the stellar evolution grid for nonrotating stars presented by Ekström et al. (2012) to relate the stellar mass M (t) at a given time with its initial mass m.

Our parametrizations for the wind parameters such as M˙ , ν, Lw, and Ew should be taken as approximations with substantial uncertainties, since theoretical models for the massive star evolution and the wind dynamics are not yet fully understood. With these caveats, we attempt to evaluate the relative importance of stellar winds at different stages. The main results of this study can be summarized as follows:

• 1. The wind luminosity Lw,k(m) at different stages increases with increasing initial mass m of stars. For a given star, the wind luminosity is strongest during the WR stage with Lw,WR ≈ 1036 − 1039 erg s−1(see Figure 5).
• 2. The time-integrated wind energy deposition, Ew,k(m), increases with increasing m. For stars with m ≳ 40 M, O-type MS winds with Ew,MS ≈ 1050 − 1051 ergs provide the greater energy than WR winds, because the MS lifetime is longer than the WR lifetime.
• 3. The galaxy-wide wind mechanical luminosity from stars at the MS, RSG and WR phases is estimated to be Lw,MSB ≈ 3.2×1036 erg s−1, Lw,MSO ≈ 7.3 × 1040 erg s−1, Lw,RSG ≈ 7.5 × 1036 erg s−1, and Lw,WR ≈ 4.1 × 1040 erg s−1, respectively. So Otype MS winds provide the greatest amount of wind mechanical power.
• 4. The galaxy-wide wind luminosity Lw ≈ 1.1 × 1041 erg s−1 is about 1/4 of the SN luminosity LSN ≈ 4.8 × 1041 erg s−1, based on 1.5 SN explosions per century in the Galaxy.
• 5. It is well established that about 10% of SN explosion energy can be transferred to CRs via strong blast waves (Caprioli & Spitkovsky 2014; Caprioli 2015). On the other hand, the CR conversion efficiency for wind mechanical energy from massive stars, through termination shocks, PACWBs, and bow-shocks of massive runaways, has not yet been estimated quantitatively (see De Becker & Raucq 2013; del Valle & Romero 2012). If we adopt ~ 1 − 10% as a somewhat conservative but educated guess, this study confirms SN explosions as the primary origin of GCRs, while winds from massive stars in pre-supernova stages can provide a significant and complementary contribution.

Acknowledgments

We thank the anonymous referee and the Editor, S.-C. Yoon, for constructive comments and suggestions. H.K. was supported by the Basic Science Research Program of the National Research Foundation of Korea (NRF) through grant 2017R1D1A1A09000567. D.R. was supported by the NRF through grants 2016R1A5A1013277 and 2017R1A2A1A05071429.

References
 1. Ackermann, M., Ajello, M., Allafort, A., et al, (2011), A Cocoon of Freshly Accelerated Cosmic Rays Detected by Fermi in the Cygnus Superbubble, Science, 334, p1103. 2. Amato, E., (2014), The Theory of Pulsar Wind Nebulae, IJMPS, 28, p1460160. 3. Binns, W. R., Wiedenbeck, M. E., Arnould, M., et al, (2005), Cosmic-Ray Neon, Wolf-Rayet Stars, and the Superbubble Origin of Galactic Cosmic Rays, ApJ, 634, p1. 4. Blandford, R. D., & Eichler, D., (1987), Particle Acceleration at Astrophysical Shocks -- a Theory of Cosmic-Ray Origin, Phys. Rep, 154, p1. 5. Blasi, P., (2013), The Origin of Galactic Cosmic Rays, A&A Rev, 21, p70. 6. Bykov, A. M., (2014), Nonthermal Particles and Photons in Starburst Regions and Superbubbles, A&A Rev, 22, p77. 7. Casse, M., & Paul, J. A., (1980), Local Gamma Rays and Cosmic-Ray Acceleration by Supersonic Stellar Winds, ApJ, 237, p236. 8. Caprioli, D., (2015), Cosmic-Ray Acceleration and Propagation, Proceedings of the 34th International Cosmic Ray Conference (ICRC2015), 34, p8. 9. Caprioli, D., & Spitkovsky, A., (2014), Simulations of Ion Acceleration at Non-Relativistic Shocks. I. Acceleration Efficiency, ApJ, 783, p91. 10. De Becker, M., (2007), Non-Thermal Emission Processes in Massive Binaries, A&A Rev, 14, p171. 11. De Becker, M., Benaglia, P., Romero, G. E., & Peri, C. S., (2017), An Investigation into the Fraction of Particle Accelerators among Colliding-wind Binaries. Towards an Extension of the Catalogue, A&A, 600, pA47. 12. De Becker, M., & Raucq, F., (2013), Catalogue of Particle-Accelerating Colliding-Wind Binaries, A&A, 558, pA28. 13. de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A., (1988), Mass Loss Rates in the Hertzsprung-Russell Diagram, A&AS, 72, p259. 14. del Valle, M. V., & Romero, G. E.., (2012), Non-Thermal Processes in Bowshocks of Runaway Stars. Application to Zeta Ophiuchi, A&A, 543, pA56. 15. del Valle, M. V., Romero, G. E., & Santos-Lima, R., (2015), Runaway Stars as Cosmic Ray Injectors inside Molecular Clouds, MNRAS, 448, p207. 16. Drury, L. O'C.., (1983), An Introduction to the Theory of Diffusive Shock Acceleration of Energetic Particles in Tenuous Plasmas, Rep. Prog. Phys, 46, p973. 17. Drury, L. O'C., (2012), Origin of Cosmic Rays, Astropart. Phys, 39, p52. 18. Dupree, A. K., (1986), Mass Loss from Cool Stars, ARA&A, 24, p377. 19. Ekström, S., Georgy, C., Eggenberger, P., et al, (2012), Grids of Stellar Models with Rotation I. Models from 0.8 to 120 ⊙ at Solar Metallicity (Z = 0.014), A&A, 537, pA146. 20. Freyer, T., Hensler, G., & Yorke, H. W., (2003), Massive Stars and the Energy Balance of the Interstellar Medium. I. The Impact of an Isolated 60 ⊙ Star, ApJ, 594, p888. 21. Garcia-Segura, G., Langer, N., & Mac Low, M.-M., (1996a), The Hydrodynamic Evolution of Circumstellar Gas around Massive Stars. II. The Impact of the Time Sequence O Star → RSG → WR Star, A&A, 316, p133. 22. Garcia-Segura, G., Mac Low, M.-M., & Langer, N., (1996b), The Dynamical Evolution of Circumstellar Gas around Massive Stars. I. The Impact of the Time Sequence O Star → LBV → WR Star, A&A, 305, p229. 23. Georgy, C., Ekström, S., Meynet, G., et al, (2012), Grids of Stellar Models with Rotation II. WR Populations and Supernovae/GRB Progenitors at Z = 0.014, A&A, 542, pA29. 24. Georgy, C., Walder, R., Folini, D., et al, (2013), Circumstellar Medium around Rotating Massive Stars at Solar Metallicity, A&A, 559, pA69. 25. Higdon, J. C., Lingenfelter, R. E., & Ramaty, R., (1998), Cosmic-Ray Acceleration from Supernova Ejecta in Superbubbles, ApJL, 509, pL33. 26. Higdon, J. C., & Lingenfelter, R. E., (2013), The Galactic Spatial Distribution of OB Associations and Their Surrounding Supernova-Generated Superbubbles, ApJ, 775, p110. 27. Hillas, A. M., (2005), Can Diffusive Shock Acceleration in Supernova Remnants Account for High Energy Galactic Cosmic Rays?, J. Phys. G, 31, pR95. 28. Jura, M., & Kleinmann, S. G., (1990), Mass-Losing M Supergiants in the Solar Neighborhood, ApJS, 73, p769. 29. Kroupa, P., & Boily, C. M., (2002), On the Mass Function of Star Clusters, MNRAS, 336, p1188. 30. Kroupa, P., Weidner, C., Pflamm-Altenburg, J., et al, (2013), The Stellar and Sub-Stellar Initial Mass Function of Simple and Composite Populations, in Planets, Stars and Stellar Systems, Volume 5: Galactic Structure and Stellar Populations, Oswalt, T. D., & Gilmore, G., Gilmore, Dordrecht, Springer, p115. 31. Krtic̆ka, J., & Kubát, J., (2010), Comoving Frame Models of Hot Star Winds. I. Test of the Sobolev Approximation in the Case of Pure Line Transitions, A&A, 519, pA50. 32. Krtic̆ka, J., (2014), Mass Loss in Main-Sequence B Stars, A&A, 564, pA70. 33. Lamers, H. J. G. L. M., & Leitherer, C., (1993), What Are the Mass-Loss Rates of O Stars?, ApJ, 412, p771. 34. Lamers, H. J. G. L. M., Snow, T. P., & Lindholm, D. M., (1995), Terminal Velocities and the Bi-Stability of Stellar Winds, ApJ, 455, p269. 35. Mauron, N., & Josselin, E., (2011), The Mass-Loss Rates of Red Supergiants and the de Jager Prescription, A&A, 526, pA156. 36. McKee, C. F., & Ostriker, E. C., (2007), Theory of Star Formation, ARA&A, 45, p565. 37. Meyer, D. M.-A., Mackey, J., Langer, N., et al, (2014), Models of the Circumstellar Medium of Evolving, Massive Runaway Stars Moving through the Galactic Plane, MNRAS, 444, p2754. 38. Miller, G., & Scalo, J. M., (1979), The Initial Mass Function and Stellar Birthrate in the Solar Neighborhood, ApJS, 41, p513. 39. Muijres, L. E., Jorick Vink, S., de Koter, A., Mller, P. E., & Langer, N., (2012), Predictions for Mass-Loss Rates and Terminal Wind Velocities of Massive O-Type Stars, A&A, 537, pA37. 40. Nieuwenhuijzen, H., & de Jager, C., (1990), Parametrization of Stellar Rates of Mass Loss as Functions of the Fundamental Stellar Parameters M, L, and R, A&A, 231, p134. 41. Nugis, T., & Lamers, H. J. G. L. M., (2000), Mass-Loss Rates of Wolf-Rayet Stars as a Function of Stellar Parameters, A&A, 360, p227. 42. Puls, J., Vink, J. S., & Najarro, F., (2008), Mass Loss from Hot Massive Stars, A&A Rv, 16, p209. 43. Reed, B. C., (2005), New Estimates of the Solar-Neighborhood Massive Star Birthrate and the Galactic Supernova Rate, AJ, 130, p1652. 44. Riquelme, M. A., & Spitkovsky, A., (2011), Electron Injection by Whistler Waves in Non-Relativistic Shocks, ApJ, 733, p63. 45. Salpeter, E. E., (1955), The Luminosity Function and Stellar Evolution, ApJ, 121, p161. 46. Scalo, J. M., (1986), The Stellar Initial Mass Function, FCPh, 11, p1. 47. Schaller, G., Schaerer, D., Meynet, G., & Maeder, A., (1992), New Grids of Stellar Models from 0.8 to 120 Solar Masses at Z = 0.020 and Z = 0.001, A&AS, 96, p269. 48. Schmidt, M., (1959), Derivation of the Initial Luminosity Function and the Past Rate of Star Formation, International Astronomical Union, Symposium, 10, p99. 49. Sironi, L., & Cerutti, B., (2017), Particle Acceleration in Pulsar Wind Nebulae: PIC Modelling, in Modelling Pulsar Wind Nebulae, Astrophysics and Space Science Library, ed. D. F. Torres, 446, p247. 50. Smith, N., (2014), Mass Loss: Its Effect on the Evolution and Fate of High-Mass Stars, ARA&A, 52, p487. 51. Strong, A, W., Porter, T. A., Digel, S. W., et al, (2010), Global Cosmic-Ray-Related Luminosity and Energy Budget of the Milky Way, ApJL, 722, pL57. 52. Treumann, R. A., (2009), Fundamentals of Collisionless Shocks for Astrophysical Application, 1. Non-Relativistic Shocks, A&A Rv, 174, p409. 53. van Marle, A. J., Meliani, Z., & Marcowith, A., (2012), A Hydrodynamical Model of the Circumstellar Bubble Created by Two Massive Stars, A&A, 541, pL8. 54. Vink, J. S., (2015), Mass-Loss Rates of Very Massive Stars, in Very Massive Stars in the Local Universe, Astrophysics and Space Science Library, ed J. S.. Vink, 412, p77. 55. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M., (1999), On the Nature of the Bi-Stability Jump in the Winds of Early-Type Supergiants, A&A, 380, p181. 56. Vink, J. S., de Koter, A., & Lamers, H. J. G .L. M., (2000), New Theoretical Mass-Loss Rates of O and B Stars, A&A, 362, p295. 57. Vink, J. S., de Koter, A., & Lamers, H. J. G .L. M., (2001), Mass-Loss Predictions for O and B Stars as a Function of Metallicity, A&A, 369, p574 (VKL01). 58. Völk, H. J., & Forman, M., (1982), Cosmic Rays and Gamma-Rays from OB Stars, ApJ, 253, p188. 59. Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R., (1977), Interstellar Bubbles. II -- Structure and Evolution, ApJ, 218, p377. 60. Weidner, C., Kroupa, P., Pflamm-Altenburg, J., & Vazdekis, A., (2013), The Galaxy-Wide Initial Mass Function of Dwarf Late-Type to Massive Early-Type Galaxies, MNRAS, 436, p3309. 61. Yoon, S.-C., (2015), Evolutionary Models for Type Ib/c Supernova Progenitors, PASA, 32, p15. 62. Yoon, S.-C., Woosley, S. E., & Langer, N., (2010), Type Ib/c Supernovae in Binary Systems. I. Evolution and Properties of the Progenitor Stars, ApJ, 725, p940. 63. Zakhozhay, V. A., (2013), Lifetimes of Stars in the Main Sequence and the Maximum Mass of Stars in the Galactic Disk, Kinematics and Physics of Celestial Bodies, 29, p195. 64. Zinnecker, H., & Yorke, H. W., (2007), Toward Understanding Massive Star Formation, ARA&A, 45, p481.