A Systematic Survey of the Effects of Wind Mass Loss Algorithms on the Evolution of Single Massive Stars

TL;DR

This study systematically compares various wind mass loss algorithms in stellar evolution models to quantify their impact on the final properties of massive stars, revealing significant uncertainties in mass predictions but less effect on observable parameters.

Contribution

It provides a comprehensive analysis of how different wind mass loss algorithms and efficiencies affect the evolution and final structure of massive stars using MESA.

Findings

Initial to final mass mapping varies by ~50% with different algorithms.

Final observable parameters are less sensitive to mass loss uncertainties for stars with ZAMS mass ≤ 30 M_.

Internal structure, including compactness, varies significantly with mass loss treatment.

Abstract

Mass loss is a key uncertainty in the evolution of massive stars. Stellar evolution calculations must employ parametric algorithms for mass loss, and usually only include stellar winds. We carry out a parameter study of the effects of wind mass loss on massive star evolution using the open-source stellar evolution code MESA. We provide a systematic comparison of wind mass loss algorithms for solar-metallicity, nonrotating, single stars in the initial mass range of $15-35\,M_\odot$. We consider combinations drawn from two hot phase algorithms, three cool phase algorithms, and two Wolf-Rayet algorithms. We consider linear wind efficiency scale factors of $1$, $0.33$, and $0.1$ to account for reductions in mass loss rates due to wind inhomogeneities. We find that the initial to final mass mapping for each zero-age main-sequence (ZAMS) mass has a $\sim 50\%$ uncertainty if all algorithm combinations and wind efficiencies are considered. The ad-hoc efficiency scale factor dominates this uncertainty. While the final total mass and internal structure of our models vary tremendously with mass loss treatment, final observable parameters are much less sensitive for ZAMS mass $\lesssim 30\,M_\odot$. This indicates that uncertainty in wind mass loss does not negatively affect estimates of the ZAMS mass of most single-star supernova progenitors from pre-explosion observations. Furthermore, we show that the internal structure of presupernova stars is sensitive to variations in both main sequence and post main-sequence mass loss. We find that the compactness parameter $\xi\propto M/R(M)$ varies by as much as $30\%$ for a given ZAMS mass evolved with different wind efficiencies and mass loss algorithm combinations. [abridged]

Tables

Table 1. Table 1: Approximate scaling of the mass loss rate with luminosity and effective temperature, log 10 ⁡ ( M ˙ ) ∝ α ​ log 10 ⁡ ( L ) + β ​ log 10 ⁡ ( T eff ) proportional-to subscript 10 ˙ 𝑀 𝛼 subscript 10 𝐿 𝛽 subscript 10 subscript 𝑇 eff \log_{10}(\dot{M})\propto\alpha\log_{10}(L)+\beta\log_{10}(T_{\mathrm{eff}}) , predicted by the considered mass loss algorithms. Appendix  A provides the full functional form of the algorithms implemented in MESA and Tab.  9 lists the scaling with all physical quantities.

	ID	$\alpha$	$\beta$
Hot	V	$2.2$	$1.0$
	K	$1.2$	$0.6$
Cool	dJ	$1.8$	$-1.7$
	NJ	$1.6$	$-1.6$
	vL	$1.0$	$-6.3$
WR	NL	$1.3$	$0.0$
	H ($L>4.5L_{\odot }$)	$1.5$	$0.0$
	H ($L\le 4.5L_{\odot }$)	$6.8$	$0.0$

Table 2. Table 2: Combinations of wind mass loss algorithms employed in this study. The temperature threshold separating the hot phase and the cool phase is T th = 15 000 subscript 𝑇 th 15000 T_{\mathrm{th}}=15\,000 (10 000) K when using the Kudritzki et al. (Vink et al. ) algorithm. The WR phase is defined using the surface (outermost computational cell) hydrogen mass fraction X s subscript 𝑋 𝑠 X_{s} , without constraints on T eff subscript 𝑇 eff T_{\mathrm{eff}} . In the text, we do not mention the WR phase algorithm if the model discussed does not enter this phase. For a description of the algorithms, see Sec.  2.1 and the appendices listed in this table. We discuss the definition of each evolutionary phase in Sec.  2.2 .

ID	Hot phase		Cool phase		WR phase
	$T_{\mathrm{eff}}\gtrsim T_{\mathrm{th}}$		$T_{\mathrm{eff}}\lesssim T_{\mathrm{th}}$
V-dJ-NL	Vink et al.	A.1	de Jager et al.	A.3	Nugis & Lamers	A.6
V-dJ-H	Vink et al.	A.1	de Jager et al.	A.3	Hamann et al.	A.7
V-NJ-NL	Vink et al.	A.1	Nieuwenhuijzen & de Jager	A.4	Nugis & Lamers	A.6
V-NJ-H	Vink et al.	A.1	Nieuwenhuijzen & de Jager	A.4	Hamann et al.	A.7
V-vL-H	Vink et al.	A.1	van Loon et al.	A.5	Hamann et al.	A.7
V-vL-NL	Vink et al.	A.1	van Loon et al.	A.5	Nugis & Lamers	A.6
K-dJ-NL	Kudritzki et al.	A.2	de Jager et al.	A.3	Nugis & Lamers	A.6
K-dJ-H	Kudritzki et al.	A.2	de Jager et al.	A.3	Hamann et al.	A.7
K-NJ-NL	Kudritzki et al.	A.2	Nieuwenhuijzen & de Jager	A.4	Nugis & Lamers	A.6
K-NJ-H	Kudritzki et al.	A.2	Nieuwenhuijzen & de Jager	A.4	Hamann et al.	A.7
K-vL-NL	Kudritzki et al.	A.2	van Loon et al.	A.5	Nugis & Lamers	A.6
K-vL-H	Kudritzki et al.	A.2	van Loon et al.	A.5	Hamann et al.	A.7

Table 3. Table 3: Model summary at the end of the mass loss phase (when T c ≥ 10 9 ​ K subscript 𝑇 𝑐 superscript 10 9 K T_{c}\geq 10^{9}\,\mathrm{K} , roughly corresponding to neon core ignition). We provide total mass M 𝑀 M , helium core mass M He subscript 𝑀 He M_{\mathrm{He}} , and carbon-oxygen (CO) core mass M CO subscript 𝑀 CO M_{\mathrm{CO}} . We omit models differing only by the WR wind scheme if they do not reach the WR stage (i.e.  X s > 0.4 subscript 𝑋 𝑠 0.4 X_{s}>0.4 at all times). We define the edge of the CO core as first location going inward where X ( 4 He ) < 0.01 X(^{4}\mathrm{He})<0.01 , without requiring a minimum mass fraction of carbon or oxygen. Analogously, we define the outer edge of the helium core as the first location going inward where X ( 1 H ) < 0.01 X(^{1}\mathrm{H})<0.01 .

End of mass loss phase: $T_c\ge {10}^{9}K$
$M_{\mathrm{ZAMS}}=15M_{\odot }$					$M_{\mathrm{ZAMS}}=20M_{\odot }$
$\eta$	ID	$M\left[M_{\odot }\right]$	$M_{\mathrm{He}}\left[M_{\odot }\right]$	$M_{\mathrm{CO}}\left[M_{\odot }\right]$	$\eta$	ID	$M\left[M_{\odot }\right]$	$M_{\mathrm{He}}\left[M_{\odot }\right]$	$M_{\mathrm{CO}}\left[M_{\odot }\right]$
0.1	V-dJ	14.66	04.99	03.20	0.1	V-dJ	19.23	07.04	4.91
	V-NJ	14.64	04.99	03.20		V-NJ	19.23	07.04	4.91
	V-vL	13.60	04.98	03.20		V-vL	18.10	07.03	4.91
	K-dJ	14.67	05.00	03.22		K-dJ	19.38	07.01	4.90
	K-NJ	14.66	05.00	03.22		K-NJ	19.37	07.01	4.90
	K-vL	13.61	05.00	03.21		K-vL	18.67	07.01	4.90
0.33	V-dJ	13.94	04.93	03.16	0.33	V-dJ	17.47	07.01	4.88
	V-NJ	13.90	04.93	03.16		V-NJ	17.48	07.01	4.88
	V-vL	10.39	04.93	03.15		V-vL	13.47	06.99	4.87
	K-dJ	13.92	04.98	03.19		K-dJ	17.62	07.00	4.87
	K-NJ	13.88	04.98	03.19		K-NJ	17.62	07.00	4.88
	K-vL	10.11	04.97	03.19		K-vL	13.90	06.98	4.87
1.0	V-dJ	12.86	04.65	02.95	1.0	V-dJ	11.81	07.06	4.92
	V-NJ	12.74	04.65	02.95		V-NJ	12.04	07.06	4.92
	V-vL	05.25	04.64	02.94		V-vL	08.80	07.02	4.89
	K-dJ	11.97	04.92	03.15		K-dJ	12.62	06.95	4.84
	K-NJ	11.87	04.92	03.14		K-NJ	12.77	06.95	4.84
	K-vL	05.70	04.90	03.13		K-vL	08.81	06.90	4.80
$M_{\mathrm{ZAMS}}=25M_{\odot }$					$M_{\mathrm{ZAMS}}=30M_{\odot }$
$\eta$	ID	$M\left[M_{\odot }\right]$	$M_{\mathrm{He}}\left[M_{\odot }\right]$	$M_{\mathrm{CO}}\left[M_{\odot }\right]$	$\eta$	ID	$M\left[M_{\odot }\right]$	$M_{\mathrm{He}}\left[M_{\odot }\right]$	$M_{\mathrm{CO}}\left[M_{\odot }\right]$
0.1	V-dJ	23.85	09.14	06.58	0.1	V-dJ	28.11	10.97	7.97
	V-NJ	23.86	09.14	06.59		V-NJ	28.17	10.97	7.94
	V-vL	23.56	09.14	06.59		V-vL	28.70	10.96	7.95
	K-dJ	23.74	09.23	06.66		K-dJ	28.13	10.91	7.90
	K-NJ	23.76	09.24	06.67		K-NJ	28.19	10.91	7.86
	K-vL	23.24	09.23	06.67		K-vL	28.75	10.90	7.91
0.33	V-dJ	21.98	08.87	06.37	0.33	V-dJ	24.16	11.11	8.13
	V-NJ	22.02	08.87	06.38		V-NJ	24.34	11.11	8.13
	V-vL	22.13	08.87	06.36		V-vL	26.02	11.10	8.12
	K-dJ	21.27	09.11	06.57		K-dJ	25.03	10.89	7.88
	K-NJ	21.34	09.11	06.57		K-NJ	25.18	10.89	7.87
	K-vL	20.26	09.11	06.55		K-vL	27.12	10.87	7.80
1.0	V-dJ	13.29	09.05	06.50	1.0	V-dJ	15.36	11.28	8.23
	V-NJ	13.67	09.05	06.48		V-NJ	15.78	11.29	8.26
	V-vL	12.97	09.04	06.51		V-vL	18.29	11.31	8.27
	K-dJ	15.57	08.89	06.38		K-dJ	18.51	10.92	7.96
	K-NJ	15.84	08.89	06.37		K-NJ	18.80	10.92	7.94
	K-vL	14.69	08.90	06.38		K-vL	22.53	10.89	7.94
$M_{\mathrm{ZAMS}}=35M_{\odot }$					a
$\eta$	ID	$M\left[M_{\odot }\right]$	$M_{\mathrm{He}}\left[M_{\odot }\right]$	$M_{\mathrm{CO}}\left[M_{\odot }\right]$
0.1	V-dJ-NL	32.61	10.25	07.19
	V-NJ-NL	32.69	10.24	07.19
	V-vL-NL	33.99	10.03	07.03
	K-dJ-NL	32.43	11.87	08.54
	K-NJ-NL	32.53	11.87	08.56
	K-vL-NL	33.97	11.83	08.51
0.33	V-dJ-NL	27.36	12.72	09.32
	V-NJ-NL	27.62	12.72	09.32
	V-vL-NL	31.40	12.67	09.23
	K-dJ-NL	27.50	12.68	09.21
	K-NJ-NL	27.78	12.68	09.28
	K-vL-NL	31.76	12.63	09.21
1.0	V-dJ-NL	20.03	13.30	09.81
	V-NJ-NL	20.20	13.30	09.94
	V-vL-NL	24.85	13.30	09.97
	K-dJ-NL	20.32	13.82	10.33
	K-NJ-NL	20.45	13.82	10.36
	K-vL-NL	25.13	13.82	10.41
	V-dJ-H	19.73	13.30	09.93
	V-NJ-H	19.93	13.30	09.99
	K-dJ-H	19.48	13.81	10.41
	K-NJ-H	19.69	13.82	10.41

Table 4. Table 4: Final color and surface properties of the computed models. We follow the definitions of Georgy ( 2012 ) to classify models as RSGs, YSGs, or BSGs (see also text). WR stars have X s ≤ 0.4 subscript 𝑋 𝑠 0.4 X_{s}\leq 0.4 , regardless of their surface temperature. The first and second columns indicate the wind efficiency and the mass loss algorithm combination, respectively. The WR stars are computed twice, once with the NL mass loss algorithm (see Sec.  A.6 ) and once with the H algorithm (see Sec.  A.7 ).

End of mass loss phase: $T_c\ge {10}^{9}K$
$M_{\mathrm{ZAMS}}=15M_{\odot }$						$M_{\mathrm{ZAMS}}=20M_{\odot }$
$\eta$	ID	$R\left[R_{\odot }\right]$	${\log }_{10}\left(L/L_{\odot }\right)$	${\log }_{10}\left(T_{\mathrm{eff}}/\left[K\right]\right)$	color	$\eta$	ID	$R\left[R_{\odot }\right]$	${\log }_{10}\left(L/L_{\odot }\right)$	${\log }_{10}\left(T_{\mathrm{eff}}/\left[K\right]\right)$	color
0.1	V-dJ	911	5.06	3.55	RSG	0.1	V-dJ	992	5.27	3.58	RSG
	V-NJ	911	5.06	3.55	RSG		V-NJ	992	5.27	3.58	RSG
	V-vL	927	5.06	3.54	RSG		V-vL	996	5.27	3.58	RSG
	K-dJ	914	5.06	3.55	RSG		K-dJ	987	5.27	3.58	RSG
	K-NJ	914	5.06	3.55	RSG		K-NJ	989	5.27	3.58	RSG
	K-vL	930	5.06	3.54	RSG		K-vL	991	5.27	3.58	RSG
0.33	V-dJ	916	5.05	3.54	RSG	0.33	V-dJ	998	5.27	3.58	RSG
	V-NJ	916	5.05	3.54	RSG		V-NJ	999	5.27	3.58	RSG
	V-vL	968	5.05	3.53	RSG		V-vL	985	5.27	3.58	RSG
	K-dJ	920	5.06	3.54	RSG		K-dJ	998	5.27	3.58	RSG
	K-NJ	921	5.06	3.54	RSG		K-NJ	998	5.27	3.58	RSG
	K-vL	959	5.06	3.54	RSG		K-vL	991	5.27	3.58	RSG
1.0	V-dJ	895	5.02	3.54	RSG	1.0	V-dJ	950	5.28	3.59	RSG
	V-NJ	896	5.02	3.54	RSG		V-NJ	956	5.28	3.59	RSG
	V-vL	630	5.02	3.62	YSG		V-vL	672	5.28	3.67	YSG
	K-dJ	947	5.05	3.54	RSG		K-dJ	972	5.27	3.59	RSG
	K-NJ	949	5.05	3.54	RSG		K-NJ	976	5.27	3.58	RSG
	K-vL	644	5.05	3.62	YSG		K-vL	712	5.27	3.65	YSG
$M_{\mathrm{ZAMS}}=25M_{\odot }$						$M_{\mathrm{ZAMS}}=30M_{\odot }$
$\eta$	ID	$R\left[R_{\odot }\right]$	${\log }_{10}\left(L/L_{\odot }\right)$	${\log }_{10}\left(T_{\mathrm{eff}}/\left[K\right]\right)$	color	$\eta$	ID	$R\left[R_{\odot }\right]$	${\log }_{10}\left(L/L_{\odot }\right)$	${\log }_{10}\left(T_{\mathrm{eff}}/\left[K\right]\right)$	color
0.1	V-dJ	899	5.41	3.64	YSG	0.1	V-dJ	698	5.50	3.72	YSG
	V-NJ	900	5.42	3.64	YSG		V-NJ	697	5.51	3.72	YSG
	V-vL	902	5.42	3.64	YSG		V-vL	696	5.51	3.72	YSG
	K-dJ	891	5.43	3.64	YSG		K-dJ	705	5.53	3.72	YSG
	K-NJ	891	5.44	3.65	YSG		K-NJ	704	5.53	3.72	YSG
	K-vL	893	5.43	3.64	YSG		K-vL	708	5.52	3.72	YSG
0.33	V-dJ	926	5.42	3.63	YSG	0.33	V-dJ	670	5.54	3.73	YSG
	V-NJ	924	5.41	3.63	YSG		V-NJ	670	5.54	3.73	YSG
	V-vL	932	5.42	3.63	YSG		V-vL	675	5.55	3.73	YSG
	K-dJ	896	5.42	3.64	YSG		K-dJ	719	5.53	3.72	YSG
	K-NJ	897	5.42	3.64	YSG		K-NJ	717	5.52	3.71	YSG
	K-vL	899	5.43	3.64	YSG		K-vL	711	5.53	3.72	YSG
1.0	V-dJ	797	5.43	3.67	YSG	1.0	V-dJ	614	5.56	3.76	YSG
	V-NJ	811	5.43	3.67	YSG		V-NJ	628	5.56	3.75	YSG
	V-vL	782	5.43	3.67	YSG		V-vL	610	5.55	3.76	YSG
	K-dJ	874	5.42	3.65	YSG		K-dJ	755	5.54	3.71	YSG
	K-NJ	878	5.42	3.65	YSG		K-NJ	758	5.55	3.71	YSG
	K-vL	857	5.42	3.65	YSG		K-vL	723	5.53	3.72	YSG
$M_{\mathrm{ZAMS}}=35M_{\odot }$						a
$\eta$	ID	$R\left[R_{\odot }\right]$	${\log }_{10}\left(L/L_{\odot }\right)$	${\log }_{10}\left(T_{\mathrm{eff}}/\left[K\right]\right)$	color
0.1	V-dJ	860	5.53	3.68	YSG
	V-NJ	861	5.53	3.68	YSG
	V-vL	920	5.52	3.66	YSG
	K-dJ	726	5.59	3.73	YSG
	K-NJ	730	5.59	3.73	YSG
	K-vL	755	5.58	3.72	YSG
0.33	V-dJ	529	5.62	3.80	BSG
	V-NJ	531	5.62	3.80	BSG
	V-vL	569	5.61	3.79	YSG
	K-dJ	539	5.62	3.80	BSG
	K-NJ	541	5.62	3.80	BSG
	K-vL	595	5.62	3.78	YSG
1.0	V-dJ-NL	258	5.53	3.94	WR
	V-NJ-NL	255	5.66	3.97	WR
	V-vL	398	5.66	3.88	BSG
	K-dJ-NL	167	5.75	4.09	WR
	K-NJ-NL	168	5.73	4.08	WR
	K-vL	313	5.67	3.93	BSG
	V-dJ-H	253	5.60	3.96	WR
	V-NJ-H	260	5.65	3.97	WR
	K-dJ-H	157	5.72	4.09	WR
	K-NJ-H	158	5.72	4.09	WR

Table 5. Table 5: Maximum spread in total mass, core masses, and radius at the end of the mass loss phase, for models differing in mass loss algorithm combination and efficiency η 𝜂 \eta . We also list the maximum and minimum total mass for each M ZAMS subscript 𝑀 ZAMS M_{\mathrm{ZAMS}} .

$M_{\mathrm{ZAMS}}$	$\max \Delta M$	$\max M$	$\min M$	$\max \Delta M_{\mathrm{He}}$	$\max \Delta M_{\mathrm{CO}}$	$\max \Delta R$
$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	$\left[R_{\odot }\right]$
15	09.42	14.66	05.25	0.36	0.28	338
20	10.58	19.38	08.80	0.16	0.12	327
25	10.89	23.86	12.97	0.37	0.31	150
30	13.39	28.75	15.36	0.44	0.47	148
35	14.51	33.99	19.48	3.79	3.38	763

Table 6. Table 6: Stellar properties at the end of hot phase, when T eff subscript 𝑇 eff T_{\mathrm{eff}} decreases below 15 000 ​ K 15000 K 15\,000\,\mathrm{K} for the first time. For the descriptions of the V and K algorithms see Sec.  A.1 and Sec.  A.2 , respectively.

Hot wind algorithm	$M_{\mathrm{ZAMS}}\left[M_{\odot }\right]$	$\eta$	$R\left[R_{\odot }\right]$	$L\left[{10}^4{L}_{\odot }\right]$	$M\left[M_{\odot }\right]$	$M_{\mathrm{He}}\left[M_{\odot }\right]$	age [Myr]
V	15	1.00	034.85	05.44	14.54	03.69	12.7917
K	15	1.00	035.44	05.64	14.73	03.71	12.8239
V	15	0.33	035.91	05.76	14.85	03.76	12.7503
K	15	0.33	035.97	05.82	14.91	03.77	12.7612
V	15	0.10	036.27	05.85	14.95	03.78	12.7363
K	15	0.10	035.96	05.83	14.97	03.78	12.7359
V	20	1.00	051.97	12.20	19.06	05.83	08.8937
K	20	1.00	052.75	12.32	19.42	05.87	08.8947
V	20	0.33	053.10	12.67	19.69	05.95	08.8373
K	20	0.33	053.35	12.76	19.81	05.95	08.8374
V	20	0.10	053.46	12.80	19.91	05.99	08.8176
K	20	0.10	054.39	13.10	19.94	05.98	08.8183
V	25	1.00	067.58	20.63	23.24	08.03	07.0367
K	25	1.00	068.70	21.41	24.00	07.91	07.0094
V	25	0.33	070.59	22.55	24.44	07.97	06.9653
K	25	0.33	070.05	22.04	24.66	08.07	06.9546
V	25	0.10	070.88	22.39	24.84	08.10	06.9390
K	25	0.10	070.50	22.19	24.90	08.13	06.9355
V	30	1.00	080.82	29.52	27.00	10.18	05.9700
K	30	1.00	084.32	31.82	28.46	10.00	05.9215
V	30	0.33	084.64	32.11	29.10	10.16	05.8853
K	30	0.33	085.14	32.83	29.48	10.11	05.8689
V	30	0.10	085.58	32.85	29.77	09.96	05.8558
K	30	0.10	085.34	32.93	29.84	09.89	05.8503
V	35	1.00	095.27	40.28	30.32	12.29	05.2807
K	35	1.00	097.93	41.83	32.83	12.51	05.2164
V	35	0.33	098.15	43.52	33.62	12.19	05.1879
K	35	0.33	099.84	44.24	34.27	12.09	05.1661
V	35	0.10	100.81	45.59	34.63	11.58	05.1559
K	35	0.10	100.35	45.03	34.78	11.78	05.1484

Table 7. Table 7: Stellar properties at core oxygen depletion, i.e.  X c ​ ( O 16 ) < 0.04 subscript 𝑋 𝑐 superscript O 16 0.04 X_{c}(\mathrm{{}^{16}O})<0.04 and X c ​ ( Si 28 ) > 0.01 subscript 𝑋 𝑐 superscript Si 28 0.01 X_{c}(\mathrm{{}^{28}Si})>0.01 . The last column shows the maximum difference in the compactness parameter ξ 2.5 O ​ depl superscript subscript 𝜉 2.5 O depl \xi_{2.5}^{\mathrm{O\ depl}} for each choice of η 𝜂 \eta and M ZAMS subscript 𝑀 ZAMS M_{\mathrm{ZAMS}} . These runs are re-started from the corresponding MESA models saved at the end of the mass loss phase, when T c ≥ 10 9 ​ K subscript 𝑇 𝑐 superscript 10 9 K T_{c}\geq 10^{9}\,\mathrm{K} .

Oxygen depletion: and $X_c\left({}_{}{}^{28}\mathrm{Si}\right)>0.01$
$M_{\mathrm{ZAMS}}\left[M_{\odot }\right]$	$\eta$	ID	$R\left[R_{\odot }\right]$	$M_{\mathrm{tot}}\left[M_{\odot }\right]$	$M_{\mathrm{He}}\left[M_{\odot }\right]$	$M_{\mathrm{CO}}\left[M_{\odot }\right]$	${\xi }_{2.5}^{O\mathrm{depl}}$	$\max \Delta {\xi }_{2.5}^{O\mathrm{depl}}$
15	0.13	V-dJ	0908	14.66	04.99	3.17	0.155	0.001
		V-vL	0924	13.60	04.98	3.17	0.156
		K-NJ	0911	14.66	05.00	3.17	0.156
	0.33	V-NJ	0914	13.90	04.93	3.13	0.153	0.001
		V-vL	0967	10.39	04.93	3.13	0.152
	1.03	V-NJ	0895	12.74	04.65	2.95	0.141	0.014
		V-vL	0629	05.25	04.64	2.94	0.139
		K-NJ	0946	11.87	04.92	3.12	0.153
		K-vL	0643	05.70	04.90	3.11	0.152
20	0.13	V-dJ	0994	19.23	07.04	4.83	0.182	0.007
		V-vL	0998	18.10	07.03	4.83	0.178
		K-NJ	0991	19.37	07.01	4.82	0.179
		K-vL	0991	18.67	07.01	4.82	0.175
	0.33	V-NJ	1001	17.48	07.01	4.81	0.161	0.052
		V-vL	0987	13.47	06.99	4.80	0.175
		K-dJ	0999	17.62	07.00	4.80	0.213
		K-vL	0993	13.90	06.98	4.80	0.175
	1.03	V-dJ	0951	11.81	07.06	4.85	0.182	0.010
		V-vL	0673	08.80	07.02	4.82	0.176
		K-NJ	0978	12.77	06.95	4.77	0.176
		K-vL	0712	08.81	06.90	4.73	0.172
25	0.13	V-dJ	0898	23.85	09.14	6.46	0.180	0.023
		V-vL	0899	23.56	09.14	6.44	0.179
		K-NJ	0888	23.76	09.24	6.52	0.159
		K-vL	0889	23.24	09.23	6.50	0.157
	0.33	V-NJ	0922	22.02	08.87	6.25	0.161	0.049
		V-vL	0929	22.13	08.87	6.20	0.164
		K-dJ	0894	21.27	09.11	6.42	0.210
	1.03	V-NJ	0812	13.67	09.05	6.32	0.164	0.050
		V-vL	0786	12.97	09.04	6.29	0.211
		K-dJ	0875	15.57	08.89	6.21	0.161
		K-vL	0860	14.69	08.90	6.23	0.200
30	0.13	V-dJ	0701	28.11	10.97	7.89	0.242	0.001
		V-NJ	0697	28.17	10.97	7.90	0.243
		K-NJ	0706	28.19	10.91	7.83	0.242
		K-vL	0719	28.75	10.90	7.74	0.242
	0.33	V-dJ	0668	24.16	11.11	7.99	0.243	0.023
		V-vL	0723	26.02	11.10	8.00	0.243
		K-NJ	0716	25.18	10.89	7.85	0.228
		K-vL	0709	27.12	10.87	7.79	0.220
	1.03	V-dJ	0610	15.36	11.28	8.20	0.179	0.065
		K-dJ	0756	18.51	10.92	7.85	0.244
		K-vL	0723	22.53	10.89	7.74	0.243

Table 8. Table 8: Properties of the subset of models run to the onset of core collapse, defined as the time when max ​ { | v | } ≥ 10 3 ​ km ​ s − 1 max 𝑣 superscript 10 3 km superscript s 1 \mathrm{max}\{|v|\}\geq 10^{3}\ \mathrm{km\ s^{-1}} . M 4 subscript 𝑀 4 M_{4} and μ 4 subscript 𝜇 4 \mu_{4} are the parameters used to predict the SN outcome of a stellar model in Ertl et al. ( 2016 ) , see also text. M ρ 6 = def M ​ ( ρ = 10 6 ​ g ​ cm − 3 ) superscript def subscript 𝑀 subscript 𝜌 6 𝑀 𝜌 superscript 10 6 g superscript cm 3 M_{\rho_{6}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}M(\rho=10^{6}\mathrm{g\ cm^{-3})} is the mass enclosed in the location where the density drops below 10 6 ​ g ​ cm − 3 superscript 10 6 g superscript cm 3 10^{6}\ \mathrm{g\ cm^{-3}} , M CO subscript 𝑀 CO M_{\mathrm{CO}} and M Fe subscript 𝑀 Fe M_{\mathrm{Fe}} are the carbon-oxygen and iron core masses, respectively.

Onset of core collapse: $\max \left\{\left|v\right|\right\}\ge {10}^3\mathrm{km}{s}^{-1}$
$M_{\mathrm{ZAMS}}\left[M_{\odot }\right]$	$\eta$	ID	${\xi }_{2.5}^{\mathrm{pre}-\mathrm{SN}}$	$M_4\left[M_{\odot }\right]$	${\mu }_4$	$M_{{\rho }_6}\left[M_{\odot }\right]$	$M_{\mathrm{CO}}\left[M_{\odot }\right]$	$M_{\mathrm{Fe}}\left[M_{\odot }\right]$
15	1.0	V-NJ	0.103	1.71	0.045	1.68	2.91	1.39
		K-vL	0.132	1.78	0.051	1.79	3.07	1.50
25	0.33	V-vL	0.227	1.73	0.084	1.84	6.38	1.51
		K-dJ	0.308	2.05	0.100	2.19	6.40	1.63
30	0.33	V-dJ	0.358	1.60	0.163	2.21	7.98	1.56
		K-NJ	0.276	1.82	0.100	1.98	7.90	1.58

Table 9. Table 9: Functional dependence of the mass loss rate M ˙ ˙ 𝑀 \dot{M} on stellar and wind parameters for each algorithm combination and mass loss phase. L 𝐿 L is the luminosity, M 𝑀 M is the mass, v esc subscript 𝑣 esc v_{\mathrm{esc}} is the escape velocity, v ∞ subscript 𝑣 v_{\infty} is the final velocity of the wind, v th subscript 𝑣 th v_{\mathrm{th}} is the thermal velocity, T eff subscript 𝑇 eff T_{\mathrm{eff}} is the effective temperature, R 𝑅 R the radius, Γ E = L / L Edd subscript Γ E 𝐿 subscript 𝐿 Edd \Gamma_{\mathrm{E}}=L/L_{\mathrm{Edd}} the Eddington ratio, and X s subscript 𝑋 𝑠 X_{s} , Y s subscript 𝑌 𝑠 Y_{s} , and Z s subscript 𝑍 𝑠 Z_{s} are the surface hydrogen abundance, helium abundance, and metallicity, respectively. The scalings are obtained from the algorithms described in Appendix  A , and errors and overall multiplying factors are omitted for the sake of brevity. See Tab.  2 for the naming convention for mass loss algorithm combinations.

ID	V-dJ-NL	V-NJ-NL	V-vL-NL	V-dJ-H	V-NJ-H	V-vL-H
Hot	$L^{2.210}{M}^{-1.339}{\left(\frac{v_{\infty }}{2v_{\mathrm{esc}}}\right)}^{-1.601}{T}_{\mathrm{eff}}^{1.07}{Z}^{0.85}$aaaaaaaaaaaaaaaa if
	$L^{2.194}{M}^{-1.313}{\left(\frac{v_{\infty }}{2v_{\mathrm{esc}}}\right)}^{-1.226}{Z}^{0.85}{T}_{\mathrm{eff}}^{\left(0.933-10.92{\log }_{10}\left(T_{\mathrm{eff}}/40kK\right)\right)}$ if $T_{\mathrm{eff}}>27.5\mathrm{kK}$
Cool	$L^{1.769}{T}_{\mathrm{eff}}^{-1.676}$	$L^{1.24}{M}^{0.16}{R}^{0.81}$	$L^{1.05}{T}_{\mathrm{eff}}^{-6.3}$	$L^{1.769}{T}_{\mathrm{eff}}^{-1.676}$	$L^{1.24}{M}^{0.16}{R}^{0.81}$	$L^{1.05}{T}_{\mathrm{eff}}^{-6.3}$
WR	$L^{1.29}{Y}_s^{1.73}{Z}_s^{0.47}$			$L^{1.5}{10}^{-2.85X_s}$ if ${\log }_{10}\left(L/L_{\odot }\right)>4.5$
				$L^{6.8}$aaaaaaaa if ${\log }_{10}\left(L/L_{\odot }\right)\le 4.5$
ID	K-dJ-NL	K-NJ-NL	K-vL-NL	K-dJ-H	K-NJ-H	K-vL-H
Hot	$L^{1.779}{v}_{\mathrm{th}}^{-1.169}{\left[GM\left(1-{\Gamma }_E\right)\right]}^{0.610}$
Cool	$L^{1.769}{T}_{\mathrm{eff}}^{-1.676}$	$L^{1.24}{M}^{0.16}{R}^{0.81}$	$L^{1.05}{T}_{\mathrm{eff}}^{-6.3}$	$L^{1.769}{T}_{\mathrm{eff}}^{-1.676}$	$L^{1.24}{M}^{0.16}{R}^{0.81}$	$L^{1.05}{T}_{\mathrm{eff}}^{-6.3}$
WR	$L^{1.29}{Y}_s^{1.73}{Z}_s^{0.47}$			$L^{1.5}{10}^{-2.85X_s}$ if ${\log }_{10}\left(L/L_{\odot }\right)>4.5$
				$L^{6.8}$aaaaaaaa if ${\log }_{10}\left(L/L_{\odot }\right)\le 4.5$

Equations

$$\dot{M}\propto Z^a,$$

$$\max \{|v|\}\ge 10^3\mathrm{km}{\mathrm{s}}^{-1},$$

$${\xi }_{\mathcal{M}}\stackrel{\mathrm{def}}{=}\frac{\mathcal{M}/M_{\odot }}{R(\mathcal{M})/1000\mathrm{km}}.$$

$${\log }_{10}(-\dot{M})=-6.697(61)+2.194(21){\log }_{10}\left(\frac{L}{10^5{L}_{\odot }}\right)+$$

$$-1.313(46){\log }_{10}\left(\frac{M}{30M_{\odot }}\right)-1.226(37){\log }_{10}\left(\frac{v_{\infty }}{2v_{\mathrm{esc}}}\right)+$$

$$+0.933(64){\log }_{10}\left(\frac{T_{\mathrm{eff}}}{40000\mathrm{K}}\right)-10.92(90){\log }_{10}^2\left(\frac{T_{\mathrm{eff}}}{40000\mathrm{K}}\right)+$$

$$+0.85(10){\log }_{10}\left(\frac{Z}{Z_{\odot }}\right),$$

$${\log }_{10}(-\dot{M})=-6.668(80)+2.210(31){\log }_{10}\left(\frac{L}{10^5{L}_{\odot }}\right)+$$

$$-1.339(68){\log }_{10}\left(\frac{M}{30M_{\odot }}\right)-1.601(55){\log }_{10}\left(\frac{v_{\infty }}{2v_{\mathrm{esc}}}\right)+$$

$$+1.07(10){\log }_{10}\left(\frac{T_{\mathrm{eff}}}{40000\mathrm{K}}\right)+0.85(10){\log }_{10}\left(\frac{Z}{Z_{\odot }}\right),$$

$$g_{\mathrm{ph}}=g_{\mathrm{rad}}^{\mathrm{Th}}\left(1+k{\left(\frac{{\sigma }_{\mathrm{Th}}{c}_s}{dv/dr}\right)}^{-\alpha }{\left(\frac{2n_e}{\sqrt{1-(1-(R/r{)}^2}}\right)}^{\delta }\mathrm{CF}\left(r,v,\frac{dv}{dr},\alpha \right)\right),$$

$$\mathrm{CF}\left(r,v,\frac{dv}{dr},\alpha\right)=\frac{1}{\alpha+1}\frac{x^{2}}{1-h}\left(1-\left(1-\frac{1}{x^{2}}+\frac{h}{x^{2}}\right)^{\alpha+1}\right)\ \ .\$$

$$v(r)=v_{\infty }{\left(1-\frac{R}{r}\right)}^{\beta },$$

$$\dot{M}\equiv \dot{M}(\alpha ,\delta ,k,M,L,v_{th})=\tilde{D}(\alpha ,\delta ,v){\left(\frac{{\sigma }_{\mathrm{Th}}kL}{4\pi c}\right)}^{1/(\alpha -\delta )}\times$$

$$\times {\left(\frac{4\pi \alpha }{{\sigma }_{\mathrm{Th}}{v}_{\mathrm{th}}}\right)}^{\alpha /(\alpha -\delta )}{\left(\frac{1-\alpha }{GM(1-\Gamma )}\right)}^{(1-\alpha )/(\alpha -\delta )},$$

$${\log }_{10}(-\dot{M})=\sum _{n=0}^N\sum _{\stackrel{i=0}{j=n-i}}^{i=n}{a}_{ij}{T}_i({\log }_{10}(T_{\mathrm{eff}}))\cdot T_j({\log }_{10}(L)).$$

$${\log }_{10}(-\dot{M})=1.769{\log }_{10}(L/L_{\odot })-1.676{\log }_{10}(T_{\mathrm{eff}}/[\mathrm{K}])-8.158.$$

$$t^{(d)}\stackrel{\mathrm{def}}{=}\frac{\delta t}{\sqrt{[\delta {\log }_{10}(T_{\mathrm{eff}}/[\mathrm{K}]){]}^2+[\delta {\log }_{10}(L/L_{\odot }){]}^2}},$$

$$\bar{M}=\frac{{\sum }_{n=1}^N\Psi (M_n)\frac{d{M}_n}{d{\log }_{10}(L/L_{\odot })}{t}_n^{(d)}{M}_n}{{\sum }_{n=1}^N\Psi (M_n)\frac{d{M}_n}{d{\log }_{10}(L/L_{\odot })}{t}_n^{(d)}},$$

$${\log }_{10}(-\dot{M})=-14.02+1.24{\log }_{10}(L/L_{\odot })+$$

$$+0.16{\log }_{10}(M/M_{\odot })+0.81{\log }_{10}(R/R_{\odot }).$$

$${\log }_{10}(-\dot{M})=-5.65(15)+1.05(14){\log }_{10}(L/10^4{L}_{\odot })+$$

$$-6.3(1.2){\log }_{10}(T_{\mathrm{eff}}/3500\mathrm{K}),$$

$${\log }_{10}(-\dot{M})=-11.0+1.29(14){\log }_{10}(L/L_{\odot })+$$

$$1.73(42){\log }_{10}(Y)+0.47(09){\log }_{10}(Z),$$

$$\log_{10}(-\dot{M})=\begin{cases}-12.25+1.5\log_{10}\left(\frac{L}{L_{\odot}}\right)-2.85X_{s}&\mbox{if }\mbox{$L>4.5L_{\odot}$},\\ -35.8+6.8\log_{10}\left(\frac{L}{L_{\odot}}\right)&\mbox{if }\mbox{$L\leq 4.5L_{\odot}$}.\end{cases}$$

$$T^4=\frac{3}{4}{T}_{\mathrm{eff}}^4\left(\tau +\frac{2}{3}\right),$$

$$\Delta t_{n+1}\le \min \{t_{\mathrm{KH}},t_{\dot{M}}\},$$

$$t_{\mathrm{KH}}\stackrel{\mathrm{def}}{=}\frac{3}{4}\frac{G{M}^2}{RL},t_{\dot{M}}\stackrel{\mathrm{def}}{=}\frac{|\dot{M}|}{M},$$

$$\mathcal{L}\sim \left((N_{\mathrm{iso}}+5)\cdot N_{\mathrm{z}}\right)\cdot (3N_{\mathrm{iso}}+9),$$

Full text

††offprints: m.renzo@uva.nl

1 Introduction

Mass loss is a key phenomenon for the co-evolution of massive stars ($M\gtrsim 8\,M_{\odot}$) and their environment, yet it is poorly understood. It plays an important role throughout the stellar evolution and it may have a deciding influence on the outcome of core collapse. Mass loss is responsible for a large part of the chemical enrichment of the interstellar medium, and its momentum input can trigger star formation, but it can also sweep away gas from stellar clusters, preventing further star formation.

In the standard picture of single massive star evolution, mass loss influences the duration of different evolutionary phases (e.g., Meynet et al., 2015), especially the amount of time spent on the red supergiant (RSG) branch. It has been suggested to be important for the solution of the so-called “red supergiant problem” (Smartt et al., 2009; Smith et al., 2011), i.e. the discrepancy between the observed maximum mass for type IIP supernovae (SNe) and the theoretical predictions for the core collapse of RSG stars. This discrepancy indicates an incomplete theoretical understanding of the evolution and explosion of massive stars, especially in the mass range \sim$$16-30\,M_{\odot}.

Mass loss also plays an essential role in the formation of Wolf-Rayet (WR) stars. Two competing scenarios exist for the removal of (most of) the hydrogen-rich envelope of stars: the single star picture (so-called “Conti scenario”, Conti, 1975; Maeder & Conti, 1994; Lamers, 2013) and the binary formation channel (e.g., Woosley et al., 1995; Wellstein & Langer, 1999; Smith & Tombleson, 2015; Shara et al., 2017). Understanding mass loss phenomena is necessary to discriminate between these two scenarios. Mass loss is invoked to explain the variety of core-collapse SNe (e.g., Eldridge & Tout, 2004; Smith et al., 2011; Georgy, 2012; Groh et al., 2013; Smith, 2014), because it can modify the surface composition of the pre-SN star, possibly removing the hydrogen-rich envelope (and perhaps some of the helium-rich shell) and leading to type Ib/c SNe. Observations of SNe IIn and superluminous SNe may perhaps be explained by the strong interaction between the SN shock and shells of material ejected from the star in late stages of its evolution, (e.g., Smith et al., 2011; Smith, 2014; Shiode & Quataert, 2014). It has also been proposed that the coupling of mass loss and rotation might prevent pair instability SNe for very massive, low metallicity stars (Ekström et al., 2008; Woosley, 2016).

Finally, mass loss plays a key role in shaping the (hydrostatic) internal structure of massive stars at the pre-SN stage, which can influence the expected SN outcome (see, e.g., Belczynski et al., 2010; O’Connor & Ott, 2011, 2013; Ugliano et al., 2012; Sukhbold et al., 2016): will the star successfully explode and leave a neutron star (NS) remnant? Will the explosion be highly asymmetric or weak, leading to fallback accretion and black hole (BH) formation? Or will the explosion fail completely, leaving a BH with little (Lovegrove & Woosley, 2013) or no electromagnetic counterpart?

Depending on the amount and geometry of the ejecta, which are governed by the mass loss during the stellar lifetime and the SN energetics, also the SN kick can vary (Zwicky, 1957; Blaauw, 1961; Boersma, 1961; Janka, 2013, 2016). The kick can change the post-explosion orbital parameters if the star is in a binary system and, therefore, our incomplete understanding of massive star mass loss affects the predicted populations of sources for gravitational wave astronomy (Abbott et al., 2016).

There exist three main channels of mass loss in the evolution of massive stars:

Steady-state winds. These are radiatively driven in hot stars (i.e., on the main sequence; see Puls et al. 2008 for a review). In the supergiant phase, the driving mechanism is uncertain. Winds could be radiatively driven via lines or dust, or by other mechanisms, e.g., wave energy deposition (see, e.g., Bennett 2010 for a brief review). 2. 2.

Impulsive, pulsational and/or eruptive mass loss, e.g. disk shedding at criticial rotation or giant eruptions such as those of Luminous Blue Variables (LBVs) (e.g., Puls et al., 2008; Smith, 2014; Puls et al., 2015); 3. 3.

Roche lobe overflow (RLOF) and possibly common envelope ejection in binary systems, which can result in mass loss from the donor star and also mass loss from the system as a whole in non-conservative cases.

Which of these processes dominates in terms of the total mass shed has been a matter of debate in the literature (see Smith, 2014, and references therein).

Mass loss is an intrinsically dynamical phenomenon which involves bulk acceleration of matter for escaping the star. Because of the dynamical nature of mass loss, it is difficult to include in stellar evolution codes: most simulations focus on single stars, and are carried out with hydrostatic codes which cannot account for dynamical or impulsive events in a physical way. Even hydrodynamical codes do not permit a self-consistent development of impulsive events, since the physical processes triggering them are currently poorly understood or even unknown. Impulsive outbursts of mass loss can, however, be included using physically plausible prescriptions (e.g., Morozova et al., 2015).

Most massive stars are found in binary systems (e.g., Mason et al., 2009; Sana & Evans, 2011; Sana et al., 2012; Kiminki & Kobulnicky, 2012; Chini et al., 2012; Kobulnicky et al., 2014; Almeida et al., 2016), where mass loss also determines the angular momentum losses, thus the orbital evolution, and ultimately the binary evolution path and its end point (merger, disruption of the binary system, double compact object binary, etc.).

In this paper, our focus is on steady, radiatively-driven wind mass loss (Lucy & Solomon, 1970) and our goal is to understand how different treatments of this process affect massive star evolution and pre-SN structure. We consider single massive stars and do not address the problem of mass loss in binaries directly.

Most evolutionary calculations of massive stars only include wind mass loss, which can be treated in the steady state approximation, although, strictly speaking, a wind is dynamical as well.

The wind is in fact driven by a radiative acceleration which formally enters into the momentum equation for the stellar plasma (Castor et al., 1975), the stellar structure responds secularly, since wind mass loss rates are low and change only slowly compared to impulsive mass loss events. At solar metallicity, the wind mass loss rate is $10^{-10}\,M_{\odot}\,\mathrm{yr^{-1}}\lesssim\dot{M}_{\mathrm{wind}}\lesssim 10^{-5}\,M_{\odot}\,\mathrm{yr^{-1}}$, while RLOF and LBV eruptions yield $10^{-6}\,M_{\odot}\,\mathrm{yr^{-1}}\lesssim\dot{M}_{\mathrm{RLOF}}\lesssim\dot{M}_{LBV}\lesssim 10^{-2}\,M_{\odot}\,\mathrm{yr^{-1}}$, (e.g., de Jager et al., 1988; Vink et al., 2001; van Loon et al., 2005; Smith & Owocki, 2006; Crowther, 2007; Puls et al., 2008; Langer, 2012; Smith, 2014). These mass loss rates correspond to timescales $\tau_{\mathrm{wind}}=M/\dot{M}_{\mathrm{wind}}\gg\tau_{\mathrm{RLOF}},\tau_{\mathrm{LBV}}$. Therefore, wind mass loss can be included in stellar evolution calculations using parametric algorithms111These are often called “recipes” in jargon. However, we prefer the term “algorithm” (Fibonacci, 1202) because it underlines that these are mathematical representations of physical phenomena relying on specific sets of assumptions..

Stellar winds are radiatively driven by the interaction of photons with metallic ions (line-driven mass loss) or dust grains (dust-driven mass loss). Therefore, they depend on the opacity and thus chemical composition, ionization state, and density stratification, so indirectly also on the equation of state of the outermost stellar layers. Metals effectively provide all the opacity in stellar atmospheres because of their large number of lines. Photons come out of the photosphere with well defined direction, interact with metallic atoms/ions via bound-bound processes (absorption and line scattering) and cede their momentum to the atoms/ions. When these de-excite, they emit photons isotropically. The momentum of de-excitation photons averages out, and the result is a net gain of momentum in the direction of the initial photon (see, e.g., Puls et al., 2008).

Note that the radiation field in the stellar atmosphere, i.e. above the photosphere, is not isotropic. If there were not a net radial flux of photons, also the momentum of the incoming photons would average out. In a nutshell, the incident photons push metals outward, and metals drag hydrogen and helium through collisional Coulomb coupling (see Puls et al., 2008; Vink, 2015, and references therein).

This simple theoretical picture of line-driven stellar winds is complicated by two phenomena: the high nonlinearity of the driving mechanism, and the possible presence of inhomogeneities (so called “clumpiness”) in the stellar atmosphere. The high nonlinearity arises because the outflow of mass is driven by the rate of interaction of photons with metals, but this in turn depends on the local opacity, and therefore on the outflow properties, such as density, and velocity (which can Doppler shift the lines), see Lamers & Cassinelli (1999); Puls et al. (2008) and references therein. The presence of inhomogeneities in the stellar atmosphere is both theoretically expected (see, e.g., Owocki & Rybicki, 1984; Owocki et al., 1988; Feldmeier, 1995; Owocki & Puls, 1999; Dessart & Owocki, 2005; Puls et al., 2008; Smith, 2014) and observed comparing diagnostic spectral lines sensitive to the density $\rho$ (e.g., lines with P Cygni profiles) and $\rho^{2}$ (e.g., recombination lines, such as H$\alpha$) in the same stellar wind (see e.g. Fullerton et al., 2006; Bouret et al., 2005; Evans et al., 2004). This comparison shows that the averaged $\langle\rho^{2}\rangle>\langle\rho\rangle^{2}$ (for reviews see Puls et al. 2008 and Smith 2014). Therefore, the presence of over-dense clumps in the wind causes an overestimation of the density inferred from observed spectral features sensitive to $\rho^{2}$. This is not taken into account in most wind mass loss algorithms in the literature. The overestimation of the density directly results in an overestimated mass loss rate. Work by Crowther et al. (2002); Hillier et al. (2003); Bouret et al. (2005); Fullerton et al. (2006); Puls et al. (2008); Smith (2014); Abbott et al. (2016), suggests that the algorithms used in stellar evolution calculations may yield mass loss rates that are a factor of 2 to 10 too high. Puls et al. (2008) and Smith (2014) suggest a factor of 3 as the most realistic overestimate. We refer the interested reader to Puls et al. (2008), Smith (2014), and Renzo (2015) for more details.

To date, there has been no systematic comparison of the various wind mass loss algorithms, with varying corrections for clumpiness, and their combinations and effects throughout the evolution and on the final structure of massive stars. However, Eldridge & Tout (2004) compared different combinations of semi-empirical mass loss rates to find the threshold in mass between type Ibc and type II SNe. Yoon et al. (2010) discussed consequences of the revision of mass loss rates because of the clumpiness during the WR phase.

In this study, we employ the open-source stellar evolution code MESA (Paxton et al., 2011, 2013, 2015) to compare various combinations of mass loss algorithms, using different efficiency factors to account for the inhomogeneities in the wind (albeit in an ad-hoc fashion). Our aim is to understand the systematics of massive star evolution and pre-SN structure caused by variations in the treatment of wind mass loss, focusing on the differences in the evolution and pre-SN structure (effective temperature, total mass, core masses, and interior structure).

The remainder of this paper is structured as follows. In Sec. 2, we discuss some general aspects of the implementation of wind mass loss in stellar evolution codes and give a very brief overview of the physical bases of the wind mass loss algorithms we compare. A longer review of these, including the limitations and the formulae implemented in MESA, can be found in Appendix A. The more technical points not relevant to the physics of stellar winds are discussed in Appendix B, with the explicit aim of making our result reproducible. We compare our models when their final pre-SN mass is determined in Sec. 3.1. We discuss the impact of winds during the hot, cool, and WR phases (if reached) separately in Sec. 3.2, 3.3, and 3.4, respectively. We compare a subset of our models at oxygen depletion in Sec. 3.5. The evolution and pre-SN structure of the core is also sensitive to the mass loss history of the stellar model, including the early mass loss during the main sequence, as we show in Sec. 3.6. In Sec. 3.7, we discuss the evolution of a subset of our models from oxygen depletion to the onset of core collapse. We discuss the implications of uncertainties in wind mass loss and potential observational constraints in Sec. 4, before concluding in Sec. 5.

2 Methods

2.1 Overview of the mass loss algorithms

Stellar evolution codes do not explicitly compute the acceleration of the gas unbound in mass loss processes. The usual approach for including wind mass loss is to use parametric algorithms prescribing a mass loss rate averaged over each timestep. The time averaging is needed to compute each timestep with a constant mass loss rate. Homogeneity of the wind is implicit in the standard formalism, which is known to be a poor approximation and could cause a significant overestimate of the mass loss rate. Stellar wind algorithms are either parametric fits to observed mass loss rates, or theoretically derived models with free parameters chosen empirically or heuristically. Each algorithm gives a formula for the mass loss rate as a function of some quantities characterizing the star $\dot{M}\equiv\dot{M}(L,T_{\mathrm{eff}},Z,\dots)$. The precise set of variables assumed to be independent varies between the algorithms. Most algorithms do not include an explicit metallicity dependence, because they either assume a specific chemical composition of the stellar atmosphere or are based on observed samples with a specific metallicity $Z$. It is common practice to impose a smooth scaling with $Z$ such as

[TABLE]

with $a\simeq 0.5$ (e.g., Vink et al., 2000; Woosley et al., 2002). Eq. 1 is in reasonable agreement with more sophisticated mass loss rate determinations (see, e.g., Vink et al., 2001), but deviations should be expected both at very low $Z$ (because of the lack of metal lines to drive the wind) and very high $Z$ (because of line saturation preventing further driving of the wind). In this study, we only consider solar metallicity. Note also that stellar evolution codes usually neglect the errors on the coefficients of the algorithms obtained as fits to observations.

Since mass loss rates have large uncertainties, it is common practice to employ a linear efficiency factor $\eta$ to rescale rates to account for various physical uncertainties. For example, $\eta\lesssim 1$ can be used to account for wind clumpiness (e.g., Hamann & Koesterke, 1998; Woosley et al., 2002; Woosley & Heger, 2007). Mass loss channels other than winds exist (e.g., LBV eruptions, binary interactions, etc.) and some authors (e.g., Dessart et al., 2013; Meynet et al., 2015) use $\eta>1$ to explore the effects of increased total mass loss. However, an averaged steady wind approximation might not properly capture the readjustment of the stellar structure to non-wind mass loss events, which may be sensitive to mass loss timing, and/or the physical process triggering them, and may not be radiatively driven.

Most wind mass loss algorithms are tailored to a specific evolutionary stage. To carry out a simulation of the entire evolution of the star, several mass loss algorithms are commonly combined using computational definitions of the evolutionary phases. This may introduce somewhat arbitrary switching points in the evolution. Below, we list the physical basis and the abbreviations for the two hot phase, three cool phase, and two WR phase wind mass loss algorithms that we combine and compare here. We define each phase of the evolution in Sec. 2.2.

Vink et al. (2000, 2001)**** (V): This wind scheme is a theoretical algorithm obtained with numerical simulation of the line-driven process. It explicitly includes the metallicity dependence and applies to OB stars during their hot evolutionary phase. It also includes a detailed treatment of the so-called “bistability jump”. This corresponds to a non-monotonic behavior of the mass loss rate $\dot{M}$ as a function of effective temperature $T_{\mathrm{eff}}$ in certain temperature ranges (e.g., $T_{\mathrm{eff}}^{\mathrm{jump}}\simeq 25\,000\,\mathrm{K}$) because of the recombination of certain ions, which provides more lines in the spectral domain relevant to drive the wind, cf. Appendix A.1.

Kudritzki et al. (1989)**** (K): This analytical mass loss rate is obtained using the Castor et al. (1975) model for line-driven acceleration. It assumes that the wind is stationary, isothermal, spherically symmetric, and without viscosity and heat conduction. The analytical solution is obtained assuming a velocity structure $v\equiv v(r)$ as a function of the radius of the wind, and solving self-consistently for density and radiative acceleration, cf. Appendix A.2.

de Jager et al. (1988)**** (dJ): This empirical mass loss rate describes the “averaged statistical behavior” of stars (excluding WR and Be stars) in the HR diagram. It is commonly used for the cool (giant) phase of the evolution of massive stars, cf. also Appendix A.3.

Nieuwenhuijzen & de Jager (1990)**** (NJ): This algorithm is also an empirical mass loss rate drawn from the same sample of stars used by de Jager et al. (1988). The two algorithms differ in the physical quantities the mass loss rate is assumed to depend on: NJ used pre-computed stellar models to add a dependence on the total mass ($\dot{M}\equiv\dot{M}(M)$), which is not a directly observable quantity for single stars. It is also usually adopted for the cool phase of stellar evolution. See Appendix A.4 for more details.

van Loon et al. (2005)**** (vL): This empirical mass loss rate is derived from a sample of oxygen rich asymptotic giant branch (AGB) and red supergiants (RSG) stars in the Large Magellanic Cloud. It assumes a dust-driven wind, i.e. mass loss is driven by photons impinging on dust grains instead of metallic ions. Note, however, that the presence of dust and its role as a wind driving agent in supergiant stars is still debated in the literature (van Loon et al., 2005; Ferrarotti & Gail, 2006), cf. Appendix A.5.

Nugis & Lamers (2000)**** (NL): This empirical mass loss rate is for WR stars. Fitting the data for two populations of WR stars (one of known distance and one for which they carry out a distance determination), they provide an algorithm, which depends strongly on the surface chemical composition of the star, cf. Appendix A.6.

Hamann et al. (1982, 1995)**** (H): This is a theoretical mass loss rate for WR stars. It is derived assuming a spherically symmetric, homogeneous, and stationary wind. They avoid solving for the dynamics with a complicated radiative acceleration term by imposing a velocity structure $v\equiv v(r)$. In this way, they are able to produce synthetic spectra, which they then fit to observed WR stars to infer $\dot{M}$. Hamann & Koesterke (1998) suggest to reduce the mass loss rate by a factor between 2 and 3 to account for wind clumpiness, cf. Appendix A.7.

We report in Tab. 1 an approximate scaling of the mass loss rate with the luminosity ($L$), which quantifies the amount of photons available to drive the wind (neglecting the frequency dependence of the line transitions), and the effective temperature ($T_{\mathrm{eff}}$), which can be considered as a rough parametrization of the ionization state at the base of the wind and therefore the opacity. We strongly recommend to consult Tab. 9 for the scaling of $\dot{M}$ with physical stellar quantities for each of these algorithms for anything beyond simple order of magnitude estimates.

2.2 Combination of mass loss algorithms

The wind mass loss of massive stars is usually divided into three separate phases, whose definition is somewhat arbitrary. When using the algorithm K for the hot phase, we adopt the following thresholds based on the effective temperature $T_{\mathrm{eff}}$ and surface (i.e., outermost computational cell) hydrogen mass fraction $X_{s}$:

• $\bullet$

Hot phase: $T_{\mathrm{eff}}\geq 15\,000\ \mathrm{K}$;

• $\bullet$

Cool phase: $T_{\mathrm{eff}}<15\,000\ \mathrm{K}$;

• $\bullet$

WR phase: $X_{s}<0.4\ \mathrm{regardless\ of}\ T_{\mathrm{eff}}$.

To follow the suggestions of Glebbeek et al. (2009), and to have a smoother transition between the hot and cool phase wind algorithm, we use a slightly different definition of the cool and hot evolutionary stages when using the V mass loss algorithm for the hot phase:

• $\bullet$

Hot phase: $T_{\mathrm{eff}}\geq 11\,000\ \mathrm{K}$;

• $\bullet$

Cool phase: $T_{\mathrm{eff}}\leq 10\,000\ \mathrm{K}$;

and we use a linear interpolation between the hot phase wind and cool phase wind in between. We choose this different threshold when using V to match how the cool and hot phases are defined for the “Dutch” wind scheme in the MESA code. Note that the interval from $T_{\mathrm{eff}}\simeq 10,\mathrm{kK}$ to $\sim 15\,\mathrm{kK}$ is covered during a fraction of the Hertzsprung gap, in a very short time.

The threshold dividing the cool and hot phases is qualitatively justified as follows: the radiation pressure is determined by the product of opacity and flux, which peaks between $10000$–$15000\ \mathrm{K}$ because of iron recombination. Therefore, an effective temperature (i.e., in other words, a radius, for each given luminosity) in this range is a physically meaningful threshold to switch between wind mass loss algorithms.

The third phase is the WR phase and our criterion $X_{s}<0.4$ just requires a hydrogen-poor stellar surface. This has very little in common with the observational definition of what a WR star is, which is based on spectral features that are not tracked by stellar evolution codes. Specifically, WR stars are identified by their surface hydrogen depletion (Schmutz & Drissen, 1999) and the presence of broad emission lines (van der Hucht, 2001; Marchenko et al., 2010), indicating the presence of a wind with a steep density and velocity gradient. Moreover, typical WR stars have high $T_{\mathrm{eff}}$, and our definition might, in principle, produce unrealistically cold (and red) WR stars. However, in absence of strong mixing processes (e.g., due to rotation), the required surface hydrogen depletion can only be reached by removing mass from the surface and revealing deep and hot stellar layers. We do not find in our calculations cool but hydrogen depleted models. WR stars are further subdivided into classes (WNH, WN, WC, WO, etc.) based on the relative flux of specific lines. Here, we do not attempt to distinguish between different WR sub-classes, because our simulations do not produce the stellar spectra that would be necessary to distinguish these sub-classes (see however Meynet & Maeder, 2003; Groh et al., 2014, and references therein).

Our definitions of the evolutionary phases are commonly used in the literature (see, e.g., Limongi & Chieffi, 2006; Eggenberger et al., 2007; Woosley & Heger, 2007). We list in Tab. 2 the algorithm combinations explored here: each combination is labeled by combining the abbreviations for each wind algorithm introduced above.

To study the effects of the possible overestimate of the mass loss rate caused, for example, by wind inhomogeneities (i.e., “clumpiness”), we use three different values of the wind efficiency $\eta=1.0,\,0.33,\,0.1$. These values span the range of observational estimates of the volume filling factor of clumps (see Smith 2014 for a review). For simplicity, we use the same $\eta$ during the entire evolution of a given model.

2.3 The Grid of Stellar Models

We employ release version 7624 of the open-source stellar evolution code MESA (Paxton et al., 2011, 2013, 2015) and compute a grid of nonrotating solar metallicity models. We choose $Z_{\odot}=0.019$ to match precisely the value adopted in Vink et al. (2001). We consider initial masses of $15$, $20$, $25$, $30$, and $35\,$ $M_{\odot}$. Higher initial mass models are more strongly affected by numerical (and possibly physical) instabilities (see Appendix B and references therein), which makes them less reliable for our purpose. Appendix B describes the details of our MESA simulations that are not directly related to mass loss. Here, we only mention that we use a 45-isotope nuclear reaction network (mesa_45.net) until oxygen depletion, and a mixing length parameter $\alpha_{\mathrm{mlt}}=2.0$ with exponential overshooting.

We run our grid of models in three steps and at each checkpoint we make a selection of models to run at the next step222The models saved at each checkpoint are available at https://zenodo.org/record/292924#.WK_eENWi60i . The input files and customized routines are available at https://stellarcollapse.org/renzo2017.. This selection is necessary to reduce the total computational cost of our grid of models. In the first step, we evolve our models from the zero age main sequence (ZAMS) to when the temperature in the central computational cell rises above $T_{c}\geq 10^{9}\,\mathrm{K}$ (“end of the mass loss phase”). At this point, the star has little time left to live (a few years) and lose mass through winds (for a typical mass loss rate the star loses less than $10^{-4}\,M_{\odot}$). Also, when $T_{c}\geq 10^{9}\,\mathrm{K}$ is reached, MESA starts to artificially damp the mass loss rate for stability reasons. Mass loss is completely shut off when $T_{c}\geq 2\times 10^{9}\,\mathrm{K}$.

Note, however, that the observation of SN impostors, type IIn SNe, and intra-night flash spectroscopy of normal type II SNe suggests that non-wind mass loss phenomena (neglected here) may occur in the very late phase of the stellar life (see, e.g., Quataert & Shiode, 2012; Smith & Arnett, 2014; Quataert et al., 2016; Khazov et al., 2016).

In the second step, we restart a subset of MESA models and run to oxygen depletion (see Sec. 3.5). This is defined as the first time when the mass fraction of ${}^{16}\mathrm{O}$ in the central cell drops below 0.04 and the mass fraction of ${}^{28}\mathrm{Si}$ is higher than 0.01 (indicating that some oxygen burning has already occurred, following Sukhbold & Woosley, 2014). Note that these thresholds are an artificial choice, since the evolution of the star is continuous throughout its lifetime.

Finally, we restart a select subset of stars at oxygen depletion and run to the onset of core collapse, defined by

[TABLE]

where $v$ is the radial infall velocity (see, e.g., Heger et al., 2005; Sukhbold & Woosley, 2014). For this last phase of the evolution, we switch from the 45-isotopes nuclear reaction network to a customized 203-isotopes nuclear reaction network (see also Appendix B), to capture the details of the weak interaction physics (electron captures and $\beta$-decays) occurring during silicon burning and before collapse. This physics determines the final number of electrons per nucleon $Y_{e}$, thus the effective Chandrasekhar mass of the stellar core, and ultimately the core structure at the onset of core collapse. Experiments with a single-zone model for silicon burning show that at least $\sim 200$ isotopes are required to obtain a converged final value of $Y_{e}$ in the core (i.e., one that is independent of the size of the nuclear reaction network; see also Farmer et al., 2016). Since MESA solves the fully coupled set of equations for the chemical composition and structure of the star, the increase in the number of isotopes forces us to reduce the number of spatial mesh points in each model because of memory limitations (see Appendix B.1.3).

3 Results

3.1 Overview at the end of the mass loss phase, $T_{c}\geq 10^{9}\,\mathrm{K}$

The lifetime remaining for the star after $T_{c}\geq 10^{9}\,\mathrm{K}$ is short ($\sim 15$ years for $15\,M_{\odot}$, $\sim 4$ years for $30\,M_{\odot}$), and the photosphere of the star remains frozen at the same effective temperature $T_{\mathrm{eff}}$ and luminosity $L$ as long as no impulsive mass loss events or other instabilities take place. The subsequent evolution will lead to changes of the internal structure only. Therefore, we discuss the final mass and pre-SN appearance of our models already at this stage.

Figure 1 shows the overall impact of the uncertainty in stellar wind mass loss rates on the final mass. We plot the final mass relative to the initial mass for each considered $M_{\mathrm{ZAMS}}$ and the large vertical spread illustrates the uncertainty. As an example, we can consider a star of $M_{\mathrm{ZAMS}}=20\,M_{\odot}$. It can, in principle, reach the onset of collapse with $M=19.38\,M_{\odot}$ if evolved with the V-NJ mass loss combination with reduced efficiency $\eta=0.1$. But it might also evolve to $M=8.81\,M_{\odot}$ with the K-vL combination and efficiency $\eta=1.0$. If these are truly limiting cases and anything in between is unconstrained, then the uncertainty in the final mass is greater than $50\%$. For more quantitative results, see Tab. 5 where we report the maximum spreads of the ZAMS to pre-SN mass mapping for our entire grid of models. Tab. 3 and Tab. 4 list the main physical quantities for all the stars in our grid at the end of the mass loss phase.

It is important to note that the vertical spread in the ZAMS to pre-SN mass mapping shown in Fig. 1 is dominated by the highly uncertain wind efficiency $\eta$. At fixed $\eta$, variations due to different wind algorithm combinations are minor. This makes $\eta$ the most important free parameter for wind mass loss in stellar evolution calculations.

In Fig. 1, the spread in $M_{\mathrm{ZAMS}}$ to pre-SN mass decreases for higher $M_{\mathrm{ZAMS}}$. This is, however, only because we show the relative final mass, i.e. at higher $M_{\mathrm{ZAMS}}$ we divide the final mass by a larger number. In absolute numbers, the uncertainty in the final to initial mass relation increases for more massive stars. We summarize this in Tab. 3. As an example, the maximum spread between the final masses of $35\,M_{\odot}$ models is $\max{\Delta M}\equiv 33.99\,M_{\odot}-19.48\,M_{\odot}=14.51\,M_{\odot}$. In the $15\,M_{\odot}$ case, it is $\max{\Delta M}\equiv 14.66\,M_{\odot}-5.25\,M_{\odot}=9.42\,M_{\odot}$. As expected, all wind mass loss algorithm combinations yield a higher mass loss rate for more massive (and thus more luminous) stars. For stars with $M_{\mathrm{ZAMS}}\lesssim 20\,M_{\odot}$, the models using the vL algorithm (dust driven mass loss in the cool evolutionary phase) produce much higher mass loss than all other algorithms (see also Sec. 3.3): for example, the $15\,M_{\odot}$ model using the V-vL combination with full efficiency $\eta=1.0$ results in a pre-SN mass of only $5.25\,M_{\odot}$, while using the combination V-dJ or V-NJ, we obtain a final mass of $\sim 12.7\,M_{\odot}$. We discuss this effect in more detail in Sec. 3.3.

From Fig. 1, we also note that for $M_{\mathrm{ZAMS}}\gtrsim 20\,M_{\odot}$ we obtain a range $1.0\lesssim M/M_{\mathrm{ZAMS}}\lesssim 0.5$. For any given $\eta$, the range is smaller. In any case, the size of this interval does not decrease going to higher masses, indicating that the different mass loss prescriptions and efficiencies do not appear to be converging with $M_{\mathrm{ZAMS}}$ in the mass range we consider here.

In Fig. 2, we show, as a representative example, the time evolution of the total mass for our $15\,M_{\odot}$ models. The overall qualitative behavior is the same for all the other considered $M_{\mathrm{ZAMS}}$. Fig. 2 shows that the amount of mass lost during the main sequence (i.e. prior to the vertical dot-dashed line) is relatively small, only a few percent of the total mass for $M_{\mathrm{ZAMS}}=15\,M_{\odot}$, even when using $\eta=1.0$. Most of the mass is lost after hydrogen core burning. Both the mass loss rate and the spread between the predictions of different algorithms increases dramatically after the end of the main sequence. Hence, uncertainties in the post-main-sequence cool phase mass loss rates have by far the greatest effect on the final mass of the star.

Using our nonrotating single-star models at the end of the mass loss phase and ignoring any potential subsequent non-wind mass loss events, we can attempt to classify their pre-SN color. We follow Georgy (2012) by assuming $\log_{10}(T_{\mathrm{eff}}/\mathrm{[K]})\lesssim 3.6$ for RSG, $3.6\lesssim\log_{10}(T_{\mathrm{eff}}/\mathrm{[K]})\lesssim 3.8$ for YSG, $\log_{10}(T_{\mathrm{eff}}/\mathrm{[K]})\gtrsim 3.8$ for BSG, and we only require the surface abundance of hydrogen to be $X_{s}<0.4$ without temperature thresholds for WR stars. We list the outcome of this classification in Tab. 4. Almost all of our $15\,M_{\odot}$ and $20\,M_{\odot}$ models end their evolution as RSGs. The exceptions are those using the vL cool mass loss algorithm and $\eta=1.0$, which end their lives as YSGs. All of our $25\,M_{\odot}$ and $30\,M_{\odot}$ models end as YSGs. The $35\,M_{\odot}$ models computed with $\eta=0.1$ also end as YSGs. With $\eta=0.33$ they instead become BSGs, because of the higher mass loss rate, except when using the vL cool mass loss rate (cf. Sec. 3.3), which produces YSGs. With $\eta=1.0$ we find WR pre-SN models, unless vL is used during the cool phase, in which case our $35\,M_{\odot}$ models would explode as BSGs. However, these results are highly dependent on the somewhat arbitrary temperature thresholds assumed to divide the categories. For example, assuming $\log_{10}(T_{\mathrm{eff}}/\mathrm{[K]})\leq 3.68$ as the threshold dividing RSG and YSG, all models with $M_{\mathrm{ZAMS}}\leq 25\,M_{\odot}$ would be RSGs.

3.2 The “hot phase” mass loss

“Hot phase” evolution, i.e. $T_{\mathrm{eff}}\gtrsim 15\,000\,\mathrm{K}$ if using the Kudritzki et al. (1989) (K) mass loss algorithm, or $T_{\mathrm{eff}}\gtrsim 11\,000\,\mathrm{K}$ if using Vink et al. (2000, 2001) (V), roughly covers the main sequence evolution, the subsequent overall contraction caused by hydrogen depletion (known as the main sequence “hook”), and the initial part of the Hertzsprung gap.

In Tab. 6, we summarize key quantitative results for our models at the end of the hot phase. While this is the longest phase of the evolution, covering $\gtrsim 90\%$ of the stellar lifetime, the amount of mass lost during this phase is relatively small regardless of the algorithm used. For example, initially $35\,M_{\odot}$ stars computed with $\eta=1.0$ lose $\sim 15\%$ of their mass, while initially $15\,M_{\odot}$ stars only lose a few percent

We can also infer from Tab. 6 that the amount of mass lost during the hot phase is always higher with V than with K, and the difference increases with increasing $\eta$ and $M_{\mathrm{ZAMS}}$: it is only $\sim 0.02\,M_{\odot}$ for $\eta=0.1$ and $M_{\mathrm{ZAMS}}=15\,M_{\odot}$, and grows to $\sim 0.2\,M_{\odot}$ for $\eta=1.0$ and the same initial mass. For the $35\,M_{\odot}$ models and $\eta=1.0$, the difference between the total mass shed using V or K reaches $\sim 0.5\,M_{\odot}$ with $\eta=1.0$.

In Fig. 3, we plot the mass loss rates $\dot{M}$ given by V and K as a function of time in the hot phase for models with $\eta=1$. The reason for the higher total mass lost with V is that this algorithm includes a detailed treatment of the bistability jump, which is an increase in the cross section for photon interactions caused by the recombination of ions driving the mass loss. This enhancement of the cross section happens when the effective temperature drops below $T_{\mathrm{eff}}^{\mathrm{jump}}\simeq 25\,000\,\mathrm{K}$, Vink et al. (2000). This is what causes the sudden tremendous increase in the mass loss rate in V models seen in Fig. 3. The subsequent drop in models with masses higher than $20\,M_{\odot}$ happens because these cross the bistability jump region twice during the contraction following core hydrogen depletion. Overall, the average mass loss rate $\langle|\dot{M}_{\mathrm{hot}}|\rangle$ of V models is driven up and surpasses that of K models on the main sequence, resulting in higher mass loss in V models. Note that the use of two different thresholds to separate the hot and cool phase of evolution for V and K does not significantly influence the duration of the hot phase: the gap between the two thresholds is covered in a fraction of the Herzsprung gap duration.

Different values of $\eta$ produce small ($\lesssim 2\%$) age differences: models remaining more massive (i.e. computed with lower $\eta$) evolve slightly faster. However, these differences are too small to potentially be used as observational tests for the wind efficiency.

Figure 3 also shows that at any point in time, the mass loss rate of more massive stars is higher, because they produce a higher photon flux to drive the wind. A factor of $\sim 2.3\left(\simeq 35/15\right)$ in initial mass translates to a difference of almost two orders of magnitude in the mass loss rate. The difference between the V and K rates is a non-monotonic function of the mass, because of the different functional dependencies of the two algorithms: for $20\,M_{\odot}$ models they are roughly equal before the surface cools enough for the bistability jump to occur; for $15\,M_{\odot}$ models, the K algorithm gives an initially higher mass loss rate, and, for $30\,M_{\odot}$ models, K mass loss is instead lower.

Although only a small amount of mass is lost, the hot phase mass loss can significantly influence the core evolution. This is because no shell sources decouple the surface from the convective core during most of this phase. The use of different algorithms during the hot phase can thus create small (seed) differences in the core structure, which may then be amplified by the subsequent evolution of the star and contraction of the core. These differences are small at the end of the hot phase, and we will discuss them at later stages in the evolution in Sec. 3.6 and 3.7.

In Tab. 6, we also list the helium core masses $M_{\mathrm{He}}$ of our models at the end of the hot phase. The effects of hot phase wind mass loss on $M_{\mathrm{He}}$ are less straightforward to interpret than its influence on the total mass. First, the value of $M_{\mathrm{He}}$ depends on the definition of the helium core. We define the outer edge of the helium core as the first location going inward where $X(^{1}\mathrm{H})<0.01$. Second, where this interface is at the end of the hot phase is very sensitive to mixing: depending on the mass and metallicity of the star (and the convective stability criterion adopted), deep convective shells can develop at the beginning of the Hertzsprung gap, and they shape the chemical composition profile and determine $M_{\mathrm{He}}$. Tab. 6 shows that the maximum spread in $M_{\mathrm{He}}$ increases with $M_{\mathrm{ZAMS}}$, starting from $\max(\Delta M_{\mathrm{He}})\simeq 0.1\,M_{\odot}$ for $15\,M_{\odot}$ models, up to $\max(\Delta M_{\mathrm{He}})\simeq 0.9\,M_{\odot}$ for $35\,M_{\odot}$ models. The spread in these values is almost entirely due to variations in $\eta$. The difference in $M_{\mathrm{He}}$ between models of same mass and $\eta$ (thus differing only in the use of V or K) is also increasing with increasing $M_{\mathrm{ZAMS}}$ but remains below $\sim 0.2\,M_{\odot}$. This trend directly reflects the larger uncertainties in the modeling of winds from more massive stars.

3.3 The “cool phase” mass loss

Regardless of efficiency $\eta$ and mass loss algorithm, most of the mass loss through stellar winds happens during the cool phase of the evolution. Fig. 2 demonstrates this clearly for the $15\,M_{\odot}$ models, which lose $2\%-60\%$ of their total initial mass during this phase, depending on the algorithm combination and wind efficiency. The increase in the mass loss rate from the hot phase can be understood in terms of the effective gravity of the star (although we stress that the algorithms compared here do not depend explicitly on it): for any given luminosity of a massive star, if the stellar surface is cool, necessarily its radius must be large, and thus it will be easier for matter to leave the gravitational potential well of the star. Moreover, at lower temperature the opacity tends to be higher because of recombination of ions and possibly dust formation, thus enhancing the wind driving.

Table 3 summarizes key quantitative results of our models at the end of the mass loss phase, including the total mass and the core masses. One striking result is that the dust-driven van Loon et al. (2005) (vL) mass loss algorithm results in significantly different total masses and core masses than the de Jager et al. (1988) (dJ) and the Nieuwenhuijzen & de Jager (1990) (NJ) algorithms. This can also be seen in Fig. 1, where wind combinations using vL for the cool phase produce different vertical spreads. These differences are strongest for $\eta=1.0$. The most extreme example are $15\,M_{\odot}$ models computed with $\eta=1.0$: regardless of the hot phase mass loss algorithm, they end their evolution with masses of $\sim 5-6\,M_{\odot}$ with vL in the cool phase, while they remain as massive as $\sim 11-12\,M_{\odot}$ with dJ or NJ. We find the opposite trend at the upper end of our mass range: a $30\,M_{\odot}$ star computed with $\eta=1.0$ and V during the hot phase reaches the end of the mass loss phase with a total mass of $\sim 15M_{\odot}$ if using either dJ or NJ in the cool phase, while it ends its life with $\sim 18\,M_{\odot}$ with vL.

On the one hand, the similarities between the dJ and the NJ rates are expected (see also Mauron & Josselin, 2011; Eldridge & Tout, 2004): both are semi-empirical rates derived from the same sample of observed stars. They differ only in the choice of the stellar variables used to parametrize $\dot{M}$. On the other hand, the vL algorithm is also semi-empirical, but based on the analysis of a different sample of stars assuming a dust-driven model of the wind, i.e. wind mass loss is not driven by photons impinging on metallic ions, but rather on dust particles. If a dust-driven (instead of line-driven) model of the wind is assumed, the resulting mass loss rate is generally higher, and much more $T_{\mathrm{eff}}$-dependent (see Tabs. 1 and 9).

The very steep dependence of the vL rate on $T_{\mathrm{eff}}$ causes the different evolution of models above and below $\sim 25\,M_{\odot}$. During the early stage of the cool phase, the vL rate is always much higher than the others, and the stellar wind described by this algorithm reveals the deeper and hotter layers of the star (see also Tab. 4). As $T_{\mathrm{eff}}$ increases, the vL mass loss rate decreases rapidly ($\propto T^{-6.3}$, see Tabs. 1 and 9), which is attributed to the temperature sensitivity of the microscopic dust formation processes (Wachter et al., 2002) that we of course do not track explicitly in our calculations. Also, for any mass loss process at a given luminosity, higher $T_{\mathrm{eff}}$ correspond to smaller radii, and thus higher effective gravity at the stellar surface. For $M_{\mathrm{ZAMS}}\gtrsim 25\,M_{\odot}$, the steep $T_{\mathrm{eff}}$-dependence leads to a self regulation of the vL rate. For lower initial masses, the vL rate is also initially higher than the dJ or NJ rate, but not high enough to reach the self-regulating regime: the vL rate remains higher than the dJ and NJ rates for the whole evolution and produces pre-SN structures of a much lower final mass than when dJ or NJ are used. This is summarized in Tab. 3.

The comparison of $M_{\mathrm{He}}$ listed in Tab. 3 for models with the same $\eta$ and hot phase mass loss (either V or K) reveals that the effect of the cool phase mass loss on the helium core mass is very small and almost negligible. We find the only appreciable differences when using $\eta=1.0$ and the vL cool mass loss rate, and they are only of order $\sim 0.01\,M_{\odot}$.

3.4 Models reaching the WR stage

Out of the 94 models computed to $T_{c}\geq 10^{9}\,\mathrm{K}$, only 8 reach the conditions to switch to the WR wind scheme. These are all $35\,M_{\odot}$ models computed with $\eta=1.0$, and none use the vL algorithm in the cool phase that precedes the WR phase (cf. Tab. 4). The lack of WR models using vL is explained by the self-damping of this mass loss scheme for more massive and thus more luminous stars (see Sec. 3.3).

The typical duration of the WR phase is $\sim 0.02-0.05\,\mathrm{Myr}$. The differences in duration of the WR phase with the Nugis & Lamers (2000) (NL) and Hamann et al. (1982, 1995); Hamann & Koesterke (1998) (H) algorithms are negligible. However, models computed with the dJ algorithm in the cool phase have WR phases that are systematically longer by a few ten thousand years than the corresponding models computed with the NJ algorithm. Moreover, the duration of the WR phase of models computed with the K algorithm in the hot phase is about a factor of $\sim 2.3$ longer than in models using the V algorithm. This is because $M_{\mathrm{ZAMS}}=35\,M_{\odot}$ models computed with $\eta=1.0$ and the V algorithm reach the end of the hot phase with a helium core that is $0.22\,M_{\odot}$ less massive than models computed with the K algorithm (cf. Tab. 6). Therefore, the subsequent evolution is slowed down, and the WR phase is reached slightly later.

The NL WR phase algorithm produces higher final masses than the H algorithm: the difference is $\sim 0.3$ $(\sim 0.8)\,M_{\odot}$ for models using the V (K) hot phase algorithm and does not depend strongly on cool phase mass loss (cf. Tab. 3). Note that these differences are very small fractions of the initial mass $M_{\mathrm{ZAMS}}=35\,M_{\odot}$ of these models (cf. Fig. 1).

The WR wind does not have a strong influence on $M_{\mathrm{He}}$ at the end of the mass loss phase. For stars with $M_{\mathrm{ZAMS}}\lesssim 40M_{\odot}$, the He core mass is determined well before the beginning of the WR phase, and the wind mass loss is not strong enough to dig into the He core directly. More massive stars with stronger winds, may become hydrogen depleted already during the core hydrogen burning phase. In that case, the WR mass loss rate might have an impact on $M_{\mathrm{He}}$ through the quasi-static response of the convective core to mass loss (Meynet et al., 1994; de Koter et al., 1997; Crowther et al., 2010; Bestenlehner et al., 2011). Also, in even more massive stars, the entire hydrogen-rich envelope can be lost to winds, making $M_{\mathrm{He}}$ the total mass of the star, and winds can then further reduce it (e.g. Woosley, 2016).

While we find no systematic effect of WR mass loss on $M_{\mathrm{He}}$, the situation is different for $M_{\mathrm{CO}}$ (cf. Tab. 3). NL models yield $M_{\mathrm{CO}}$ that are systematically lower by $\sim 0.05-0.1\,M_{\odot}$ than H models and the largest differences are between models that also use different hot phase mass loss algorithms. We find that models differing only in the WR algorithm have lower $M_{\mathrm{CO}}$ for higher final masses. For example, the V-dJ-NL model has $M_{\mathrm{CO}}=9.81\,M_{\odot}$ and final mass $M=20.03\,M_{\odot}$. In contrast, the V-dJ-H model has a higher $M_{\mathrm{CO}}=9.93\,M_{\odot}$, but a lower final mass of $M=19.73\,M_{\odot}$. However, we note that the trend that lower final masses correspond to higher $M_{\mathrm{CO}}$ holds for most of our $35\,M_{\odot}$ models, independent of if they become WR stars or not.

The differences in both total mass and $M_{\mathrm{CO}}$ found varying the WR mass loss algorithm are more sensitive to the previously employed hot phase mass loss algorithm than to the cool phase algorithm. Moreover, $M_{\mathrm{He}}$ is also almost insensitive to the cool phase (cf. Sec. 3.3) and WR phase mass loss. Therefore, the differences found here are most likely related to the differences in $M_{\mathrm{He}}$ at the end of the hot phase (and in the total mass for a given $M_{\mathrm{He}}$), and consequently the position of the hydrogen burning shell, which indirectly influences the helium burning, the mixing processes shaping the composition profile, and ultimately the resulting $M_{\mathrm{CO}}$ and the amount of mass lost during the WR phase.

3.5 Models at oxygen depletion

To reduce the computational cost of our model grid, we select a subset of 44 stars to continue until oxygen depletion, which we define as the time when $X_{c}(^{16}\mathrm{O})\leq 0.04$ and $X_{c}(^{28}\mathrm{Si})\geq 0.01$. These models span the range of properties found at the end of the mass loss phase, and are listed in Tab. 7 together with their properties at oxygen depletion. This selection allows us to avoid running multiple models that have very similar evolutionary paths from the end of the mass loss phase onward. The duration of the evolution between the end of the mass loss phase ($T_{c}\geq 10^{9}\,\mathrm{K}$) and oxygen depletion is of order years to decades, depending on the total mass and core masses at the end of the mass loss phase.

In the very short time to oxygen depletion, neither total mass nor helium core mass change appreciably. This can be inferred by comparing the entries of Tab. 3 and Tab. 7, which also reveals that the stellar radii vary only within $\pm 3\,R_{\odot}$ in most models.

The CO core masses at oxygen depletion summarized in Tab. 7 are systematically a few percent lower than those listed in Tab. 3 at the end of the mass loss phase. This seems counter-intuitive, since one would expect the CO core to grow in mass because of the ashes of helium shell burning. However, the boundary location for the CO core is determined by convective mixing within and above the He burning shell, which brings helium-rich material inward and moves the CO core boundary to a smaller mass coordinate. This implies that the core mass is very sensitive to the mixing parameters. The maximum spread in $M_{\mathrm{CO}}$ obtained varying the wind algorithm is of order $\sim 0.1\,M_{\odot}$ and it increases with $M_{\mathrm{ZAMS}}$ up to about $\max(\Delta M_{\mathrm{CO}})\simeq 0.5\,M_{\odot}$ for $30\,M_{\odot}$ models (see also Tab. 5). Note that the $15\,M_{\odot}$ models are outliers in that they show a larger spread of CO core masses of up to $0.28\,M_{\odot}$ between models using the V and K hot mass loss algorithm. This is a consequence of a combination of (i) the V algorithm leading to a lower total mass and a lower He core mass at the end of the hot phase and (ii) the relatively low mass loss during the cool phase (compared to more massive stars), which results in deeper convective episodes. Together, these lead to small CO cores: the two $15M_{\odot}$ models using the V algorithm with $\eta=1.0$ have $M_{\mathrm{CO}}\simeq 2.95M_{\odot}$ at the end of the mass loss phase, which is about $0.2\,M_{\odot}$ smaller than the average for $15\,M_{\odot}$ models.

3.6 The Compactness Parameter

$\xi_{\mathcal{M}}$

Although the internal structure of our models is not yet final at oxygen depletion, some quantities (e.g., core masses) are already close to their final values, and it becomes possible to draw first connections between internal structure and the potential final outcome of core collapse (Sukhbold & Woosley, 2014). For this, we include the compactness parameter $\xi_{2.5}^{\mathrm{O\,depl}}$ in Tab. 7. O’Connor & Ott (2011) define the compactness parameter as

[TABLE]

This parameter provides a single measure for the complex inner core structure (mass coordinate smaller than $\mathcal{M}$) of a star, allowing for a simplified discussion of the differences in the internal structure produced by the various mass loss algorithms we compare.

We set $\mathcal{M}=2.5\,M_{\odot}$ because this is the typical mass above which the proto-NS that will form during core collapse will become a BH (O’Connor & Ott, 2011). This mass cut remains well outside of the typical iron core mass and includes the layers of the star that the shock will encounter after core bounce. These layers determine the accretion ram pressure that the shock has to overcome for a successful explosion. The qualitative behavior of the supernova dynamics and outcome with $\xi_{\mathcal{M}}$ is known to be robust against different choices of $\mathcal{M}$ (O’Connor & Ott, 2011, 2013; Ugliano et al., 2012; Sukhbold & Woosley, 2014).

One-dimensional parametric core-collapse SN explosion simulations (O’Connor & Ott, 2011, 2013; Ugliano et al., 2012; Ertl et al., 2016; Sukhbold et al., 2016) show that the value of $\xi_{2.5}$ at the onset of core collapse indicates the most probable remnant. High values of $\xi_{2.5}$ indicate a more compact pre-SN structure that is harder to explode and that will more likely result in a BH remnant. Conversely, low values of $\xi_{2.5}$ indicate a steeper density gradient and an easier to explode structure, suggesting that the remnant will more likely be a NS (O’Connor & Ott, 2011; Ugliano et al., 2012; Clausen et al., 2015). Sukhbold & Woosley (2014) suggest that the value of the compactness parameter evaluated at oxygen depletion, $\xi_{2.5}^{\mathrm{O\ depl}}$, can already be used to infer the most likely outcome of the core collapse event. The evolution from oxygen depletion to core collapse tends to increase the compactness and amplify the differences between different stellar models, but the key features that determine the interpretation of $\xi_{2.5}$ appear to be set already at oxygen depletion. Other parameters to relate the pre-SN structure to the most-likely remnant can be defined in the context of neutrino-driven explosions (see Ertl et al., 2016), but they rely on physical quantities that are not at all set at earlier stages of the evolution (e.g., the entropy profile throughout the silicon layer and iron core), and therefore they are not useful diagnostics before the onset of core collapse.

3.6.1 Evolution of $\xi_{2.5}$ until oxygen depletion

The compactness parameter is a function of time, $\xi_{\mathcal{M}}\equiv\xi_{\mathcal{M}}(t)$, because of the changes in the radius of a given mass coordinate $\mathcal{M}$. These can be caused by contraction of the core, onset of partial electron degeneracy within the mass coordinate $\mathcal{M}$, and by episodes of convective mixing and shell burning (Sukhbold & Woosley, 2014). The top panel of Fig. 4 shows examples of the evolution of $\xi_{2.5}$ until oxygen depletion in our $25\,M_{\odot}$ models. Note that we use a reversed logarithmic scale on the x axis to emphasize the late evolutionary stages.

Figure 4 shows that the compactness parameter is constant during the main sequence evolution. During this phase, it is also almost independent of $M_{\mathrm{ZAMS}}$ and mass loss algorithm because all stars considered here have convective main-sequence cores that are always much larger than the mass coordinate at which we evaluate the compactness. After core hydrogen exhaustion, $\xi_{2.5}$ increases because of the overall contraction, reaching $\xi_{2.5}\simeq 0.02$ in our $25\,M_{\odot}$ models. Then it slowly continues to increase during the hydrogen-shell burning and helium core burning phases. The increase speeds up significantly during core carbon burning, reaching values of $\xi_{2.5}\simeq 0.1$ in our $25\,M_{\odot}$ models. Neon core burning ignition and the onset of carbon shell burning mark a critical point in the evolution of $\xi_{2.5}$ at which the various curves in Fig. 4 begin to diverge. Sukhbold & Woosley (2014) find the same and point out that the subsequent evolution of $\xi_{2.5}$ is highly sensitive to the details of carbon shell burning (i.e. the number, locations, and durations of shell burning episodes). It is important to note from Fig. 4 that the effects of mass loss on core structure (represented by $\xi_{2.5}$) are delayed: At the time mass loss ends in our models (at $T_{c}>10^{9}\,\mathrm{K}$), differences in $\xi_{2.5}$ are minute. These seed differences grow and become substantial only in the last decade before core collapse.

The bottom panel of Fig. 4 depicts central density–temperature tracks for our $25\,M_{\odot}$ models that are evolved to oxygen depletion. The tracks start roughly at neon core ignition and show that the mass-loss history (i.e. the choice of wind mass loss algorithm combination) also influences the innermost core thermodynamics and structure. This is because the nuclear burning processes in the core are regulated by the amount of mass that needs to be sustained by the core itself, i.e. the mass below the innermost burning shell above the core. This in turn depends on the location and luminosity of the shell burning regions and therefore on the total mass of the star.

In Fig. 5, we show the values of $\xi_{2.5}^{\mathrm{O\ depl}}$ for all models that we run to oxygen depletion. The spread in each panel is due to the different algorithmic treatments of wind mass loss: for example, $25\,M_{\odot}$ models show values ranging between 0.210 and 0.157. Generally speaking, the spread in $\xi_{2.5}^{\mathrm{O\ depl}}$ increases with increasing $\eta$ and $M_{\mathrm{ZAMS}}$, i.e. the stronger the stellar wind, the more it influences the core structure. We emphasize that a few percent variation of $\xi_{2.5}^{\mathrm{O\ depl}}$ can result in important differences in the core structure at the onset of core collapse: the subsequent contraction of the core and the details of carbon, oxygen, and silicon shell burning, amplify the differences between models that are still relatively similar at oxygen depletion (see Sec. 3.7 and Sukhbold & Woosley, 2014).

3.6.2 Effects of the wind efficiency on $\xi_{2.5}^{\mathrm{O\ depl}}$

The effects of varying $\eta$ on $\xi_{2.5}^{\mathrm{O\ depl}}$ can be inferred from the comparison of models in Tab. 7 of the same $M_{\mathrm{ZAMS}}$ using the same wind algorithm combinations but different efficiencies. The variations of $\xi_{2.5}^{\mathrm{O\ depl}}$ with $\eta$ are typically non-monotonic: for example, the $25\,M_{\odot}$ model computed using the V-vL combination reaches oxygen depletion with $\xi_{2.5}^{\mathrm{O\ depl}}=0.179$, 0.164, and 0.211 for $\eta=0.1$, 0.33, and 1.0, respectively. Within the framework of each wind algorithm, higher values of $\eta$ correspond to higher mass loss rates and thus to a progressive shift of the evolution toward that of lower initial mass. However, the compactness parameter is known to be a highly non-monotonic function of $M_{\mathrm{ZAMS}}$ (Sukhbold & Woosley, 2014). Therefore, a higher mass loss rate can sometimes result in a decrease, an increase, or even almost no variation of $\xi_{2.5}^{\mathrm{O\ depl}}$. For example, in $15\,M_{\odot}$ models computed with V-vL, $\xi_{2.5}^{\mathrm{O\ depl}}$ decreases when going from $\eta=0.1$ to $\eta=1.0$. In $25\,M_{\odot}$ models computed with V-vL, we instead find that $\xi_{2.5}^{\mathrm{O\ depl}}$ increases with the same change in $\eta$. And in the case of $20\,M_{\odot}$ models computed with V-dJ, we find only tiny variations of $\xi_{2.5}^{\mathrm{O\ depl}}$ when $\eta$ is varied from $0.1$ to $1.0$. See Tab. 7 for more examples and details.

3.6.3 Effects of varying the mass loss algorithm on

$\xi_{2.5}^{\mathrm{O\ depl}}$

The last column of Tab. 7 shows that the effects of different wind mass loss algorithms on $\xi_{2.5}$ are in most cases small until oxygen depletion for $\eta=0.1$ and $\eta=0.33$. For $\eta=0.1$, the V algorithm generally results in higher values of $\xi_{2.5}^{\mathrm{O\ depl}}$ than the K algorithm. This holds for all studied ZAMS masses with the exception of the $30\,M_{\odot}$ models. These do not exhibit this trend because all their burning shells are outside the mass coordinate $\mathcal{M}=2.5\,M_{\odot}$, as can be seen from the values of $M_{\mathrm{CO}}$ listed in Tab. 7.

The spread in $\xi_{2.5}^{\mathrm{O\ depl}}$ between different mass loss algorithms increases for models with $\eta=0.33$. For example, Fig. 5 shows that the $25\,M_{\odot}$ model with the K-dJ combination reaches $\xi_{2.5}^{\mathrm{O\ depl}}\simeq 0.21$ (similar to its $20\,M_{\odot}$ counterpart, cf. Tab. 7), which is $\sim 30\%$ higher than with other algorithm combinations.

Models computed with $\eta=1.0$ are most suitable to discuss the effect of different wind mass loss algorithm combinations. Both the hot phase and the cool phase mass loss algorithms influence $\xi_{2.5}^{\mathrm{O\ depl}}$, but their detailed effect varies with $M_{\mathrm{ZAMS}}$ and is strongest in the $25\,M_{\odot}$ and $30\,M_{\odot}$ models. For example, from Tab. 7, we find $\xi_{2.5}^{\mathrm{O\ depl}}\sim 0.16$, for the $\eta=1.0$, $25\,M_{\odot}$ models with the V-NJ and K-dJ combinations, while models with V-vL and K-vL result in $\xi_{2.5}^{\mathrm{O\ depl}}\sim 0.2$. At lower $M_{\mathrm{ZAMS}}$, even for $\eta=1.0$, differences in $\xi_{2.5}^{\mathrm{O\ depl}}$ due to the choice of mass loss algorithm combination are overall (with few exceptions) rather small and typically at the few percent level. While these differences will be amplified by the subsequent evolution toward core collapse, they are small compared to the tremendous differences in total mass resulting from the different algorithm combinations (cf. Tab. 7).

An interesting question to address is the relative importance of hot phase (i.e. main sequence) and cool phase (i.e. post main sequence) mass loss for $\xi_{2.5}^{\mathrm{O\ depl}}$. Naively, one would think that by the time the core and envelope are essentially decoupled, loss of envelope mass should have limited impact on the subsequent evolution of the core. Our results suggest that this is not generally the case.

The limited set of models run to oxygen depletion and listed in Tab. 7 give a complex, but necessarily incomplete picture of the relative importance of each mass loss phase for the core structure. Which phase is most relevant depends on $M_{\mathrm{ZAMS}}$ and $\eta$. For brevity, we focus here on the $\eta=1$ case and compare models evolved with the same hot phase mass loss algorithm (V or K) and different cool phase algorithms.

For $15\,M_{\odot}$ models with $\eta=1$, the tremendous mass loss with the vL algorithm in the cool phase has little effect on He and CO core masses and on $\xi_{2.5}^{\mathrm{O\ depl}}$. For example, the final masses of K-vL and K-NJ are $5.70\,M_{\odot}$ and $11.87\,M_{\odot}$, respectively. Yet their $\xi_{2.5}^{\mathrm{O\ depl}}$ are very close to each other, $0.152$ and $0.153$, respectively. Qualitatively, the same is true for the V-NJ and V-vL combinations. On the other hand there is a larger spread between combinations using V and K, suggesting that the small differences seeded by hot phase mass loss dominate in the $15\,M_{\odot}$ $\eta=1$ case. Similarly, we find for $30\,M_{\odot}$ models with $\eta=1$ that K-dJ and K-vL lead to final masses of $18.51\,M_{\odot}$ and $22.53\,M_{\odot}$, respectively, but their $\xi_{2.5}^{\mathrm{O\ depl}}$ are $0.244$ and $0.243$. The situation is more complicated for $25\,M_{\odot}$ and $20\,M_{\odot}$ models that straddle the $M_{\mathrm{ZAMS}}$ range where $\xi_{2.5}$ varies chaotically (Sukhbold & Woosley, 2014). From Tab. 7 we see that for the $20\,M_{\odot}$, $\eta=1$ models all considered K- and V- combinations yield roughly the same $\xi_{2.5}^{\mathrm{O\ depl}}$. In $25\,M_{\odot}$, $\eta=1$ models, on the other hand, cool phase mass loss has the dominant impact on $\xi_{2.5}^{\mathrm{O\ depl}}$. For example, K-dJ and K-vL have compactness of $0.161$ and $0.200$, respectively, although their final masses differ by only $\sim 1\,M_{\odot}$ due to the self-regulation of the vL algorithm in this mass range.

3.6.4 Comparison with Sukhbold & Woosley (2014)

Sukhbold & Woosley (2014) employed the NJ mass loss algorithm without any efficiency scaling factor (i.e. at efficiency $\eta=1.0$) throughout the entire evolution of their models. Comparing our Tab. 7 with their Fig. 23, we find that the compactness parameter values at oxygen depletion of our models lie in the same range as theirs, with a tendency toward slightly higher values.

Our $15\,M_{\odot}$ models produce values of $\xi_{2.5}^{\mathrm{O\ depl}}\simeq 0.15$, which is slightly higher than their value of $\sim 0.11-0.13$, especially for reduced wind mass loss rates (i.e. $\eta<1.0$). For this initial mass, increasing $\eta$ decreases the compactness of the core and reduces the difference between our models and those of Sukhbold & Woosley (2014).

Most of our $20\,M_{\odot}$ models have $\xi_{2.5}^{\mathrm{O\ depl}}\simeq 0.18$, close to the corresponding models of Sukhbold & Woosley (2014). For these models, the maximum difference varying the wind mass loss algorithm is $\Delta\xi_{2.5}^{\mathrm{O\ depl}}\lesssim 0.05$. This is not surprising, because large variations of $\xi_{2.5}$ are expected because of the transition from convective to neutrino-cooled and radiative carbon shell burning, which happens around $M_{\mathrm{ZAMS}}\simeq 20\,M_{\odot}$ (Sukhbold & Woosley, 2014). Therefore, in this mass range, changing the wind mass loss algorithm can substantially change $\xi_{2.5}^{\mathrm{O\ depl}}$ by shifting the evolutionary track of the star in the slightest way.

Our $25\,M_{\odot}$ models have values of $\xi_{2.5}^{\mathrm{O\ depl}}$ similar to those of the $20\,M_{\odot}$ models, in agreement with Sukhbold & Woosley (2014). However, once again, we obtain a large variation of $\xi_{2.5}^{\mathrm{O\ depl}}$ changing the treatment of mass loss ($\Delta\xi_{2.5}^{\mathrm{O\ depl}}\lesssim 0.05$).

Most of our $30\,M_{\odot}$ models have $\xi_{2.5}^{\mathrm{O\ depl}}\simeq 0.23$, which is significantly larger than the corresponding value of $\sim 0.16$ found by Sukhbold & Woosley (2014). The variations of the compactness parameter with mass loss algorithm combination and efficiency are smaller for the $30M_{\odot}$ models than for lower $M_{\mathrm{ZAMS}}$, except for the model computed with $\eta=1.0$ and the V-dJ algorithm. This model has $\xi_{2.5}^{\mathrm{O\ depl}}=0.179$, which is much closer to the values of Sukhbold & Woosley (2014). We expect that the similar V-NJ algorithm combination with $\eta=1.0$ would produce a structure close to this model at oxygen depletion, based on the similarities between the models at the end of the mass loss phase. The relatively low compactness of this model is therefore likely a consequence of the V hot phase mass loss algorithm. For $M_{\mathrm{ZAMS}}\gtrsim 30\,M_{\odot}$, V with $\eta=1.0$ produces substantially more mass loss than K (cf. Tab. 6), and therefore, the subsequent evolution is closer to the path of less massive stars – which are also expected to reach oxygen depletion with a lower compactness.

We speculate that the quantitative differences between our findings for the compactness parameter at oxygen depletion and those of Sukhbold & Woosley (2014) are in fact due primarily to their choice of mass loss algorithm: they employ the NJ algorithm in both the hot and the cold phase, which results in overall greater early mass loss than the K algorithm (and the V algorithm until the bistability jump). This notion is corroborated by our finding that our $30\,M_{\odot}$ model closest to theirs is the one that loses the most mass during the hot phase (using $\eta=1$ and the V-dJ combination; cf. Tab. 6).

3.7 Models at the onset of core collapse

We select a subset of six of our models at oxygen depletion for continuation to the onset of core collapse. Reducing the model set is necessary to limit the computational cost of this study. We choose two $15\,M_{\odot}$ models with efficiency $\eta=1.0$ and mass loss algorithm combination V-NJ and K-VL, two $20\,M_{\odot}$ models with $\eta=0.33$ computed using the V-vL and the K-dJ combination, and two $30\,M_{\odot}$ models with $\eta=0.33$ and the V-dJ or the K-NJ combination. When restarting our models from oxygen depletion, we switch from the 45-isotopes nuclear reaction network used so far to a larger customized network with 203-isotopes (see Appendix B). This is necessary to capture core deleptonization due to electron capture during and after silicon burning. Continuing the evolution with a larger nuclear reaction network also requires reducing the number of computational mesh points (from $\sim 10^{4}$ to $\sim 10^{3}$) to run the simulations within the memory constraints of MESA (see Appendix B for details). By the time oxygen depletion is reached, the effect of wind mass loss on the core structure is already pronounced, and our reduced-resolution models still have spatial resolution that is comparable to that of published models (e.g., Woosley et al., 2002; Woosley & Heger, 2007). Furthermore, resolution tests in Appendix B.1.3 give us confidence that the presently unavoidable reduction in resolution does not affect our overall results and conclusions.

Figure 6 shows the time evolution of $\xi_{2.5}$ (top panel) and the central temperature–central density evolutionary tracks (bottom panel) from oxygen depletion to the onset of core collapse for our pre-SN model set. The top panel shows that $\xi_{2.5}$ settles onto its final value before the criterion for the onset of core collapse is reached. The bottom panel clearly shows a hook at $\log_{10}(T_{c}/\mathrm{[K]})\simeq 9.55$, where $\rho_{c}$ decreases at roughly constant temperature, indicating the point of silicon core ignition. This is also the more “noisy” part of these tracks, indicating that this phase of nuclear burning with very high and nearly balancing reaction rates is the most challenging to simulate (Hix & Thielemann, 1996; Hix et al., 2007).

The evolution of the compactness parameter shows maxima and minima, which are related to oxygen shell ignition, around $\log_{10}(\{t_{\mathrm{pre-SN}}-t\}/\mathrm{[yr]})\simeq-2$, silicon core ignition around $\log_{10}(\{t_{\mathrm{pre-SN}}-t\}/\mathrm{[yr]})\simeq-3$, and silicon shell ignition at about $\log_{10}(\{t_{\mathrm{pre-SN}}-t\}/\mathrm{[yr]})\simeq-4.5$. However, note that the ignition times and the durations of these burning phases are mass-dependent. Figure 7 shows the corresponding increase in the neutrino emission from nuclear reactions for the $15\,M_{\odot}$ model computed with the V-NJ combination and $\eta=1.0$. Similar features are present in the $T_{c}(\rho_{c})$ evolutionary tracks. However, while the $T_{c}(\rho_{c})$ track only probes the innermost part of the stellar core, the $\xi_{2.5}$ evolution is determined by the interplay between silicon, oxygen, and carbon burning shells, core contraction, and onset of electron degeneracy.

Interestingly, the final compactness of the $15\,M_{\odot}$ models is lower at the onset of core collapse than at oxygen depletion: for example, the $15\,M_{\odot}$ model computed with K-VL and $\eta=1.0$ has $\xi_{2.5}^{\mathrm{pre-SN}}=0.132<\xi_{2.5}^{\mathrm{O\ depl}}=0.152$. This is because of the presence of nuclear burning shells within $\mathcal{M}=2.5\,M_{\odot}$ whose energy generation tends to expand the material in the layers above them (the same process happens during the Hertzsprung gap for hydrogen-shell-burning stars). The location of the shells can be estimated using the core masses listed in Tab. 8 (and Tab. 7): the $15\,M_{\odot}$ models have two shells of nuclear burning (Si and O, respectively) within the $2.5\,M_{\odot}$ mass coordinate, while only one shell exists in this region at oxygen depletion. Models with $M_{\mathrm{ZAMS}}>15\,M_{\odot}$ settle on a pre-SN compactness that is higher than the corresponding $\xi_{2.5}^{\mathrm{O\ depl}}$, because only the silicon burning shell is within $\mathcal{M}=2.5\,M_{\odot}$.

Table 8 lists the properties of the six models that we run to the onset of core collapse. $\xi_{2.5}^{\mathrm{pre-SN}}$, and ($M_{4}$, $\mu_{4}$), where $M_{4}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}M(s=4)$ is the mass location where the specific entropy is $s=4\,\mathrm{k_{b}\,baryon^{-1}}$ and $\mu_{4}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}dm/dr|_{s=4}$ is the mass gradient at that location, offer two different ways to estimate how hard it will be for the SN shock to unbind the stellar mantle and leave a NS remnant (see Ertl et al., 2016). The total mass at density higher than $10^{6}\,\mathrm{g\ cm^{-3}}$ ($M_{\rho 6}$), the carbon-oxygen core mass ($M_{\mathrm{CO}}$), and the iron core mass ($M_{\mathrm{Fe}}$), defined as the location where $X(^{28}\mathrm{Si})<0.01$) can be used to estimate the nickel yields of the possible SN explosion and the remnant mass (Fryer et al., 2012; Sukhbold et al., 2016).

The final $M_{\mathrm{CO}}$ and $M_{\mathrm{Fe}}$ depend, although weakly, on the mass loss algorithm adopted during the hot evolutionary phase (V or K). The V algorithm yields slightly smaller cores (and total masses at the end of the hot phase, cf. Sec. 3.2). The $15\,M_{\odot}$ models with $\eta=1.0$ reach core collapse with $M_{\mathrm{CO}}=2.91\,(3.07)\,M_{\odot}$, and $M_{\mathrm{Fe}}=1.39\,(1.50)\,M_{\odot}$ when using the combination V-NJ (K-vL). For $25\,M_{\odot}$ models with $\eta=0.33$, we find $M_{\mathrm{CO}}=6.38\,(6.40)\,M_{\odot}$, and $M_{\mathrm{Fe}}=1.51\,(1.63)\,M_{\odot}$ for V-vL (K-dJ). Finally, the $30\,M_{\odot}$ models with $\eta=0.33$ yield $M_{\mathrm{CO}}=7.98\,(7.90)\,M_{\odot}$, and $M_{\mathrm{Fe}}=1.56\,(1.58)\,M_{\odot}$ for the combination V-dJ (K-NJ). The differences in $M_{\mathrm{CO}}$ (and to a lesser extent $M_{\mathrm{Fe}}$) decrease with $M_{\mathrm{ZAMS}}$. However, note that the decreasing difference is most likely caused by the lower wind efficiency $\eta=0.33$ for the $25\,M_{\odot}$ and $30\,M_{\odot}$ models, while our $15M_{\odot}$ models use full efficiency, i.e. $\eta=1.0$.

As anticipated in Sec. 3.5, the spread in $\xi_{2.5}$ increases between oxygen depletion and the onset of core collapse: the final variations are about $\sim 30\%$ for models with the same initial mass (cf. Tab. 8). The two $15\,M_{\odot}$ models have $\Delta\xi_{2.5}^{\mathrm{O\ depl}}=0.011$, and $\Delta\xi_{2.5}^{\mathrm{pre-SN}}=0.029$. For the $25\,M_{\odot}$ models, the spread at oxygen depletion is $\Delta\xi_{2.5}^{\mathrm{O\ depl}}=0.046$, while it is $\Delta\xi_{2.5}^{\mathrm{pre-SN}}=0.081$ at the onset of core collapse. The two $30\,M_{\odot}$ models go from $\Delta\xi_{2.5}^{\mathrm{O\ depl}}=0.015$ to $\Delta\xi_{2.5}^{\mathrm{pre-SN}}=0.082$.

The pre-SN chemical abundances in the core are also affected by the choice of the mass loss algorithm combination. Fig. 8 shows a comparison of the distribution of the dominant isotopes for our pre-SN models. The composition of the oxygen rich layer is sensitive to the early (hot phase) mass loss, with the V scheme producing a lower ratio of $X(^{16}\mathrm{O})/X(^{20}\mathrm{Ne})$, owing to the higher early mass loss (cf. Sec. 3.2), and thus lower core temperature during the late phases (cf. bottom panel of Fig. 6). The distribution of the chemical elements in mass coordinate is also sensitive to the adopted mass loss algorithm combination. However, it is likely to depend more strongly on the treatment of mixing processes (which we do not vary here), mainly convection and overshooting, which are the only processes fast enough to have an effect inside the CO core of a star before the onset of collapse.

4 Discussion

Sensitivity of wind mass loss to evolving stellar properties.

For a given $M_{\mathrm{ZAMS}}$, the luminosity varies little between models that experience different mass loss rates and the effective temperature (and radius) varies much more. Once mass-loss induced differences in $T_{\mathrm{eff}}$ appear, they feed back on mass loss to amplify these differences. Hence, the $T_{\mathrm{eff}}$ evolution is most important for governing the mass loss of models of a given $M_{\mathrm{ZAMS}}$. The dependence of $T_{\mathrm{eff}}$ is particularly strong if the dust-driven vL algorithm is included.

Role of dust in stellar winds.

The role of dust as a driver of RSG winds is a subject of debate (e.g., van Loon et al., 2005; Ferrarotti & Gail, 2006; Bennett, 2010). Dust-driven mass loss might occur in only a part of RSG evolution, and possibly even under only rather specific conditions. Even if dust exists in the envelopes of cool, evolved, massive stars, uncertainties in the grain properties result in highly uncertain mass loss rates (van Loon et al., 2005). The extreme dust-driven vL mass loss seen in our models is just one example (see Sec. 3.3). In this context, it is important to mention that our models switch to the vL algorithm already at $T_{\mathrm{eff}}\lesssim 15\,000\,\mathrm{K}$ ($\lesssim 11\,000\,\mathrm{K}$) when K (V) is used in the hot phase. These temperatures are clearly far too hot for dust to form. However, the strong temperature dependence $(\sim T^{-6.3})$ in combination with the extremely short time spent in the Hertzsprung gap prevent this from having a substantial consequence for the stellar mass. Hence, the extreme mass loss we find for $M_{\mathrm{ZAMS}}\lesssim 20\,M_{\odot}$ with vL and $\eta=1$ is not an artifact of how we use the vL algorithm. Within its framework the predicted very low final masses for $15\,M_{\odot}$ and $20\,M_{\odot}$ stars, $5.25-5.70\,M_{\odot}$ and $8.8\,M_{\odot}$, respectively, are robust. It may thus be possible to rule out the extreme case of vL with $\eta=1$ using pre-explosion imaging and SN ejecta mass estimates.

We also note that in stars with $M_{\mathrm{ZAMS}}\gtrsim 20\,M_{\odot}$ the vL algorithm self-regulates since its extreme mass loss uncovers hotter layers of the stars. The steep temperature dependence of the vL algorithm then results in lower mass loss rates than dJ and NJ.

The role of the efficiency parameter $\eta$.

The wind efficiency $\eta$ is the parameter that has the greatest impact on the evolution of initially identical models (cf. Fig. 1 and Tab. 3). However, it lacks an interpretation from first principles. Here, we investigate only values $\eta\leq 1$, focusing on reduced wind mass loss motivated by possible inhomogeneities in the wind structure (also called “clumpiness”, Smith, 2014; Puls et al., 2008, and references therein). However, starting from first principles, it cannot be excluded that clumps might actually enhance the mass loss rate (Lucy & White, 1980). This might be the case if the overdense clumps are efficiently pushed outward by impinging photons.

Furthermore, we assume a constant efficiency factor throughout the evolution, but in principle the “clumpiness” of the wind might evolve (possibly even in a stochastic way) and may require different $\eta$ in different evolutionary phases. The mass loss routines in MESA333These are available at https://stellarcollapse.org/renzo2017 are already adapted for varying this parameter in different evolutionary phases. It is also possible that during each single evolutionary phase we define, the efficiency of mass loss may vary significantly, producing enhanced mass loss episodes separated by reduced mass loss phases.

Although changing $\eta$ induces substantial changes in the mass loss rate and final mass, the appearance (i.e., luminosity, effective temperature) is less affected, making it difficult to use observed stellar populations to constrain $\eta$ (Renzo, 2015).

Metallicity effects.

We do not investigate the effects of decreasing the metallicity in this study. However, an approximate picture can be drawn by considering models with reduced efficiency $\eta$ as proxy for low metallicity models. Most stellar evolution codes implement the metallicity dependence of the wind mass loss by just rescaling the mass loss rate at solar metallicity (cf. Eq. 1), which is exactly the purpose of $\eta$. While the metallicity at the surface of the star can change throughout the evolution, the main element driving a wind is iron (Vink et al., 2001; Tramper et al., 2016), and its abundance is unlikely to change because of upward mixing from the stellar interior. Therefore, it is not unrealistic to consider a constant metallicity-related reduction factor for the entire evolution. Nevertheless, the approach of considering reduced $\eta$ as a proxy for lower metallicity does not take into account metallicity effects on stellar radius and nuclear burning. These could indirectly affect wind mass loss.

WR stars.

We emphasize in Sec. 2.2 the shortcomings of the computational definition of WR stars adopted in stellar evolution codes ($X_{s}<0.4$). Although this definition is artificial, we stress that the mass loss algorithms used during this phase are derived from the observation of real WR stars. In the mass range considered here, few WR stars are expected, and indeed our results show that only $M_{\mathrm{ZAMS}}=35\,M_{\odot}$ stars with $\eta=1.0$ can become sufficiently hydrogen-depleted at their surface to switch to a WR mass loss algorithm. We compare only two WR mass loss algorithms (NL and H), and neglect algorithms obtained by fitting either very luminous (i.e. more massive) or hydrogen-free WR stars (e.g. Gräfener & Hamann, 2008; Tramper et al., 2016), since none of our WR models reach the corresponding regions of the parameter space. Therefore, in the framework of single nonrotating stars with wind efficiency $\eta=1.0$, the minimum ZAMS mass to obtain a WR model is somewhere between 30 and 35$\,M_{\odot}$. Although on the high end, this is in relatively good agreement with the results obtained with other stellar evolution codes (e.g. Woosley et al., 2002; Limongi & Chieffi, 2006; Eldridge & Vink, 2006; Georgy et al., 2015). Even with full efficiency of the wind before the WR phase, our models would likely underestimate the number of observable WR stars.

Lowering $\eta$ has the obvious consequence of decreasing the mass loss rate, and thus increasing the minimum ZAMS mass for single nonrotating WR stars. However, the standard picture of a single nonrotating star misses pieces of physics of great importance for the formation of WR stars (see, e.g., Maeder, 1996; Meynet & Maeder, 2003; Eldridge & Vink, 2006). These include, but are not necessarily limited to, rotationally-enhanced mass loss, rotational mixing processes (which can help depleting hydrogen from the surface, Meynet & Maeder, 2003; Maeder, 1996), and binary interactions. Binarity can lead to the formation of WR stars via envelope stripping in RLOF. Alternatively, accretion or merger with a companion could increase the mass (and luminosity) of the star sufficiently to enhance the wind mass loss and remove the hydrogen-rich envelope. Also note that the envelope hydrogen depletion needed to switch to a WR mass loss algorithm can be reached also because of upward mixing of thermonuclearly processed material, e.g., because of efficient rotational mixing, and not only because of mass loss. The choice for the algorithmic representation of mixing processes (see Appendix B) influences directly the surface mass fraction of hydrogen, but also the core mass, and consequently the luminosity. Indirectly, these effects can change the mass loss rate and consequently the fate of a model from/to WR.

Nucleosynthetic yields.

The nucleosynthetic yields of massive stars are mass loss (and angular momentum loss) dependent (see, e.g., Maeder, 1992; Frischknecht et al., 2016). Processes such as rotational mixing can bring thermonuclearly processed material upward that is then lost through winds. This is especially relevant for s-process elements which are synthesized during the hydrostatic lifetime of massive stars and the ratio of carbon to oxygen abundance in the stellar wind yields. On top of this, the success or failure of the SN explosion, and the details of the explosive nucleosynthesis, depend on the interior structure of the exploding star and thus on its mass loss history (e.g., Sukhbold et al., 2016).

Consequences for SN explosions and compact remnants.

We find that mass loss affects the core structure and the burning shells surrounding the core. Mass loss during the hot phase of the evolution (i.e. the main sequence, roughly speaking) is important for the core structure, because the core itself re-adjusts quasi-statically to the wind from the stellar atmosphere. During this phase, different algorithms produce small variations in the core, which are then amplified by the subsequent evolution (cf. Sec. 3.2 and Fig. 4). The general trend is that a higher mass loss rate during the hot phase produces structures with lower core compactness (cf. Sec. 3.6). The cool phase mass loss also impacts the core compactness, but more indirectly, through its effect on the burning shells. Most of the mass is lost during the cool phase of the evolution (cf. Sec. 3.3), and the vigor, extent and type of mixing in the burning shells all depend on the amount and the timing of mass loss. Since our present understanding of core-collapse SN explosions strongly depends on the details of the input stellar models (e.g., Janka et al., 2012; Couch & Ott, 2015; Chatzopoulos et al., 2016), overlooking the impact of wind mass loss on the core structure might bias detailed hydrodynamical simulations of stellar explosions. This, in turn, can have significant implications for the NS/BH ratio and mass distribution, and consequently also for the inferred gravitational wave sources. We provide444Data are available at https://zenodo.org/record/292924#.WK_eENWi60i and input parameter files at https://stellarcollapse.org/renzo2017 stellar models at the end of the main sequence, at the end of the hot phase of the evolution, at the end of the mass loss phase ($T_{c}\geq 10^{9}\,\mathrm{[K]}$), at oxygen depletion, and at the onset of core collapse. These models can be used as starting points for stellar experiments during late evolutionary phases (see e.g. Couch & Ott, 2015; Chatzopoulos et al., 2016).

Observational $M_{\mathrm{ZAMS}}$

estimates of SN progenitors.

The large vertical spread in final mass caused by differences in mass loss for a given $M_{\mathrm{ZAMS}}$ (Fig. 1) suggests at first sight that these large uncertainties may map to equally large uncertainties in $M_{\mathrm{ZAMS}}$ estimates from pre-explosion observations. However, the pre-SN mass is not a direct observable. Rather, $M_{\mathrm{ZAMS}}$ estimates typically rely on observational measurements of luminosity and effective temperature that are then compared with stellar models (e.g., Smartt 2009). In Fig. 9, we plot the final luminosity $L$ and effective temperature $T_{\mathrm{eff}}$ for our entire model set. Table 4 summarizes the numerical results. These results demonstrate that wind mass loss variations have very little effect on the final luminosity for stars with $M_{\mathrm{ZAMS}}\lesssim 30\,M_{\odot}$ and luminosity variations are smaller than typical observational uncertainties (see, e.g., Smartt 2009). The similarity in luminosity of models of a given $M_{\mathrm{ZAMS}}$ is a consequence of the rather small effect that mass loss has on the core mass (cf. Tab. 3). Interestingly, $T_{\mathrm{eff}}$ variations are also small for $M_{\mathrm{ZAMS}}\lesssim 30\,M_{\odot}$. The maximum variation from the average $T_{\mathrm{eff}}$ of a ZAMS mass is only $0.08\,\mathrm{dex}$ and comes from the $\eta=1$ K-vL and V-vL $20\,M_{\odot}$ models that are YSGs at the end of their lives. Wind mass-loss dependent variations in $L$ and $T_{\mathrm{eff}}$ are much larger for $35\,M_{\odot}$ models, some of which die as BSGs and some as WR stars.

Given the above results, it appears that for massive stars with $M_{\mathrm{ZAMS}}\lesssim 30\,M_{\odot}$ wind mass loss uncertainties do not increase the overall level of uncertainty with which the SN progenitor $M_{\mathrm{ZAMS}}$ can be estimated from pre-explosion observations. However, our results do not say anything about the other possibly existing mass loss channels (e.g., binary interactions and/or impulsive phenomena, see, for example, Smith, 2014; Smith & Arnett, 2014; Morozova et al., 2015; Margutti et al., 2017) that are usually neglected in stellar evolution calculations, or included in an very simplified way using enhanced winds (see, e.g., Meynet et al., 2015).

The RSG problem.

The observed lack of RSG SN progenitors with $M_{\mathrm{ZAMS}}\gtrsim 16\,M_{\odot}$ (Smartt et al., 2009; Smartt, 2009) might be explained with the effects of mass loss on the pre-SN stellar appearance (Georgy et al., 2016). Table 4 and Fig. 9 show that different wind algorithms can change the pre-explosion appearance of a massive star. However, this effect is relatively small in luminosity and $T_{\mathrm{eff}}$, and the relative number of YSG to RSG we find (cf. Sec. 3) largely depends on the adopted definition of YSG vs. RSG. Therefore, the systematic uncertainty in the treatment of wind mass loss for $16\lesssim M_{\mathrm{ZAMS}}/M_{\odot}\lesssim 30$ does not seem sufficient to solve the RSG problem. It is likely that either (i) stars with initial mass in the range $\sim 16-30\,M_{\odot}$ do not produce a bright transient at their death, or (ii) mass loss phenomena other than winds, happening at late stages in the evolution, change the pre-SN appearance of the star. SN observations can constrain the total ejected mass and pre-explosion imaging might also provide constraints (although subject to large uncertainties) on the total mass lost and on the mass loss timing (if the surface chemical composition can be inferred, e.g., Smartt et al., 2009; Gordon et al., 2016).

Potential observational constraints.

A range of special systems, events, and phenomena offer alternative and complimentary ways to constrain mass loss to the more traditional spectral observations of massive stars. These special systems include bow shocks of runaway stars (e.g., Gull & Sofia, 1979; Meyer et al., 2016), SN shocks running into the circumstellar material (e.g., Maeda et al., 2015; Chakraborti et al., 2016; Margutti et al., 2017), flash-spectroscopy of material ejected shortly before core collapse (e.g., Khazov et al., 2016), accretion in wind-fed high mass X-ray binaries, or binary wind collisions. Constraints from special systems can be combined with those arising from observed populations of stars and their compact remnants (including gravitational wave sources). The different mass loss algorithms that we compare here give a range of mass loss timing, which suggests that the chemical composition and dust properties of the circumstellar material may also contain hints regarding the mass loss of massive stars.

Other known unknowns and uncertainties.

In this study, we have only varied the mass loss algorithms used throughout the evolution and their efficiency parameter $\eta$. However, other uncertainties (with possible degeneracies) are known, for example in the treatment of rotation, magnetic fields (e.g., Petit et al., 2017), convective mixing, and overshooting (e.g., Arnett et al., 2015; Arnett, 2015; Arnett & Meakin, 2016; Farmer et al., 2016). The coupling of these uncertainties with the wind mass loss may modify the outcomes of our numerical experiments.

5 Conclusions

Massive star mass loss is a longstanding issue in stellar evolution. Despite decades of observational and theoretical work it remains incompletely understood. Mass loss influences the lifetime and appearance of massive stars, their internal structure at the onset of core collapse, and their total nucleosynthetic yields. Through its effect on the pre-supernova (pre-SN) structure, wind mass loss can impact the outcome of core collapse and the nature of the compact remant, with potential implications for gravitational wave astronomy.

We studied the impact of a broad range of wind mass loss algorithm combinations on the evolution and pre-collapse structure of nonrotating, single, solar-metallicity stars with initial masses of $15,\,20,\,25,\,30,$ and $35\,M_{\odot}$. We compared 12 different mass loss algorithm combinations, drawing from 2 algorithms for the hot phase of the evolution (corresponding roughly to the main sequence), 3 for the cool phase of evolution, and 2 for the Wolf-Rayet phase (if it is reached). We explored the effects of reducing the mass loss rate with an efficiency scaling factor $\eta=1.0,\,0.33$, or $0.1$ to crudely account for reduced stellar mass loss that could be caused, for example, by inhomogeneities in the wind (e.g., clumpiness). The resulting differences in stellar structure and total mass at various stages of the evolution are caused by the different algorithmic representations of stellar winds.

The different mass loss efficiencies and algorithm combinations have profound effects on the evolution and the pre-SN masses of massive stars. On the one hand, this can be expected given the inherent differences of the various mass loss algorithms and the various assumptions that enter them. On the other hand, these algorithms all attempt to describe the same physical process – steady wind mass loss – and less sensitivity to the theoretical/empirical treatment of mass loss would in general be desirable.

We find that the choice of wind efficiency scaling factor $\eta$ has the greatest impact on our stellar models. It affects their total mass loss, their evolutionary path, and their pre-SN structure. $\eta$ is therefore the main uncertainty and limiting factor for our present understanding of wind mass loss and its effects. If the wind efficiency is low, the differences between various mass loss algorithms are less important.

Considering the full range of wind efficiencies $0.1\leq\eta\leq 1.0$, we find that there is a $\sim 50\%$ uncertainty in the pre-SN mass for a given $M_{\mathrm{ZAMS}}$. For fixed efficiency $\eta=1.0$, the uncertainty varies with initial mass and is $\sim 15-30\%$ in most cases. Impulsive mass loss events (eruptions, pulsational instabilities, etc., all neglected here) could only make the uncertainties in the initial to pre-SN mapping more severe.

Despite the large uncertainty in the pre-SN mass, we find that the key observables from pre-SN imaging, the luminosity and effective temperature before explosion, are only mildly affected by varying the mass loss algorithm combination and wind efficiency. The uncertainties in $L$ and $T_{\mathrm{eff}}$ from our models are within observational limits, suggesting that wind mass-loss uncertainties do not affect observational estimates of SN progenitor masses from pre-SN observations for most massive stars (assuming a single star evolution scenario). Nevertheless, the impact of stellar winds on the internal structure of the star can affect the mass and composition of the SN ejecta.

Independent of the employed algorithm during the hot phase, the amount of mass lost in this phase is only a small fraction ($\sim$ few percent) of the total mass. However, since the core can respond directly to it, mass loss during the hot phase creates seed differences that grow during the subsequent evolution, leading to changes in the pre-SN structure and composition of the stellar core. Later, burning shells re-adjust to the mass loss instead of the core itself, and the effect of mass loss is more indirect.

Most of the mass is lost during the late and short cool phase of the evolution. Wind mass loss during this phase is more uncertain. On the cool side of the Hertzsprung-Russel diagram, two different wind driving mechanism might exist: line-driving or dust-driving. The latter, assumed by the van Loon et al. (2005) (vL) algorithm, has a much stronger temperature dependence and produces much higher mass loss rates than algorithms describing line-driven mass loss. vL mass loss can be so strong that it reveals the deep and hot layers of the star. This results in self-damping of the wind itself for higher mass progenitors and thus a higher value of the final mass. The vL algorithm produces very different evolutionary tracks from those obtained with algorithms assuming line driving, with the vL algorithm driving blueward displacements on the Hertzsprung-Russel diagram. The more indirect effect of late mass loss on the pre-SN core structure is difficult to pinpoint with our limited model grid.

In our model grid, only models with $M_{\mathrm{ZAMS}}\geq 35\,M_{\odot}$ and full wind efficiency ($\eta=1.0$) develop into Wolf-Rayet stars. The absence of Wolf-Rayet models with reduced wind efficiency suggests that either (i) our initial mass range is too small to produce Wolf-Rayet stars from single, nonrotating stars, or (ii) other formation channels, such as binarity or impulsive mass loss events, might be dominant. The standard picture for the evolution of massive stars, the so-called “Conti scenario” (e.g., Maeder & Conti, 1994; Lamers, 2013; Smith & Tombleson, 2015), predicts the formation of Wolf-Rayet stars with $M_{\mathrm{ZAMS}}\lesssim 40\,M_{\odot}$, therefore, if mass loss occurs in nature with reduced efficiency ($\eta<1$), then our results with reduced wind efficiency disagree with this prediction.

From our limited sample of six models that we were able to evolve to the pre-SN stage, we find that changing the wind mass loss algorithm combination can lead to changes of the compactness parameter $\xi_{2.5}$ of up to $30\%$ for a given $M_{\mathrm{ZAMS}}$. Moreover, the pre-SN models show a spread in terms of composition profiles, density profiles, and core masses. These uncertainties add to those arising from the incomplete understanding of mixing processes (mainly convection and overshooting). They complicate the study of the core-collapse SN explosion mechanism by adding uncertainty to the initial conditions from which core-collapse SN simulations start. This finding underlines that systematic uncertainty in massive star mass loss can have important implications for the relative number of BHs and NSs resulting from core collapse events, and consequently for gravitational wave sources, and for the nucleosynthetic yields from the explosions of massive stars.

Although the present study provides new insights into the effects of wind mass loss on the evolution and pre-SN structure and appearance of massive stars, it suffers from a number of important limitations that must be addressed by future work.

Our model grid coarsely samples a limited mass range, includes only nonrotating solar-metalicity models, and we could evolve only six models to the onset of core collapse. Extending this grid to finer mass sampling, higher masses, lower metallicity, and evolving all models to the pre-SN stage will be important future work needed to infer more robustly how wind mass loss affects pre-SN structure. Furthermore, we did not consider time-dependent wind efficiency, rotation and magnetic fields, and varitations in mixing processes, which all may have important implications for wind mass loss and its effects on evolution and pre-SN structure.

As in any computational astrophysics study, numerical resolution is a major concern in our work. Initially, the effect of mass loss on the core structure is small, and impossible to resolve with a coarse spatial mesh. We tested our numerical resolution (cf. Appendix B) and ran our calculations at unprecedented spatial resolution (between $20\,000$ and $100\,000$ mesh points) until oxygen depletion to capture the delayed effect of mass loss on the core structure. However, in order to follow the core deleptonization after oxygen depletion with a large nuclear reaction network, we were forced to reduce the spatial resolution. We find that by the end of core oxygen burning, the differences in core structure due to different mass loss algorithm combinations are already pronounced. The limited resolution study that we were able to perform suggests that resolution effects are smaller than the overall effects of mass loss and do not affect our conclusions. Future work is needed to more formally demonstrate robustness and numerical convergence of simulations of the late evolutionary stages.

Acknowledgements.

We wish to thank B. Paxton, F. Timmes, and R. Farmer for invaluable help with MESA. We acknowledge helpful exchanges with W. D. Arnett, M. Cantiello, D. Clausen, A. de Koter, L. Dessart, Y. Götberg, A. L. Piro, J. Fuller, N. Smith, and E. Zapartas. We are also grateful to the referee, G. Meynet, for the helpful comments and suggested improvements. MR thanks the University of Pisa, where this study was initiated as part of his Master Thesis in Physics, and acknowledges the precious support of the “Lorentz group”. This work is supported in part by NSF under grant numbers AST-1205732, AST-1212170, PHY-1151197, and PHY-1125915 and by the Sherman Fairchild Foundation. Some of the computations used resources of NSF’s XSEDE network under allocation TG-PHY100033. Most of the simulation were carried out using the Caltech compute cluster Zwicky funded through NSF grant no. PHY-0960291 and the Sherman Fairchild Foundation. The runs to the onset of core collapse were carried out on the Dutch national e-infrastructure (Cartesius, project number 15162) with the support of the SURF Cooperative.

Appendix A Parametric Wind Algorithms

Here, we give a brief review of the physical assumptions entering in each mass loss algorithm and the corresponding prescribed rate. We summarize in Tab. 9 the scaling of the mass loss rate with physical parameters of the stellar model for each mass loss algorithm considered in this study.

A.1 Vink et al. (V)

The wind mass loss algorithm proposed by Vink et al. (2000, 2001) is based on Monte Carlo simulations of the photon transport in the stellar atmosphere to evaluate the radiative acceleration. Their algorithm is explicitly metallicity dependent and is supposed to be applied only to single OB stars with metallicity $1/30\leq Z/Z_{\odot}\leq 3$ and effective temperature $12\,500\ \mathrm{K}<T_{\mathrm{eff}}<50\,000\ \mathrm{K}$, so during the hot, blue evolutionary phase. In this temperature range, a non-monotonic behavior of $\dot{M}$ as a function of $T_{\mathrm{eff}}$ is expected (the so-called “bi-stability jumps”): normally the lower the temperature, the lower the mass loss rate (because the radiation pressure is proportional to $T^{4}$). However, when the temperature drops below $\sim 25\,000\ \mathrm{K}$, the recombination of $\mathrm{Fe\ IV}\rightarrow\mathrm{Fe\ III}$ provides a new ion with an increased number of lines to drive the wind, resulting in an increased mass loss rate at lower temperatures around $T_{\mathrm{eff}}\sim 25\,000\ \mathrm{K}$. This happens with different iron ion recombinations also at $T_{\mathrm{eff}}\sim 12\,000\ \mathrm{K}$.

Vink et al. (2000, 2001) provide two different formulae for above and below $T_{\mathrm{eff}}\sim 25\,000\ \mathrm{K}$,

[TABLE]

for $27\,500\ \mathrm{K}<T_{\mathrm{eff}}\leq 50\,000\ \mathrm{K}$, and

[TABLE]

for $12500\ \mathrm{K}<T_{\mathrm{eff}}\leq 22500\ \mathrm{K}$. In between these two temperature ranges, it is common practice to simply interpolate between the two formulae (cf. Appendix B). The numbers in parenthesis are the estimates of the error on the last digit reported, according to Vink et al. (2001). These error estimates are usually neglected in stellar evolutionary calculations.

A.2 Kudritzki et al. (K)

The mass loss algorithm proposed by Kudritzki et al. (1989) is based on an analytic solution of the equations for a stationary, isothermal, spherically-symmetric, ideal (neither viscosity nor heat conduction) gas flow with no magnetic fields and no rotation. They include the radiative acceleration $g_{\mathrm{ph}}$ in the momentum equation using the standard parametrization of Castor et al. (1975) (see also Eq. 16 in Pauldrach et al. 1986),

[TABLE]

where $g_{\mathrm{rad}}^{\mathrm{Th}}$ is the radiative acceleration due to Thomson scattering, and the second term in the brackets is the so-called “force multiplier”, i.e. the line acceleration in units of $g_{\mathrm{rad}}^{\mathrm{Th}}$. It depends on the Thomson cross section $\sigma_{\mathrm{Th}}$, the speed of sound $c_{s}$, the electron number density $n_{e}$, and three free parameters $k$, $\alpha$, and $\delta$. $k$ can be interpreted roughly as the number of lines strong enough to have an effect, and $\alpha$ as the slope of the distribution of the number of lines as a function of their strength. The parameter $\delta$ and the correction factor $\mathrm{CF}$, which is the ratio between the opacity as a function of the incoming angle and the opacity in the radial direction, are used to include the “finite cone-angle effect” to account for photons travelling in non-radial directions (e.g., Castor et al. 1975). If $\delta$ and $\mathrm{CF}$ were both equal to one, the parametrization would be valid only in the “radial streaming” limit (Abbott 1982; Pauldrach et al. 1986), i.e. considering only incoming photons from the radial direction. While this is a good approximation in the outer portion of the wind, where $r\gg R$ and $R$ is the stellar radius, it is quite poor in the inner portion, where the mass loss rate is determined and photons traveling in non-radial directions can have a relevant effect. For completeness, we report the expression for the correction factor $\mathrm{CF}$ as a function of the dimensionless radial coordinate $x\stackrel{{\scriptstyle\mathrm{def}}}{{=}}r/R$ and $h\stackrel{{\scriptstyle\mathrm{def}}}{{=}}d\log_{10}(x)/d\log_{10}(v)$, cf. Eq. (4) in Kudritzki et al. (1989),

[TABLE]

We refer the reader to Castor et al. (1975); Pauldrach et al. (1986); Vink (2015) and references therein for more details on this parametrization.

To find an analytic solution to their model, further assumptions are needed. Kudritzki et al. (1989) impose a “$\beta$-law” velocity field (a common assumption found to be close to numerical solutions, see Lamers 2013),

[TABLE]

where $\beta$ is a free parameter assumed to be $\beta=1$, $v_{\infty}$ is the asymptotic velocity of the wind and $R$ is the stellar radius. Note that this assumption means Kudritzki et al. (1989) do not solve for the dynamics of the system, but rather assume the velocity structure and solve self-consistently for the density and the acceleration. In the limit where $v\gg c_{s}$ (reasonable in the outer portion of the wind), Kudritzki et al. (1989) find a solution of the form

[TABLE]

where $v_{\mathrm{th}}$ is the thermal velocity of protons, $k$, $\alpha$, and $\delta$ are the free parameters, $L$ is the luminosity, $v_{\mathrm{th}}$ is the thermal velocity, $\Gamma=L/L_{\mathrm{Edd}}$ is the Eddington ratio, and $\tilde{D}$ is a function of the free parameters and the velocity (cf. Eqs. 47, 62, 65 in Kudritzki et al. 1989).

The main limitation of Eq. 9 is that the parameters $k$, $\alpha$ and $\delta$ are not constants, but rather depend on the optical depth. The numerical values commonly adopted are $\alpha=0.657$, $\beta=1$, $\delta=0.095$, and $k=0.085$ and they are calibrated on $\zeta$ Puppis by Pauldrach et al. (1994). These should be interpreted as values averaged over the optical depth.

A.3 de Jager et al. (dJ)

The wind algorithm proposed by de Jager et al. (1988) is an empirical relationship of the form $\dot{M}\equiv\dot{M}(T_{\mathrm{eff}},L)$. This choice of variables uses only observable quantities, making it easy to track the mass loss rate while the star moves on the HR diagram. This allows a better understanding of how mass loss changes during stellar evolution. The drawback is that no information about the physical mechanism driving the wind is considered.

To formulate a reliable mass loss algorithm, de Jager et al. (1988) collect from the literature mass loss rates observed with different techniques for a sample of galactic stars with spectral types from O to M. They determine the “average” measured mass loss rate for stars observed with multiple techniques and the deviation from this average for each available measurement and for each star. Then, they define the “average intrinsic error per determination” as the one-sigma value of the distribution of these deviations. They fit their entire data sample with a sum of Chebychev polynomials of the first kind $T_{n}(x)=\cos(n\ \mathrm{arccos}(x))$:

[TABLE]

What is commonly used in stellar evolution codes is the first-order approximation to Eq. 10,

[TABLE]

To asses the quality of this fit, they derive the distribution of the differences between the observed values of $\log_{10}(-\dot{M})$ with the result of Eq. 11. The standard deviation of this distribution is $\sim 0.45$, slightly larger than the ”averaged intrinsic error per determination” previously determined. This indicates that quantities other than $T_{\mathrm{eff}}$ and $L$ must physically enter in the determination of $\dot{M}$, and the parametrization in Eq. 11 is incomplete. The main limitation of the dJ algorithm is that it is representative of the “averaged statistical behavior” of stellar winds in the entire HR diagram, which might average over different physical regimes. It is important to note that WR and Be stars are intentionally excluded from the data sample used to derive Eq. 11.

A.4 Nieuwenhuijzen & de Jager (NJ)

The algorithm of Nieuwenhuijzen & de Jager (1990) is intended as an improvement over the de Jager et al. (1988) algorithm of Eq. 11, since it is derived from the same data sample with a similar method but includes the dependence of mass loss on the total stellar mass. The goal of this is to capture one of the missing stellar parameters entering in the mass loss determination indicated by the large standard deviation of Eq. 11. The Nieuwenhuijzen & de Jager (1990) algorithm also translates the temperature dependence into a radius dependence.

Since the total mass is not a directly observable quantity for single stars, the authors’ mass determination is based on stellar model calculations. The theoretical models used are from Maeder & Meynet (1988, 1989). However, different stellar evolution codes consider a large variety of physical processes (e.g. for mixing, mass loss, etc.), or just use different implementations of them. Hence, there is spread in the stellar masses found at the same point of an evolutionary track. This implies that the Nieuwenhuijzen & de Jager (1990) mass loss algorithm depends on the set of stellar models used to derive it. This drawback applies to all mass loss algorithms involving a functional dependence of the form $\dot{M}\equiv\dot{M}(M)$ derived from stellar evolution calculations. However, Mauron & Josselin (2011) suggests that the dependence on the total mass of the Nieuwenhuijzen & de Jager (1990) mass loss rate is so weak (see Eq. 14) that it can be averaged (by substituting $M$, which changes during the evolution, with a constant value) without dramatic consequences on the evolved stellar model.

Stars with different masses pass through the same point on the HR diagram at different stages of their evolution. Therefore, in order to include the total mass as a variable for the mass loss rate, Nieuwenhuijzen & de Jager (1990) determine an “average expected mass” $\bar{M}$ of a star at a given ($T_{\mathrm{eff}},L$) point. This value is derived as follows. The authors define a dwell time representing the time for a star to travel over a unit length track on the HR diagram, i.e. ,

[TABLE]

where $\delta t$ is the time spent to travel over the $(\delta\log_{10}(T_{\mathrm{eff}}/\mathrm{[K]}),\delta\log_{10}(L/L_{\odot}))$ distance.

For every point on the HR diagram, there are N left-ward or right-ward subtracks of the stellar evolutionary tracks crossing it. Let $t^{(d)}_{n}$ be the dwell time for the n-th subtrack. The average expected mass $\bar{M}$ is obtained as:

[TABLE]

where $\Psi$ is the initial mass function for stars on the subtrack considered, and $dM_{n}/d\log_{10}(L/L_{\odot})$ is the density of tracks over a unit $\log_{10}(L/L_{\odot})$ interval.

The authors perform a fit of the data set555Which excludes WR and Be stars. used in de Jager et al. (1988), adding the value of $\bar{M}$ from Eq. 13 to the set, and they find the interpolation formula666The formula in the abstract of Nieuwenhuijzen & de Jager (1990) has a typo, see their Eq. 2. (where $\bar{M}$ is substitued by the total mass $M$ to use the formula in a stellar evolution simulation)

[TABLE]

With the inclusion of $M$ among the parameters determining $\dot{M}$, Nieuwenhuijzen & de Jager (1990) obtain a standard deviation of $\sim 0.37$ for the distribution of the differences between the prediction of Eq. 14 and the observed values, comparable to the standard deviation of the distribution of the differences between individual mass loss determination used as input data. The main limitations of this algorithm are its dependence on the input stellar models, and its statistical-average nature (in the same sense as de Jager et al. 1988).

A.5 Van Loon et al. (vL)

The mass loss rate of van Loon et al. (2005) is empirically determined on the basis of observations of a sample of oxygen-rich AGB and RSG stars in the Large Magellanic Cloud (LMC). The van Loon et al. (2005) analysis is based on a dust-driven wind model: AGB and RSG stars have very extended and cool envelopes where dust grains might form through sublimation. The photons from the radiation field transfer momentum to these grains, pushing them away. The dust grains drag the gas with them through collisional coupling. To obtain their mass loss algorithm, van Loon et al. (2005) fit the observed IR spectra to synthetic spectra obtained with a simple model of the gas/dust mixture (identical grains and dust-to-gas ratio set to the value observed at $Z_{\odot}$ rescaled to $Z_{\mathrm{LMC}}$), using $T_{\mathrm{eff}}$ and $L$ as variables. They obtain the relationship

[TABLE]

where the numbers in parentheses indicate the estimated errors on the last digits. These errors are typically neglected in stellar evolution codes. The main limitation of this algorithm is the high uncertainty of the dust grain properties (mass fraction, opacity, when they form, etc.). Since it gives such a high mass loss rate, the van Loon et al. (2005) algorithm is sometimes combined with a high efficiency factor to ad hoc mimic non-wind mass loss phenomena (e.g., Meynet et al. 2015).

A.6 Nugis & Lamers (NL)

The mass loss algorithm derived by Nugis & Lamers (2000) applies only to WR stars. The wind mass loss rate of these stars depends strongly on their chemical composition: not only the metallicity has an important role, but also the helium mass fraction $Y$. This is because the amount of helium in the stellar atmosphere influences its temperature and therefore the ionization fraction and the level populations of all other atoms and ions. Nugis & Lamers (2000) derive a mass loss rate algorithm as a function of the luminosity and the chemical composition starting from a relevant sample of observed galactic WR stars. Their sample is made of two subsets of stars: one in which both mass and distance (i.e. luminosity) are known, thanks to binarity and membership association in open clusters; and another subset for which the intrinsic luminosity is not known. They use stars from the first subset to derive an empirical bolometric correction777The bolometric correction is the difference between the bolometric magnitude and the observed (visual) magnitude, influenced by the instrumental intrinsic band pass, $\mathrm{BC}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}M-M_{\mathrm{obs}}$. They then use a theoretical mass-luminosity relation to infer the luminosity of stars in the second subset, and correct it with the previously derived bolometric correction. Note that the mass-luminosity relation is determined using as input the age of the star and its spectral type, not its luminosity, therefore this relation can be consistently used to estimate the luminosity of stars in the second subset (Nugis & Lamers 2000).

The mass loss rate observed for the stars in the sample are then fitted as follows. The authors make two independent fits for stars of different composition and then merge them together in a single formula, valid for all WR stars:

[TABLE]

where the numbers in parentheses are the estimates of the error on the last digits reported. According to Nugis & Lamers (2000), the mass loss algorithm for WR stars cannot be expressed as a function of $T_{\mathrm{eff}}$ and/or the radius of the star $R$. This is because WR winds are so strong (i.e. dense) that they are optically thick, and thus the observed radius is a function of the wavelength: the black body relation $L=4\pi R^{2}\sigma T_{\mathrm{eff}}^{4}$, where $\sigma$ is the Stefan-Boltzmann constant, loses its meaning.

A.7 Hamann et al. (H)

This wind scheme, which applies only to WR stars, is a combination of the algorithms from Hamann et al. (1982), Hamann et al. (1995), and Hamann & Koesterke (1998). It is derived from a spherically-symmetric, homogeneous, and stationary (but not static, i.e. $\partial_{t}=0$ but $v\neq 0$) expanding WR atmosphere model. The authors assume an ad-hoc velocity structure $v\equiv v(r)$ as follows. For the supersonic part of the wind, they assume a $\beta$-law of the form of Eq. 8 with $\beta=1$, while for the subsonic part, $v(r)$ is chosen in such way that the density approaches the hydrostatic limit. These assumptions are in reasonable agreement with the numerical solutions. Note that since the velocity field is imposed, the acceleration is not computed. This allows the authors to adopt a very simple chemical composition, since they do not need to evaluate the line-driven acceleration and do not need to keep track of all atomic/ionic species and their level populations. The authors include only ions of H and He, and the radiation field is considered only to determine non-LTE populations of these species. The temperature stratification is derived with the assumption of a gray LTE model, assuming a value of $T_{\mathrm{eff}}$ at the base of the atmosphere determined by the stellar luminosity $L$ and the radius $R$ via a black body relation $L=4\pi R^{2}\sigma T_{\mathrm{eff}}^{4}$, in contrast to the suggestion of Nugis & Lamers (2000). A synthetic spectrum is derived from the simulations and a best fit to the observed line profiles for the ions of H and He is obtained via variation of the stellar parameters (i.e. the radius of the inner boundary of the atmosphere $R$ and the luminosity $L$, the surface hydrogen mass fraction $X_{s}$ and $\dot{M}$, Hamann et al. 1995). Once the stellar parameters are known from this fit, a mass loss formula is derived for high luminosity WR stars, i.e. $\log_{10}(L/L_{\odot})>4.5$. The algorithm for the low luminosity WR stars, i.e. $\log_{10}(L/L_{\odot})<4.5$, is derived with a similar technique in Hamann et al. (1982), but the spectra fitted are from a small sample of Helium stars (i.e. stars undergoing He shell burning with most of the mass in a CO core, without H lines in their spectra). The resulting formula is

[TABLE]

Improvements that take into account inhomogeneities in the wind are suggested in Hamann & Koesterke (1998). Specifically, the authors suggest to reduce the wind efficiency by a factor between 2 and 3 to account for the wind clumpiness, which strongly affects the fitted spectral lines.

Appendix B MESA input physics, customization, and resolution study

Here, we provide and discuss the MESA parameters not directly related to the physics of mass loss that we omit in Sec. 2 for brevity.

We employ the Ledoux criterion (Ledoux 1947) to determine convective stability. This means that we consider the effects of the temperature and chemical composition gradients on the stability of a region, and allow for semiconvective mixing driven by compositional gradients. We follow the suggestions of Sukhbold & Woosley (2014) for the treatment of convective mixing (relying on mixing length theory; Böhm-Vitense 1958) and overshooting, but we also include thermohaline mixing using the default MESA treatment of this process (Kippenhahn et al. 1980). Our specific MESA parameter choices are as follows: mixing length parameter $\alpha_{\mathrm{mlt}}=2$, exponential overshooting/undershooting parameters $f_{\mathrm{ov}}=0.025$ and $f_{0}=0.05$ for all regions (see Moravveji et al. 2016; Paxton et al. 2011), semiconvection efficiency of $0.2$, and thermohaline mixing with coefficient $2$ according to Kippenhahn et al. (1980).

We employ the 45-isotope nuclear reaction network mesa_45.net (in $MESA_DIR/data/net_data/nets/) for all models until oxygen depletion. We switch to a customized 203-isotope network for models evolved further, and initialize to zero the abundances of all the new isotopes introduced (i.e. adjust_abundances_for_new_isos = .false.). We obtain the 203-isotope network by taking the union of the isotope sets of mesa_45.net and mesa_201.net. The list of isotopes included in the network is available at https://stellarcollapse.org/renzo2017.

Stars more massive than $\sim 20\,M_{\odot}$ can develop a radiation-dominated, near super-Eddington convective envelope (e.g., because of the iron opacity bump) in which the convective flux is not sufficient to transport energy outward. This can lead to numerical (and possibly physical, e.g., Langer 1997) density and pressure inversions and consequently dynamical instabilities (see, e.g., Joss et al. 1973; Paxton et al. 2013, and references therein). To handle these envelope issues, MESA introduces the so-called MLT++ scheme (see Paxton et al. 2013). MLT++ arbitrarily decreases the superadiabaticity in near-to-Eddington ($L\geq 0.5L_{\mathrm{Edd}}$) radiation dominated ($P_{\mathrm{gas}}/P_{\mathrm{tot}}\leq 0.35$) convective envelopes in order to prevent numerical issues. We note that MLT++ can modify the radius and luminosity of the star by artificially enhancing the energy flux carried by convection, and thus indirectly, it can also change the mass loss rate. Furthermore, we find that MLT++ combined with the simple_atmosphere option for the outer boundary condition causes an unphysical oscillation of the mass loss rate and surface characteristics (e.g., luminosity $L$, effective temperature $T_{\mathrm{eff}}$, and radius $R$; Renzo 2015). To avoid this issue, we employ the more realistic atmosphere boundary condition provided by the Eddington_grey option in all our models. This uses the Eddington gray $T(\tau)$ relation,

[TABLE]

to find the boundary pressure on the photosphere, instead of guessing the approximate photosphere location.

At the beginning of each evolutionary step, MESA evaluates the mass loss rate $\dot{M}$ according to the algorithm specified by the user, and then removes from the surface the amount of mass $\dot{M}\times\Delta t$, where $\Delta t$ is the timestep. The mass loss algorithm to evaluate $\dot{M}$ can be chosen from the many built-in algorithms, or it can be implemented by the user using the run_star_extras.f hooks to override the default MESA routines888See also http://mesa.sourceforge.net/run_star_extras.html. We use the latter option to combine three different algorithms (one for the hot phase of the evolution, one for the cool phase, and one for the WR phase, see Sec. 2.2). When possible, we call the built-in MESA mass loss routines from our run_star_extras.f, except for the V algorithm. The default implementation of the V mass loss algorithm in MESA uses a linear interpolation between Eq. 4 and Eq. 5. This results in a jump of $\dot{M}$ (as a function of time $t$) with discontinuous derivative (i.e. $M(t)$ is not $\mathcal{C}^{2}$) as $T_{\mathrm{eff}}$ drops below $25\,000\ \mathrm{K}$ during the evolution. For a physically more realistic, smoother time dependence of $\dot{M}$, we implement a hyperbolic tangent interpolation between the two formulae. This is depicted in Fig. 10.

To save MESA “photos”999These are binary files that allow one to restart a run and obtain bit-to-bit identical results (provided that the parameter set does not change). See also http://mesa.sourceforge.net/. and to call the built-in implementation of the K algorithm from run_star_extras.f, we copy the Fortran modules write_model.mod and kuma.mod from $MESA_DIR/star/make into$MESA_DIR/include in our standard MESA installation. This is necessary since these two modules do not have a “public” interface (Paxton et al. 2011) in MESA release version 7624.

We describe the settings for spatial and temporal resultion in the next Sec. B.1. We use the default MESA settings for massive stars (see $MESA_DIR/star/inlist_massive_defaults) for anything else not explicitly mentioned.

B.1 Resolution Dependence

Any computational study in astrophysics must carefully assess the effects of numerical resolution (in both space and time) on its results. To study the sensitivity of our results to variations in temporal and spatial resolution, we carry out a resolution study using a 30$M_{\odot}$ star since stars of this $M_{\mathrm{ZAMS}}$ appeared to be the most sensitive to the resolution in our preliminary calculations using the default MESA parameters. The post core carbon depletion evolutionary tracks produced by the default MESA parameters show large amplitude oscillations (e.g., in the central temperature–central density plane) when varying the spatial discretization. It is unlikely that nature would do this. Thus, such oscillations are most likely artificial, and they are generally worse at higher initial mass. The aim of this section is to find a set of parameters that reduces and possibly eliminates these oscillations.

For each set of resolution parameters, we run our test model to oxygen depletion with two different mesh refinement parameters (mesh_delta_coeff=1.0 and mesh_delta_coeff=0.5), and two different wind mass loss algorithm combinations (V-dJ and K-NJ; both with $\eta=1.0$). The use of two different mass loss algorithm combinations allows us to check that our settings are not cherry-picked for a particular model in our grid. However, the results of our study indicate that the resolution dependence is insensitive to the mass loss algorithm combination.

To obtain numerically converged results, we use both the available MESA controls (specified in the inlists) and customized routines in run_star_extras.f. Both are available at https://stellarcollapse.org/renzo2017.

B.1.1 Timestep selection throughout the evolution

To avoid overstepping relevant physical processes101010Preliminary calculations with low resolution showed that in some rare cases, e.g. during the Herzsprung gap, MESA would try to overstep the Kelvin-Helmholtz timescale., we enforce a customized timestep control for the evolution of our models (see the routine extras_finish_step in our run_star_extras.f). In particular, we tighten the default MESA (release version 7624) timestep controls by manually enforcing

[TABLE]

where $\Delta t_{n+1}$ is the timestep proposed at the end of the $n-\mathrm{th}$ step for the next ($(n+1)-\mathrm{th}$) step, and all the timescales on the right hand side refer to the $n-\mathrm{th}$ step. On the right hand side of condition 19,

[TABLE]

are the Kelvin-Helmholtz timescale and the mass change timescale, respectively.

Note that the right-hand side of condition 19 is evaluated with the quantities of the $n-\mathrm{th}$ step, but it limits the $(n+1)-\mathrm{th}$ step. Moreover, until oxygen depletion, we limit the timestep by setting varcontrol_target = 1.0d-4 as the maximum relative variation for quantities in each computational cell. We also limit the amount of matter nuclearly processed in one single timestep by setting the maximum variation of the mass fraction of each element due to nuclear buring in each time step to dX_nuc_drop_limit = 1.0d-4.

B.1.2 Spatial meshing until oxygen depletion

Until oxygen depletion, we impose a maximum value for the fraction of the total mass in each computational cell by setting max_dq = 0.5d-4 in the inlist. This means that each of our models has at least $1/\texttt{max\_dq}=2\times 10^{4}$ computational cells. Moreover, we refine the mesh around some specific regions of the star to focus resolution there:

• $\bullet$

Regions with steep temperature gradients: to better resolve the deep interior of the star, we limit the maximum variation of $d\log_{10}(T)/dm$ across adjacent cells, where $T$ is the temperature and $m$ is the mass coordinate. Specifically, we impose a maximum variation of mesh_delta_coeff$/10$ for $\log_{10}\left(d\log_{10}(T)/dm+1\right)$. This is is achieved using the other_mesh_fcns_data routine in run_star_extras.f;

• $\bullet$

Stellar surface: to properly resolve the amount of mass lost at each timestep, we impose that the outermost $0.5\,M_{\odot}$ of the star is sampled by at least $500\times(\texttt{mesh\_delta\_coeff})^{-1}$ cells (see other_mesh_fcns_data routine in run_star_extras.f for the implementation). Note that for each timestep $\Delta t$, $\dot{M}\times\Delta t\ll 0.5M_{\odot}$;

• $\bullet$

Edges of burning regions: to resolve the edges of burning regions, we constrain the spatial variation of $d\log_{10}(\varepsilon_{\mathrm{nuc}})/d\log_{10}(P)$ by multiplying its maximum allowed variation (regulated by mesh_delta_coeff) by 0.015 (as in Dessart et al. 2013, private communication). This is done separately for each nuclear burning process, see mesh_dlog_*_dlogP_extra in the inlists;

• $\bullet$

Edges of the cores of different composition: to resolve sharp variations in the chemical composition, we impose mesh_delta_coeff$/20$ as the maximum spatial variation allowed for the mass fractions of several isotopes (${}^{1}\mathrm{H}$, ${}^{4}\mathrm{He}$, ${}^{12}\mathrm{C}$, ${}^{16}\mathrm{O}$, ${}^{20}\mathrm{Ne}$, ${}^{28}\mathrm{Si}$,${}^{24}\mathrm{Mg}$, ${}^{32}\mathrm{S}$, ${}^{54}\mathrm{Fe}$, ${}^{56}\mathrm{Fe}$). These are specified via xa_function_* in the inlists.

Our models with these settings and mesh_delta_coeff = 1.0 have typically between $\sim 50\,000$ and $\sim 100\,000$ mesh points. Figure 11 illustrates for the model using the V-dJ-NL mass loss algorithm combination and $\eta=1.0$ that our setup does not produce any appreciable variations when changing mesh_delta_coeff by a factor of 2 (corresponding to an increase of the number of spatial mesh points from 54 814 to 77 400 at oxygen depletion). Note the linear scale for the central temperature – central density evolutionary tracks.

B.1.3 Re-meshing after Oxygen Depletion

After oxygen depletion, we switch from a 45-isotope to a 203-isotope nuclear reaction network (cf. Sec. B). This switch forces us to decrease the number of spatial mesh points for the following reason: MESA solves the fully coupled set of equations for stellar structure and evolution Paxton et al. (2011). This results in the use of a work array of length $\mathcal{L}$ which scales as

[TABLE]

where $N_{\mathrm{iso}}$ is number of isotopes in the nuclear reaction network and $N_{\mathrm{z}}$ is the number of mesh points (see the routine get_newton_work_sizes in $MESA_DIR/star/private/star_newton.f90).$\mathcal{L}$is stored as a 4-byte integer in MESA, which sets the maximum adressable memory to$\sim 17$GB (R. Farmer, private communication). Therefore, for a nuclear reaction network including$\sim 200$isotopes, the maximum number of mesh points cannot exceed$\sim 17,000$. Changing the relevant variables to an 8-byte integer data type would require substantial changes throughout MESA, which we have opted to defer to future work.

Another significant limitation on the resolution that can be achieved in the very late evolutionary phases is the stability of stellar evolution codes: many highly uncertain and/or poorly understood physical phenomena take place in the cores of evolved massive stars, and the evolutionary timescale gets progressively closer to the (neutrino-)thermal timescale and, finally, the dynamical timescale.

Experiments with the highest achievable resolution resulted in frequent failures of the code to find solutions to the stellar structure equations. However, note that by the time core oxygen burning is complete, the variations in core structure caused by different mass loss algorithm combinations are already largely developed and will be amplified by the subsequent evolution.

We reduce the number of mesh points in our models down to a few thousand (the precise value varies from model to model). We do this by restarting from oxygen depletion and running MESA’s re-meshing algorithm for 100 timesteps of $\Delta t<10^{-9}\,\mathrm{s}$ with nuclear burning and neutrino cooling turned off (see inlist_remesh). At the same time, we shut down thermohaline mixing, which is a secular process that does not have time to produce any physically relevant change in the star after core oxygen burning. We also drop the spatial mesh refinement criteria described above for the remaining evolution. After the 100 re-meshing timesteps during which the star is de facto frozen in its state, we resume the evolution by again turning on nuclear burning and neutrino cooling and not imposing a maximum timestep, but we do not re-enable thermohaline mixing. We compare pre-re-meshing and post-re-meshing core structure, thermodynamics, and compactness parameter and ensure that the post-re-meshing core is still very well resolved with the reduced resolution as we proceed in the evolution toward core collapse.

As the temperature in the core increases, the nuclear burning rate accelerates. When the central temperature rises above $3\times 10^{9}\,\mathrm{K}$, we progressively relax the constraints on the timestep from nuclear burning. We do this by modifying, at the beginning of each timestep, the parameters loaded in our inlist using the routine extras_startup of our run_star_extras.f. Specifically, for $3.0\leq T_{c}/\mathrm{[10^{9}K]}\leq 3.5$ we impose dX_nuc_drop = 5d-3, and for even higher $T_{c}$ we only require dX_nuc_drop = 5d-2 (cf. Sec. B.1.1). After silicon core depletion ($X_{c}(^{28}\mathrm{Si})\leq 0.001$), we also decrease the spatial resolution of the innermost infalling core by increasing mesh_delta_coeff_for_highT=3.0 (used where $\log_{10}(T/\mathrm{[K]})\geq 9.3$, the value used during the previous evolution is 1.0).

Finally, an important question to address is the sensitivity to numerical resolution of the subsequent evolution toward core collapse and of the final pre-collapse structure and compactness parameter. Due to the memory limitations and computational difficulties described in the above, we are unable to carry out a rigorous convergence test with MESA at this time. However, in order to gain some insight into the effects of resolution, we choose the $25\,M_{\odot}$ V-vL $\eta=0.33$ model and carry out two additional simulations at reduced resolution from oxygen depletion to the onset of core collapse using our 203-isotope nuclear reaction network. Specifically, we choose (mesh_delta_coeff, mesh_delta_coeff_for_highT) = (1.75, 2.75) and (2, 3), whereas the standard setting for our models discussed in Sec. 3.7 is (1, 1). Figure 12 shows the evolution of the compactness parameter from oxygen depletion (after re-meshing) to the onset of core collapse for these models. We find that these lower-resolution models evolve qualitatively very similar to our standard model, but produce variations in the pre-SN compactness parameter of about $\sim$ 9%, which is smaller than the variations of $\sim 30\%$ caused by different wind mass loss algorithm combinations.