# The influence of transition metal solutes on dislocation core structure and
values of Peierls stress and barrier in tungsten.

G. D. Samolyuk  Materials Science and Technology Division, Oak Ridge National
Laboratory, Oak Ridge, TN 37831, USA  Y. N. Osetsky  Materials Science and
Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
R. E. Stoller  Materials Science and Technology Division, Oak Ridge National
Laboratory, Oak Ridge, TN 37831, USA

(4 February 2009)

###### Abstract

Several transition metals were examined to evaluate their potential for
improving the ductility of tungsten. The dislocation core structure and
Peierls stress and barrier of $1/2\langle 111\rangle$ screw dislocations in
binary tungsten-transition metal alloys (W1-xTMx) were investigated using
first principles electronic structure calculations. The periodic quadrupole
approach was applied to model the structure of $1/2\langle 111\rangle$
dislocation. Alloying with transition metals was modeled using the virtual
crystal approximation and the applicability of this approach was assessed by
calculating the equilibrium lattice parameter and elastic constants of the
tungsten alloys. Reasonable agreement was obtained with experimental data and
with results obtained from the conventional supercell approach. Increasing the
concentration of a transition metal from the VIIIA group, i.e. the elements in
columns headed by Fe, Co and Ni, leads to reduction of the $C^{\prime}$
elastic constant and increase of elastic anisotropy A=$C_{44}/C^{\prime}$.
Alloying W with a group VIIIA transition metal changes the structure of the
dislocation core from symmetric to asymmetric, similar to results obtained for
W1-xRex alloys in the earlier work of Romaner et al (Phys. Rev. Lett. 104,
195503 (2010)). In addition to a change in the core symmetry, the values of
the Peierls stress and barrier are reduced. The latter effect could lead to
increased ductility in a tungsten-based alloy. Our results demonstrate that
alloying with any of the transition metals from the VIIIA group should have
similar effect as alloying with Re.

###### pacs:

74.70.Dd,72.15.-v,74.25.-q

##  I Introduction

Tungsten is the prime candidate for use in the divertor of future fusion
reactors because of its high melting temperature and resistance to sputtering
El-Guebaly et al. (2011); Nygren et al. (2011). However, its lack of ductility
is an impediment to its use. The low-temperature brittleness is a common
problem for all metals from the VIA group, such as chromium, molybdenum and
tungsten Klopp (1969). The ductility of these metals can be improved by
alloying with rhenium, leading to the so-called ”Re effect” Geach and Hughes
(1956); Klopp (1969); Savitskii et al. (1965). However, Re, is a very rare and
expensive element. It is therefore desirable to find alternate elements which
provide a similar increase in ductility at lower cost. An experimental
investigation of the range of candidates in the periodic table would be rather
expensive, but computational materials science methods based on accurate
first-principles calculations provide a very promising way to narrow the range
of possible candidates.

Several possible mechanisms for the Re-effect have been discussed Edington et
al. (1966); Klopp (1969); Raffo (1969); Luo et al. (1991); Gornostyrev et al.
(1991); Trefilov et al. (1975); Kurdyumova et al. (1980), and the following
two mechanisms are selected as most promising in application to monocrystals,
i) solid solution softening in which an impurity improves mobility of 1/2(111)
screw dislocation and ii) enhancement of cross-slip in which an impurity
modifies the dislocation core structure making it easier to cross-slip and
increasing the number of possible slip planes. Both of these are related to
the effect of impurities on the dislocation core and may be amenable to
investigation by first-principles computational methods even though such
calculations cannot be used to directly estimate the mechanical properties.
Therefore, it should be possible to define one or more calculable figures of
merit that are related to a material’s elastic properties and potentially to
ductility. Among the potential figures of merit are the material’s individual
elastic constants, Poisson’s ratio Gao et al. (2008), and Peierls stress; and
there is strong evidence that these parameters are also directly related to
the electronic structure of particular impurities.

The experimentally observed strong correlation between the number of valence
electrons on the solute atom and the degree of softening Hiraoka et al.
(2004); Klopp (1975); Stephens and Witzke (1975) points to the importance of
electronic factors in solution softening Medvedeva et al. (2005, 2007). First
principles calculation Romaner et al. (2010); Li et al. (2012) indicate that
the ductilizing effect of Re could originate from ”direct” improvement of the
mobility of $1/2\langle 111\rangle$ screw dislocations Raffo (1969); Trinkle
and Woodward (2005); Medvedeva et al. (2005, 2007), either by decreasing the
Peierls stress, or enhancing cross-slip by changing the slip plane from
$\\{110\\}$ to $\\{112\\}$ which would increase the number of available slip
planes from 6 to 12 Edington et al. (1966); Garfinkle (1966). First principles
density functional (DFT) calculations demonstrated that alloying with Re in a
W1-xRex alloy leads to a transition of the $1/2\langle 111\rangle$ screw
dislocation structure from the symmetric core to asymmetric core, and to a
reduction in Peierls stress Romaner et al. (2010).  Closely related results
were obtained from first-principles calculation for Mo alloys with 5$d$
transition metals Trinkle and Woodward (2005). The authors Trinkle and
Woodward (2005) placed solutes in a row along the $1/2\langle 111\rangle$
dislocation core and calculated the change of stiffness associated with moving
the row along the $\langle 111\rangle$ direction. According to the results,
solutes having fewer $d$ electrons (Hf and Ta) increase the stiffness, which
authors infer strengthens the Mo alloy, whereas those having more $d$
electrons (Re, Os and Ir) decrease stiffness, leading to softening.

The goal of the present investigation is to study possible substitutes for Re
in alloys that will result in a similar ductilizing effect as Re. The
properties of W-transition metal (TM) alloys, were modeled using the virtual
crystal approximation Faulkner (1982) (VCA).  The applicability of this
approach to modeling elastic properties, structural stability and phonon
properties of a W-Re alloy has been demonstrated Gornostyrev et al. (1991);
Persson et al. (1999); Ekman et al. (2000). The approach we use can be
described as follows. First, a dislocation is placed in the effective media
representing the W-TM alloy. The properties of the dislocation in this medium
are calculated exactly but for a W-TM alloy described by a pseudopotentail
which is the weights average of the actual pseudopotentials of the impurity
and host atoms and the averaged number of electrons. A periodic lattice is
used with the lattice parameter chosen by obtain zero pressure condition. This
approach corresponds to a zero order contribution to the electronic system
energy expansion with respect to the difference between real atomic potentials
and the virtual atom. Such an approach allows separating the so-called ”band
structure effects”, in this particular case filling of $d$-states, from the
effect of local modification of the lattice due to W substitution. The latter
is not possible in a direct super cell calculation. It was demonstrated
Romaner et al. (2010) that dislocation properties are sensitive to
modifications of this effective media and alloying with Re leads to a sizable
reduction of Peierls stress and barrier. Within this approach, a comparison of
elastic constants calculated using the VCA and super-cell methods provides a
verification of the VCA accuracy as we demonstrate below.

We show that alloying with transition metals with a higher number of $d$
electrons (VIIIA group) reduces the Peierls stress and barrier for a
$1/2\langle 111\rangle$ screw dislocation. Similar to W-Re alloys, the
dislocation core symmetry is reduced by alloying with TM from group VIIIA. It
is demonstrated that the scale of Peierls barrier reductions are similar for
all W1-xTMx alloys with the same ratio of electrons per atom ($e/a$).
Although, the VCA approach describes an influence of band-structure on the
properties of screw dislocations, it does not describe the discrete nature of
alloy structure. However, this approach permits an assessment that can
significantly reduce the range of possible solute candidates suitable for Re
substitution. This paper is organized as follows. In section II we briefly
review computational details used in the calculations. Section III.1 describes
the electronic structure of tungsten and its alloys with transition metals
from IVA, VA and VIIA groups. This is an almost complete list of possible
binary solid solutions Lassner and Schubert (1999). The results of
modifications to the elastic constants with alloying are discussed in section
III.2. Section III.3 presents our prediction of core structure, Peierls
barrier and Peierls stress for W1-xTMx alloys. Finally, the conclusions are
given in section IV.

##  II Computational approaches

The electronic structure within the generalized gradient approximation (GGA)
of density functional theory (DFT) was calculated using the QUANTUM ESPRESSO
(QE) package Giannozzi et al. (2009). The calculation was done using a plane-
wave basis set and ultrasoft pseudo-potentials optimized in the RRKJ scheme
Rappe et al. (1990). We used the Perdew-Wang Perdew and Wang (1992) exchange-
correlation functional. The Brillouin zone (BZ) summations were carried out
over a $24\times 24\times 24$ BZ grid for the system with one unit cell and
$16\times 16\times 16$ grid for the supercell containing $2\times 2\times 2$
unit cells, with electronic smearing with a width of 0.02 Ry applied according
to the Methfessel-Paxton method. The plane wave energy cut off of 42 Ry allows
reaching an accuracy of 0.2 mRy/atom. As a realization of VCA for the pseudo
potential method, we used the scheme proposed by Ramer and Rappe Ramer and
Rappe (2000). The elastic constants were calculated from the total energies
obtained for the set of unit cell deformations Mehl et al. (1990).

We use a periodic quadrupolar arrangement for a $1/2\langle 111\rangle$ screw
dislocation Bigger et al. (1992) in the cell with basis vectors
$\vec{b}_{1}=9\vec{u}_{1}$, $\vec{b}_{2}=5\vec{u}_{2}$ and
$\vec{b}_{3}=\vec{u}_{3}$, where $\vec{u}_{1}=[\bar{1}10]$,
$\vec{u}_{2}=[\bar{1}\bar{1}2]$ and $\vec{u}_{3}=1/2[111]$. An appropriate
choice of lattice vectors Bigger et al. (1992) reduces the quadrupole cell to
half the size, and we therefore use a cell with basis vectors equal
$\vec{h}_{1}=(\vec{u}_{1}+\vec{u}_{2}+\vec{u}_{3})/2$,
$\vec{h}_{2}=\vec{b}_{2}$ and $\vec{h}_{3}=\vec{b}_{3}$ and 135 atoms. This
unit cell contains only two dislocations with opposite Burgers vectors (Fig.
1). It was demonstrated earlier that this cell size is large enough to
reproduce such characteristics of the dislocations as Peierls stress and
barrier Segall et al. (2001a); Ismail-Beigi and Arias (2000); Segall et al.
(2003); Frederiksena and Jacobsena (2003); Li et al. (2004); Ventelon and
Willaime (2007); Odbadrakh et al. (2011) reasonably well  if the elastic
interaction correction is included. We are interested in how the barriers
change with solute concentration and therefore the above correction is not
included since it’s the same for all concentrations. For the dislocation
calculation, the BZ summation was carried out over a $1\times 2\times 8$ BZ
grid Romaner et al. (2010) and the initial structure was relaxed until the
forces were smaller than 0.0005 Ry/Å.

##  III Results

###  III.1 Electronic structure

The electronic density of states (DOS) calculated using the VCA are presented
in Figure 2 (colored online) for pure tungsten (blue solid line) and two
tungsten alloys, one with 6.25 % Re (red dashed line) and the second one with
6.25 % Zr (green dash-dot line). Zero energy on this plot corresponds to
occupation of electronic bands by six electrons. Thus, for pure tungsten, zero
energy corresponds to the Fermi level. The Fermi energy is placed in the
pseudo gap between bonding and antibonding $d$-states. As can be seen from
Fig. 2, alloying with Re or Zr doesn’t produce sizable changes in the DOS, at
least in VCA approach. The substitution of tungsten with 6.25 % Re leads to an
upwards shift of the Fermi energy, shown by a red vertical line. The area
under DOS between zero energy and the W0.9375Re0.0625 Fermi energy corresponds
to 0.0625 additional electrons, i.e., occupation of each virtual atom by an
additional 0.0625 electrons. Correspondingly, alloying with atoms that have a
lower number of electrons, such as Zr, leads to a downward shift of the Fermi
energy, shown by vertical green line, and each virtual atom is occupied by
0.0625 fewer electrons.

Such a modification of electronic structure in tungsten alloys can be
described using the rigid band approximation Faulkner (1982). This approach
assumes that alloying does not change the DOS but shifts the Fermi level, so
the occupation of each virtual atom corresponds to the average number of
electrons per atom $(e/a)$ in the alloy. Thus, within this approximation, all
alloys with the same $(e/a)$ value have the same properties. The concentration
of different alloying elements which have the same $(e/a)$ value can be
calculated through the simple expression

$x=\frac{(e/a)}{Z_{TM}-Z_{W}}$ (1)

, where $Z_{TM}$ is the number of valence electrons of the alloying atom and
$Z_{W}$ is the number for tungsten. Although this approach has limited
accuracy, it can be very useful for qualitatively estimating the change in
such properties as elastic constants or Peierls barrier. It should be
mentioned that alloying of tungsten with TM from VIIA and VIIIA groups, within
their solubility limits, fills bonding $d$ states with additional electrons
and increases strength of the bonds.  It reflects in the increase of cohesive
energy and reduction of Wigner-Seitz radius for 4$d$ TM with filling of
bonding $d$ states and reduction of cohesive energy and increase of Wigner-
Seitz radius with filling of antibonding states Moruzzi et al. (1977).

The VCA approach could be quite inaccurate for describing the properties of
disordered substitutional metallic alloys, especially for transition metals
Faulkner (1982) and therefore the accuracy of this approach should be analyzed
in each particular case. In order to verify the accuracy of VCA we compared
its results with those obtained using a supercell approach. The supercell
contained $2\times 2\times 2$ cubic unit cells with one tungsten atom
substituted by a solute atom. Thus, the solute concentration corresponds to
6.25 %. In Figure 3, the DOS calculated using the supercell (blue line) and
VCA models (red line) are presented for W15Fe. The DOS calculated using the
VCA approach is very close to the supercell result. Comparisons of elastic
properties calculated using the two approaches will be discussed below.

###  III.2 Elastic constants

The modification of elastic properties in W-Re alloys has been widely
investigated both experimentally Ayres et al. (1974) and theoretically Persson
et al. (1999); Ekman et al. (2000). Here we expand these investigations to a
wide set of tungsten-based alloys. The results of equilibrium lattice
constant, $a$ , bulk modulus, $B$, and elastic constants, $C_{11}$, $C_{12}$,
$C_{44}$ and $C^{\prime}=(C_{11}-C_{12})/2$, calculated by both VCA and
supercell methods, are presented in Table LABEL:elst_cnst_vca_vs_sc. For the
elastic constants, the largest difference between calculated and experimental
values was obtained for $C_{44}$ and is about 13 % for the VCA calculations
and slightly larger for $2\times 2\times 2$ super cell calculations (SC).  The
agreement between VCA and SC results W alloys with Zr, Ta, and Re is very good
and differences for elastic constants do not exceed 2 %. Hoever this
difference increases for alloys with solute atoms with a larger number of
valence electron. For example, for the Fe-family (Fe, Ru and Os) the
difference is about 5 %, while for the Co- (Co, Rh and Ir) and for Ni-families
(Ni, Pd and Pt) the difference increases to 12-15 %. Since we investigate the
reduction of Peierls stress and barrier for alloys in which the number of
valence electrons per atom $(e/a)$ corresponds to 10 % of Re or smaller the
elastic properties of the alloys should be described with an accuracy similar
to Re of 5 % or less. As discussed above, alloying of tungsten with transition
metals with a higher number of $d$ electrons, metals from the VIIA and VIIIA
groups, fills bonding $d$ states. According to the general results obtained
for transition metals Moruzzi et al. (1977); Faulkner (1982), this filling
leads to a decrease of the lattice parameter. The calculated results for
W1-xRex alloys reproduce this general tendency (see Table LABEL:elst_cnst) and
agrees with experiment Ayres et al. (1974). The same tendency was obtained in
tungsten alloyed with other TM except Pd and Pt. For these two elements, the
lattice parameter calculated in SC slightly increases. Another general trend
for all tungsten alloys with solute atoms having a higher number of $d$
electrons is reduction of the $C^{\prime}$ elastic constant. The calculated
results in Table LABEL:elst_cnst are also in good agreement with experimental
values for the W-Re alloy Ayres et al. (1974) and reflect a reduction in the
stability of bcc structure. At concentrations higher than x=0.25% the bcc
W1-xRex transforms to the $\sigma$-phase Massalski et al. (1986). Similar to
previous results Persson et al. (1999), the $C^{\prime}$ modulus decreases
with Re concentration until it becomes negative at 85 % Re. This change in the
sign of $C^{\prime}$ corresponds to the dynamic loss of stability of the bcc
structure. The $C^{\prime}$ elastic constant corresponds to the long-wave
transversal phonon branch in $[\xi\xi 0]$ directions (T $[1\bar{1}0][\xi\xi]$)
and softening of this phonon mode provides a transition path from the bcc to
dhcp structure (see discussion in Persson et al. (1999)).

According to the rigid band model discussed in the previous section, the
parameter which determines the elastic properties of W1-xTMx alloys is the
number of electrons per atom $e/a=(1-x)Z_{W}+xZ_{TM}$ (see explanation of
Eq.1). In Fig. 4 the elastic constant $C^{\prime}$ and elastic anisotropy
$A=C_{44}/C^{\prime}$ are presented as a function of $e/a$. The $C^{\prime}$
and $A$ values are very similar for all TM from group VIIA and VIIIA except
for Fe and Ir. However, even for these two elements the deviation from the
average values is around 7 %. Additional evidence for the importance of the
$e/a$ parameter follows from the correlation between the width of the stable
bcc solid solution region of W1-xTMx and the number of valence electrons for
TM elements from groups VIIA and VIIIA Massalski et al. (1986). Alloys with
solute atoms having a higher number of valence electrons have a narrower solid
solution region.

###  III.3 Dislocation core structure, Peierls stress and Peierls barrier.

In bcc structures $1/2\langle 111\rangle$ screw dislocations have three
nonequivalent core configurations, that are referred to as ”easy”, ”hard” and
”split” Vitek (1974); Xu and Moriarty (1996); Ismail-Beigi and Arias (2000);
Takeuchi (1979) \- stable and two metastable, respectively. These
configurations can be obtained for a given Burgers vector $\vec{b}$ by placing
the core origin in specific positions. Recently, it was demonstrated that in
bcc Fe the lowest energy is the symmetric ”easy” core configuration, then
symmetric ”hard” core and the highest energy is the ”split” core configuration
Itakura et al. (2012). In a contrast to Fe, for the case of W we determined
that the symmetric ”hard” core configuration decays into the ”split” core.
Thus, in the present publication we assume that the dislocation migrates
between the ”easy” and the ”split” core configurations.

Following, Xu and Moriarty Xu and Moriarty (1998) we marked these sites as 1,
4 or 5 for ”easy” cores and site 2 or 3 for ”split” cores in Figure 5. The
calculated structure of the ”easy” core configuration for a $1/2\langle
111\rangle$ dislocation is shown in Fig. 5 for pure W, W0.75Re0.25 and
W0.88Fe0.12, where the concentrations of transition metals were chosen to give
the same $(e/a)$ value. The circles in Fig. 5 represent W atoms looking in the
$\langle 111\rangle$ direction and the dislocation structure is illustrated by
differential displacement (DD) maps Vitek (1974). On a DD map the displacement
of atoms in the [111] direction relative to a neighbor is plotted as an arrow
connecting the atom with its neighbors. The length of an arrow is normalized
so that arrows connecting two atoms corresponds to displacement by $b/3$.
Thus, summing arrows in the circuit around a dislocation gives the Burgers
vector. For tungsten alloys  with VIIIA group TM having the same $(e/a)$
value, the core configurations looks exactly the same as for the Re and Fe
alloys. On the contrary, alloying with Ta or Zr produces the same core
symmetry as pure tungsten. The same result was obtained for W-Ta by Li et al
Li et al. (2012). The core structure of pure W is symmetric, as shown in
Figure 5a, i.e., the dislocation expands equally along the six $\langle
112\rangle$ directions, similar to results obtained earlier Romaner et al.
(2010).

Alloying with Re or any other TM with a higher number of valence electrons
leads to a change in dislocation core structure from symmetric to asymmetric.
The core is spread out in three $\langle 112\rangle$ directions on $\\{110\\}$
planes Xu and Moriarty (1998). There are a six possible $\langle 112\rangle$
orientations and  the core is thus double degenerate. The transition from
symmetric to asymmetric core can change the dislocation slip plane Romaner et
al. (2010). The symmetric core dislocations glide uniformly on $\\{110\\}$
planes; asymmetric ones glide in a zigzag manner Hirth and Lothe (1982) and
the slip plane changes to $\\{112\\}$. Using the core site notation in Figure
5,  the symmetric core glide path $1\rightarrow 2\rightarrow 4$, which is
shown by the gray band in Figure 5a, will be changed to glide path
$1\rightarrow 2\rightarrow 4\rightarrow 3\rightarrow 5$, shown by the gray
band in Figure 5c Gröger (2007) . In the present work, the transition from
symmetric to asymmetric core was obtained in the VCA calculations and this
result was confirmed for the W-Re alloy in a super cell model calculation
reported earlier Romaner et al. (2010); Li et al. (2012).

Addition of Re or any TM from group VIIIA leads to the change of dislocation
slip plane and decrease in the value of Peierls stress $\sigma_{p}$ Romaner et
al. (2010). Following a widely used technique Segall et al. (2001b),  we apply
pure shear strain in the $\langle 111\rangle$ direction which results in a
stress that has $\sigma_{zx}$ as its main component, the influence of
$\sigma_{yx}$ can be neglected since $\sigma_{yx}/\sigma_{zx}=0.04$ for pure
W. The corresponding strain is produced by modification of the basis vector
$\vec{h}_{1}$ along $\vec{u}_{3}$ direction
$\vec{h}_{1}=1/2\vec{u}_{1}+1/2\vec{u}_{2}+(1/2+\varepsilon)\vec{u}_{3}$.
Figure 6 shows the dependence of total energy (see inset) per dislocation per
Burgers vector and shear stress as a function of strain, $\varepsilon$.

For small strains, i. e. in the elastic regime, the energy increases as
$\varepsilon^{2}$. At larger $\varepsilon$ the energy dependence deviates from
square dependence and abruptly drops. This drop in energy is caused by a jump
of the dislocation core to  the next stable ”easy” core neighboring position
($1\rightarrow 4$). The corresponding stress and strain are considered to be
the Peierls stress $\sigma_{p}$ and strain $\varepsilon_{p}$. Earlier it was
demonstrated Romaner et al. (2010) that the cell size of 135 atoms as used in
the present work is enough to reach convergence in $\sigma_{p}$ values.
Alloying with both Re and Fe reduces values of $\varepsilon_{p}$ and
$\sigma_{p}$. For Fe this reduction is even larger than for Re although the
solute concentrations corresponds to the same $e/a$ value in both cases. The
absolute value of $\sigma_{p}$ equal to 1.71 GPa for pure tungsten (Table
LABEL:PS)) is somewhat lower than those obtained earlier Romaner et al.
(2010), which were 2.09 GPa using QE and 2.37 GPa using VASP.  The difference
is attributed to the use of different types of pseudo-potentials, ultra-soft
in our case and norm-conserving (QE) and PAW (VASP) in Romaner et al. Romaner
et al. (2010), and the higher accuracy obtained using our larger plane wave
energy cut-off of 42 Ry compared to 30 Ry (QE calculations) and 16.4 Ry (VASP
calculations) Romaner et al. (2010). Moreover, the norm-conserving
pseudopotential would require an even larger cut-off energy than the ultra-
soft one. The relative reduction of $\sigma_{p}$ is 20 % for Re and 36 % for
Fe. A larger reduction of $\sigma_{p}$ in tungsten-iron alloys correlates with
the larger reduction of $C^{\prime}$ modulus caused by alloying with Fe
(Figure 4 and Table LABEL:elst_cnst). These results demonstrate that even if
the reduction of Peierls stress and $C^{\prime}$ values can be qualitatively
understood as arising from a filling of $d$-states in tungsten by additional
electrons from group VII and VIII solutes, some other element-specific
mechanism exists. It should be mentioned that the solubility limit of Fe in W
is around 2% and in the concentration of Fe in limited to this value the
reduction of Peierls barrier is about one half of that of the alloy with 5%
Fe.  For the W-Ru alloy with a solubility limit of 3 %, $\sigma_{p}$ is
reduced by 6 %. This result is consistent with the experimentally observed
reduction of ductile-brittle transition temperature in W-Ru alloy O’Dell et
al. (2012).

The second method used in present work to estimate Peierls stress was
originally proposed by Nabarro Nabarro (1947), and has been discussed in
publications by Gröger et al. Gröger (2007); Gröger and Vitek (2012).
Following Nabarro, we measure the system energy as a function of the
dislocation location as it moves from position 1 to position 2 in Figure 5.
It is easy to see that both $1\rightarrow 2\rightarrow 4$ ($\\{110\\}$ slip
plane) and $1\rightarrow 2\rightarrow 4\rightarrow 3\rightarrow 5$
($\\{112\\}$ slip plane) paths contain a set of dislocation core jumps between
”easy” and ”split” core configurations with the same $1\rightarrow 2$ barrier.
The energy dependence obtained is associated with the Peierls barrier and the
Peierls stress is given by the maximum gradient of this function

$\sigma_{p}=\frac{1}{b}\left(\frac{dV}{dR_{c}}\right)_{max},$ (2)

where $b$ is the Burgers vector,  $V$ is energy per unit length of a straight
dislocation and $R_{c}$ is the distance along the dislocation core path on the
move from 1 to 2. $R_{c}$ equals 0 for position 1 and $\sqrt{2}a/3$, where $a$
is the lattice parameter, for the position 2. We use the drag method, also
called the reaction coordinate method Xu and Moriarty (1998), to move the
dislocation from 1 (”easy”) to 2 (”split”) core configuration. According to
this method, the $z$-coordinates of atoms in two columns around both cores in
the modeling cell, shown by purple circles in Fig. 5a, are fixed so that
$z=r_{c}z_{0}+(1-r_{c})z_{1}$, where $z_{0}$ and $z_{1}$ correspond to the
$z$-coordinate of the specified atoms in the ”easy” and ”split” core
configurations respectively, and $r_{c}=R_{c}/(\sqrt{2}a/3$) is the so called
reaction coordinate. In the initial implementation Xu and Moriarty (1998), the
authors fix the $z$-coordinate of atom which has the largest displacement
resulting from the movement of the dislocation core. However from our
experience the atomic relaxation convergence increases if we fix the
$z$-coordinate of any two of the three atoms surrounding site 2 in Fig. 5a. In
order to avoid any changes in the elastic interaction between dislocations,
both cores are moved simultaneously. The results obtained for the Peierls
barrier changes are presented in Figure 7 for pure tungsten and alloys of
tungsten with Re, Fe, Ru and Zr. In order to demonstrate that the ”split” core
configuration corresponds to a metastable state, we calculated the potential
barrier using a small reaction coordinate step for pure tungsten around
$r_{c}=1$. As can be seen from the results presented in the upper insert in
Figure 7, the ”split” core configuration corresponds to a local minimum in the
potential energy. When normalized using the values for pure tungsten, the
shape of the Peierls barrier for the alloys is shown to be similar to tungsten
and quite insensitive to solute type and concentration, see lower insert of
Figure 7. In all these alloys, the results yield a single barrier. For alloys
with Os, Ir, Pt, Ru, Rh, Pd, Co and Ni the total energy was calculated for the
”easy” core (reaction coordinate, $r_{c}=0$) and ”split” core ($r_{c}=1$).
Since the shape of Peierls barrier is insensitive to type of solute, the
difference between these energies provides an estimate of the value of the
Peierls barrier as shown in Figure 8. The reduction of the Peierls barrier
observed in these calculations correlates well with the reduction in Peierls
stress estimated in the deformation simulation described above. Alloying with
Re or any group VIII transition metal at a concentration which corresponds to
$(e/a)=6.10$ reduces the barrier by 25 %. Fe reduces the Peierls barrier even
more significantly, whereas Zr increases the barrier. This is consistent with
the effect obtained for Mo alloys with 5$d$ transition metals in which placing
a row of solute atoms with higher number of valence $d$ electrons along the
dislocation core reduce stiffness and a row of atoms with lower number
increase stiffness Trinkle and Woodward (2005).

The Peierls stress calculated using the direct deformation method,
${{\sigma}}_{p}$, and from the Peierls barrier using Eq.2,
${{\tilde{\sigma}}}_{p}$ are presented in Table LABEL:PS. For all alloys the
Peierls stress values obtained from the Peierls barrier are lower than those
from the direct method. The difference is largest, 14 %, for W alloys with 10
% Re, or 2 % Fe. Thus, the ${{\tilde{\sigma}}}_{p}$ results, introduction of 2
% Fe to tungsten reduces the Peierls stress by 12 %. While for the pure W and
W alloy with 5 % of Fe it is below 5 %. Gröger and Vitek Gröger and Vitek
(2012) connected the discrepancy in the results of two methods with the fact
that both drag and Nudged Elastic Band Jonsson et al. (1998) (NEB) methods
garantee that images of the system are distributed uniformly along the minimum
energy path and do not implay that dislocation position is distributed along
this path. In order to overcome this problem, the authors introduced a
modified NEB Gröger and Vitek (2012) method which gives a Peierls stress that
agrees within 8 % of the directly calculated stress. However we believe that
the 14 % accuracy obtained in our calculation is reasonably good considering
that maintains the same trend in ${{\sigma}}_{p}$. It should be metioned that
problems with direct application of drag or NEB methods can be more
substentional. According to modeling result obtained in the deformed cell the
dislocation core structure is slightly modified with increase of deformation
and then jumps directly to next ”easy” core position. The dislocation core is
not observed in metastable, ”split” core position. Thus the agreemnt between
data presented in Table LABEL:PS gives a support to application of drag method
to calculation of Peierls stress.

The reduction of Peierls stress/barrier due to alloying with transition metals
with a higher number of valence electrons and the increase for lower number of
valence electrons is supported by existing experimental data. Hardening was
observed in W-Zr alloys Pod’yachev and Gavrilyuk (1975) and softening in
W-Re/Ir alloys Luo et al. (1991). Very interesting experimental results for Mo
alloys with group VIIIA $3d$ transition metals were presented by Hiraoka et al
Hiraoka et al. (2004). The authors report a decrease of yield strength with
addition of any element from group VIIIA, Fe, Co, Ni or $4d$ Pd. The author
concluded that this effect could be well understood in the terms number of
valence electrons of alloying elements, and that the effect of atomic size
mismatch is secondary and minor. Although the value of yield strength is not
directly related to ductility, the observed correlation of properties with
$e/a$ value supports the validity of our calculation using the VCA approach.

##  IV Conclusions

The influence of alloying tungsten with transition metal solutes on the
elastic properties and properties of $1/2\langle 111\rangle$ screw dislocation
based on electronic structure calculations. In comparison with a standard
supercell method, we demonstrated that the virtual crystal approximation gives
a fairly good description of such alloying, especially in the case of 4$d$ and
5$d$ transition metals. For the case of alloying with transition metals from
group VIIIA, the modification of elastic constants, Peierls stress and barrier
can be understood within the rigid band approximation. This means that alloys
with TM concentrations leading to the same number of electrons per atoms
exhibit a similar reduction of $C^{\prime}$ modulus and elastic anisotropy
$A$. Together with this modification of elastic moduli the value of the
Peierls stress and barrier are reduced by alloying with transition metals from
group VIIIA. In addition to the Peierls barrier reduction, alloying with any
metal from this group changes the ”easy” core dislocation structure from the
symmetric to asymmetric configuration. This similarity allows us to conclude
that the search for an alternative to Re could be reduced to transition metals
from the group VIIIA.

##  V Acknowledgments

Work was supported by the U.S. Department of Energy Office of Fusion Energy
Sciences. This research used resources of the National Energy Research
Scientific Computing Center, which is supported by the Office of Science of
the U.S. Department of Energy.

## References

  * El-Guebaly et al. (2011) L. El-Guebaly, R. Kurtz, M. Rieth, H. Kurishita, A. Robinson, and A. Team, Fusion Science and Technology 60, 185 (2011).
  * Nygren et al. (2011) R. Nygren, R. Raffray, D. Whyte, M. Urickson, M. Baldwin, and L. Snead, Journal of Nuclear Materials 417, 451 (2011).
  * Klopp (1969) W. D. Klopp, J. Metals 21, 23 (1969).
  * Geach and Hughes (1956) J. C. Geach and J. E. Hughes, _Plansee Proc._ (Pergamon Press, Oxford, 1956).
  * Savitskii et al. (1965) E. M. Savitskii, M. A. Tylkina, and K. V. Povarova, _Rhenium-Based Alloys_ (Izd. Nauka, Moskva, 1965).
  * Edington et al. (1966) J. W. Edington, A. Gilbert, and M. J. Klein, NASA CR-81225, 1 (1966).
  * Raffo (1969) P. Raffo, Journal of the Less Common Metals 17, 133 (1969).
  * Luo et al. (1991) A. Luo, D. Jacobson, and K. Shin, International Journal of Refractory Metals and Hard Materials 10, 107 (1991).
  * Gornostyrev et al. (1991) Y. N. Gornostyrev, M. I. Katsnelson, G. V. Peschanskikh, and A. V. Trefilov, physica status solidi (b) 164, 185 (1991).
  * Trefilov et al. (1975) V. I. Trefilov, Y. V. Milman, and S. V. Firsov, _Physics of Refractory Metal Strength_ (Izd. Naukova Dumka, Kiev, 1975).
  * Kurdyumova et al. (1980) G. G. Kurdyumova, Y. V. Milman, and V. I. Trefilov, _Electronic structure and Physico-Chemical Properties of Refractory Compounds and Alloys_ (Izd. Naukova Dumka, Kiev, 1980).
  * Gao et al. (2008) M. C. Gao, mer N. Dogan, P. King, A. D. Rollett, and M. Widom, JOM 60, 61 (2008).
  * Hiraoka et al. (2004) Y. Hiraoka, T. Ogusu, and N. Yoshizawa, Journal of Alloys and Compounds 381, 192 (2004).
  * Klopp (1975) W. D. Klopp, Journal of the Less Common Metals 42, 261 (1975).
  * Stephens and Witzke (1975) J. R. Stephens and W. R. Witzke, Journal of the Less Common Metals 41, 265 (1975).
  * Medvedeva et al. (2005) N. I. Medvedeva, Y. N. Gornostyrev, and A. J. Freeman, Phys. Rev. Lett. 94, 136402 (2005).
  * Medvedeva et al. (2007) N. I. Medvedeva, Y. N. Gornostyrev, and A. J. Freeman, Phys. Rev. B 76, 212104 (2007).
  * Romaner et al. (2010) L. Romaner, C. Ambrosch-Draxl, and R. Pippan, Phys. Rev. Lett. 104, 195503 (2010).
  * Li et al. (2012) H. Li, S. Wurster, C. Motz, L. Romaner, C. Ambrosch-Draxl, and R. Pippan, Acta Materialia 60, 748 (2012).
  * Trinkle and Woodward (2005) D. R. Trinkle and C. Woodward, Science 310, 1665 (2005).
  * Garfinkle (1966) M. Garfinkle, NASA TN D-3190, 1 (1966).
  * Faulkner (1982) J. S. Faulkner, Progress in Materials Science 27, 1 (1982).
  * Persson et al. (1999) K. Persson, M. Ekman, and G. Grimvall, Phys. Rev. B 60, 9999 (1999).
  * Ekman et al. (2000) M. Ekman, K. Persson, and G. Grimvall, J. Nucl. Mater. 278, 273 (2000).
  * Lassner and Schubert (1999) E. Lassner and W.-D. Schubert, _Tungsten - Properties, Chemistry, Technology of the Element, Alloys, and Chemical Compounds_ (Springer-Verlag, 1999).
  * Giannozzi et al. (2009) P. Giannozzi, S. Baroni, and N. B. …, J of Physics: Condensed Matter 21, 395502 (2009).
  * Rappe et al. (1990) A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos, Phys. Rev. B 41, 1227 (1990).
  * Perdew and Wang (1992) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
  * Ramer and Rappe (2000) N. J. Ramer and A. M. Rappe, Phys. Rev. B 62, R743 (2000).
  * Mehl et al. (1990) M. J. Mehl, J. E. Osburn, D. A. Papaconstantopoulos, and B. M. Klein, Phys. Rev. B 41, 10311 (1990).
  * Bigger et al. (1992) J. R. K. Bigger, D. A. McInnes, A. P. Sutton, M. C. Payne, I. Stich, R. D. King-Smith, D. M. Bird, and L. J. Clarke, Phys. Rev. Lett. 69, 2224 (1992).
  * Segall et al. (2001a) D. Segall, T. Arias, A. Strachan, and W. Goddard, Journal of Computer-Aided Materials Design 8, 161 (2001a).
  * Ismail-Beigi and Arias (2000) S. Ismail-Beigi and T. A. Arias, Phys. Rev. Lett. 84, 1499 (2000).
  * Segall et al. (2003) D. E. Segall, A. Strachan, W. A. Goddard, S. Ismail-Beigi, and T. A. Arias, Phys. Rev. B 68, 014104 (2003).
  * Frederiksena and Jacobsena (2003) S. L. Frederiksena and K. W. Jacobsena, Phil. Mag. 83, 365 (2003).
  * Li et al. (2004) J. Li, C.-Z. Wang, J.-P. Chang, W. Cai, V. V. Bulatov, K.-M. Ho, and S. Yip, Phys. Rev. B 70, 104113 (2004).
  * Ventelon and Willaime (2007) L. Ventelon and F. Willaime, Journal of Computer-Aided Materials Design 14, 85 (2007).
  * Odbadrakh et al. (2011) K. Odbadrakh, A. Rusanu, G. M. Stocks, G. D. Samolyuk, M. Eisenbach, Y. Wang, and D. M. Nicholson, Journal of Applied Physics 109, 07E159 (2011).
  * Moruzzi et al. (1977) V. L. Moruzzi, A. R. Williams, and J. F. Janak, Phys. Rev. B 15, 2854 (1977).
  * Ayres et al. (1974) R. A. Ayres, G. W. Shanette, and D. F. Stein, J. Appl. Phys 46, 1528 (1974).
  * Massalski et al. (1986) T. B. Massalski, J. L. Murray, L. H. Bennett, and H. Baker, _Binary alloy phase diagrams_ (American Society for Metals, Metals Park, Ohio, 1986).
  * Vitek (1974) V. Vitek, Cryst. Lattice Defects 5, 1 (1974).
  * Xu and Moriarty (1996) W. Xu and J. A. Moriarty, Phys. Rev. B 54, 6941 (1996).
  * Takeuchi (1979) S. Takeuchi, Philosophical Magazine A 39, 661 (1979).
  * Itakura et al. (2012) M. Itakura, H. Kaburaki, and M. Yamaguchi, Acta Materialia 60, 3698 (2012).
  * Xu and Moriarty (1998) W. Xu and J. A. Moriarty, Computational Materials Science 9, 348 (1998).
  * Hirth and Lothe (1982) J. Hirth and J. Lothe, _Theory of dislocations_ (Wiley-Interscience, New York, 1982).
  * Gröger (2007) R. Gröger, Ph.D. thesis, University of Pennsylvania, Pennsylvania (2007), eprint arXiv:cond-mat/0707.3577v1.
  * Segall et al. (2001b) D. E. Segall, T. Arias, A. Strachan, and W. Goddard III, Journal of Computer-Aided Materials Design 8, 161 (2001b).
  * O’Dell et al. (2012) S. O’Dell, T. McKechnie, A. Shchetkovskiy, and A. Smirnov, in _Workshop on Challenges to Developing W-Based Materials fo Fusion Application._ (UC Santa Barbara, Santa Barbara., 2012).
  * Nabarro (1947) F. R. N. Nabarro, Proceedings of the Physical Society 59, 256 (1947).
  * Gröger and Vitek (2012) R. Gröger and V. Vitek (2012), eprint arXiv:cond-mat/111.5709v1.
  * Jonsson et al. (1998) H. Jonsson, G. Mills, and K. W. Jacobsen, in _Classical and Quantum Dynamics in Condensed Phase Simulations_ , edited by E. B. J. Berne, G. Ciccotti, and D. F. Coker (World Scientific, 1998), p. 385\.
  * Pod’yachev and Gavrilyuk (1975) V. N. Pod’yachev and M. I. Gavrilyuk, Metal Science and Heat Treatment 17, 299 (1975).

###### List of Figures

  1. 1 The [111] projection of dislocation dipole unit cell, where atoms are shown by circles and arrows corresponds to differential displacements.
  2. 2 DOS for pure W, blue solid line, W0.9375Re0.0625, red dashed line, and W0.9375Zr0.0625, green dash-dotted line, calculated using VCA.
  3. 3 W0.9375Fe0.0625 DOS per atom calculated using a $2\times 2\times 2$ supercell model, blue dashed line, and VCA, red solid line. The Fermi energy corresponds to zero.
  4. 4 $C^{\prime}$ and $A$ as a function of number of electrons per atom $e/a$ calculated for W1-xTMx using VCA.
  5. 5  The DD map of ”easy” core dislocation in W1-xTMx calculated using the VCA approach. Some of arrows are highlighted in green to illustrate a change of core symmetry. The sites 1, 4 and 5 correspond to the stable ”easy” core dislocation center, while sites 2 and 3 indicate the metastable ”split” core center. The $z$-coordinate of atoms in columns shown by filled circles are fixed during calculation of Peierls barrier within the reaction coordinate method. The symmetric and asymmetric core glide paths are shown by gray bands.
  6. 6 Shear stress and total energy, $\Delta E$, per dislocation per Burgers vector $\vec{b}$ (shown in the inset) as a function of strain, $\varepsilon$, for pure $W$ (shown by squares), W0.90Re0.10 (shown by  circles), W0.95Fe0.05 (shown by up pointing triangles) and W0.98Fe0.02 (shown by down pointing triangles).
  7. 7 Calculated Peierls barrier, where reaction coordinate, $r_{c}$, equal 0 for ”easy” core configuration and 1 for ”split” core. The curves are normalized in the center insert. The detailed structure of the calculated Peierls barrier in pure W around maximum is shown in the insert at upper right.
  8. 8 The height of Peierls barrier for the set of W alloys.

###### List of Tables

  1. 1 Comparison of lattice parameter ($a$) in Å, bulk modulus ($B$) and elastic constants ($C_{ij}$) in GPa calculated using the VCA and supercell (SC) approaches for W1-xTMx. The concentration $x=0.0625$ corresponds to W15TM.
  2. 2 Experimental and calculated(based on VCA) values of lattice parameter ($a$) in Å, bulk modulus ($B$) and elastic constants ($C_{ij}$) in GPa for W1-xTMx.  The concentration for all TM except Re is chosen to keep number of electrons per atom $(e/a)=6.10$.
  3. 3 In order to compare the reduction of Peierls stress calculated using the direct deformation method, ${\sigma}_{p}$ in GPa, and from the Peierls barrier using Eq.2, ${\tilde{\sigma}}_{p}$ in W, W0.90Re0.10, W0.95Fe0.05, W0.98Fe0.02 and W0.97Ru0.03.

Figure 1: The [111] projection of dislocation dipole unit cell, where atoms
are shown by circles and arrows corresponds to differential displacements.

Figure 2: DOS for pure W, blue solid line, W0.9375Re0.0625, red dashed line,
and W0.9375Zr0.0625, green dash-dotted line, calculated using VCA.

Figure 3: W0.9375Fe0.0625 DOS per atom calculated using a $2\times 2\times 2$
supercell model, blue dashed line, and VCA, red solid line. The Fermi energy
corresponds to zero.

Figure 4: $C^{\prime}$ and $A$ as a function of number of electrons per atom
$e/a$ calculated for W1-xTMx using VCA.

Figure 5:  The DD map of ”easy” core dislocation in W1-xTMx calculated using
the VCA approach. Some of arrows are highlighted in green to illustrate a
change of core symmetry. The sites 1, 4 and 5 correspond to the stable ”easy”
core dislocation center, while sites 2 and 3 indicate the metastable ”split”
core center. The $z$-coordinate of atoms in columns shown by filled circles
are fixed during calculation of Peierls barrier within the reaction coordinate
method. The symmetric and asymmetric core glide paths are shown by gray bands.

Figure 6: Shear stress and total energy, $\Delta E$, per dislocation per
Burgers vector $\vec{b}$ (shown in the inset) as a function of strain,
$\varepsilon$, for pure $W$ (shown by squares), W0.90Re0.10 (shown by
circles), W0.95Fe0.05 (shown by up pointing triangles) and W0.98Fe0.02 (shown
by down pointing triangles).

Figure 7: Calculated Peierls barrier, where reaction coordinate, $r_{c}$,
equal 0 for ”easy” core configuration and 1 for ”split” core. The curves are
normalized in the center insert. The detailed structure of the calculated
Peierls barrier in pure W around maximum is shown in the insert at upper
right.

Figure 8: The height of Peierls barrier for the set of W alloys.

Table 1: Comparison of lattice parameter ($a$) in Å, bulk modulus ($B$) and elastic constants ($C_{ij}$) in GPa calculated using the VCA and supercell (SC) approaches for W1-xTMx. The concentration $x=0.0625$ corresponds to W15TM. TM | Approach | $a$ | $B$ | $C^{\prime}$ | $C_{44}$ | $C_{11}$ | $C_{12}$
---|---|---|---|---|---|---|---
pure W | VCA | 3.1903 | 304 | 160 | 141 | 518 | 197
| SC | 3.1896 | 308 | 157 | 139 | 517 | 203
Zr | VCA | 3.2095 | 287 | 153 | 126 | 491 | 185
| SC | 3.2079 | 286 | 150 | 127 | 489 | 186
Ta | VCA | 3.1946 | 298 | 161 | 136 | 513 | 191
| SC | 3.1963 | 308 | 157 | 135 | 517 | 203
Fe | VCA | 3.1802 | 284 | 132 | 128 | 459 | 196
| SC | 3.1716 | 300 | 140 | 132 | 487 | 207
Co | VCA | 3.1788 | 284 | 125 | 115 | 450 | 200
| SC | 3.1725 | 297 | 129 | 131 | 469 | 211
Ni | VCA | 3.1830 | 279 | 116 | 126 | 434 | 202
| SC | 3.1743 | 295 | 126 | 131 | 463 | 211
Ru | VCA | 3.1772 | 308 | 142 | 142 | 498 | 213
| SC | 3.1842 | 301 | 135 | 137 | 481 | 211
Rh | VCA | 3.1710 | 309 | 135 | 141 | 490 | 219
| SC | 3.1862 | 298 | 124 | 133 | 462 | 215
Pd | VCA | 3.1639 | 310 | 129 | 139 | 482 | 224
| SC | 3.1901 | 293 | 116 | 131 | 448 | 216
Re | VCA | 3.1857 | 308 | 151 | 144 | 510 | 207
| SC | 3.1859 | 308 | 150 | 143 | 508 | 208
Os | VCA | 3.1809 | 312 | 143 | 146 | 503 | 216
| SC | 3.1848 | 307 | 136 | 144 | 489 | 216
Ir | VCA | 3.1736 | 316 | 147 | 148 | 513 | 218
| SC | 3.1865 | 304 | 125 | 140 | 470 | 221
Pt | VCA | 3.1687 | 317 | 132 | 148 | 494 | 229
| SC | 3.1917 | 299 | 113 | 135 | 449 | 223

Table 2: Experimental and calculated(based on VCA) values of lattice parameter ($a$) in Å, bulk modulus ($B$) and elastic constants ($C_{ij}$) in GPa for W1-xTMx.  The concentration for all TM except Re is chosen to keep number of electrons per atom $(e/a)=6.10$. TM | $x$ | $a$ | $B$ | $C^{\prime}$ | $C_{44}$ | $C_{11}$ | $C_{12}$
---|---|---|---|---|---|---|---
experiment Ayres et al. (1974)
W | 0.0000 | 3.1659 | 314 | 164 | 163 | 533 | 205
calculated, virtual crystal approximation
Re | 0.0000 | 3.1903 | 304 | 160 | 141 | 518 | 197
| 0.0300 | 3.1880 | 306 | 156 | 142 | 514 | 202
| 0.0500 | 3.1866 | 307 | 153 | 143 | 512 | 205
| 0.1000 | 3.1831 | 310 | 146 | 145 | 504 | 212
| 0.3000 | 3.1689 | 322 | 130 | 156 | 495 | 235
| 1.0000 | 3.1317 | 350 | -30 | 163 | 309 | 370
Zr | 0.0250 | 3.1972 | 297 | 161 | 136 | 511 | 190
Ta | 0.0500 | 3.1936 | 299 | 162 | 137 | 515 | 191
Fe | 0.0250 | 3.1849 | 296 | 148 | 136 | 494 | 198
Ru | 0.0250 | 3.1846 | 306 | 153 | 141 | 510 | 204
Rh | 0.0167 | 3.1846 | 306 | 153 | 141 | 509 | 204
Os | 0.0250 | 3.1862 | 307 | 153 | 143 | 511 | 205
Ir | 0.0167 | 3.1849 | 308 | 154 | 143 | 512 | 205

Table 3: In order to compare the reduction of Peierls stress calculated using the direct deformation method, ${\sigma}_{p}$ in GPa, and from the Peierls barrier using Eq.2, ${\tilde{\sigma}}_{p}$ in W, W0.90Re0.10, W0.95Fe0.05, W0.98Fe0.02 and W0.97Ru0.03.  | $\sigma_{p}$ | ${\tilde{\sigma}}_{p}$
---|---|---
W | 1.71 | 1.64
W0.95Zr0.05 | 2.18 | 1.82
W0.90Re0.10 | 1.37 | 1.15
W0.95Fe0.05 | 1.09 | 1.05
W0.98Fe0.02 | 1.65 | 1.40
W0.97Ru0.03 | 1.60 | 1.35


# Evolutionary tracks of massive stars with different rotation and metallicity
in neutrino H-R diagram

Hao Wang1,2, Chunhua Zhu1, Helei Liu1, Sufen Guo1, Guoliang Lü2,1
1School of Physical Science and Technology, Xinjiang University, Urumqi,
830017, China
2Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 150 Science
1-Street, Urumqi, 830011, China
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Neutrino losses play a crucial role in the evolution of massive stars. We
study the neutrino luminosity of stars ranging from 20 to 90 $\rm M_{\odot}$
from Zero Age Main Sequence (ZAMS) to Fe Core Collapse (FeCC) with different
rotation and metallicity in a neutrino Hertzsprung-Russell diagram. In our
simulations, we consider $\rm\omega/\omega_{crit}=0$ and 0.7 to represent non-
rotation and high rotation, respectively, and set the metallicities to 0.014,
0.001, and 0.0001. During hydrogen burning stages, neutrino luminosity
primarily originates from CNO cycle, and increases with higher stellar mass
while decreasing with increasing metallicity. For the high metallicity models
($Z=0.014$) during the helium burning stage, the reduction of the hydrogen
envelope caused by a larger mass loss rate leads to a gradual decrease in
neutrino luminosity. The rapid rotation results in extra mixing inside massive
stars, which increases the neutrino luminosity during main sequence (MS),
while decreases the neutrino luminosity during helium burning phase.
Simultaneously, the rapid rotation also increases CO core mass, which enhances
the neutrino luminosity during C and O burning phase. We also investigate the
effect of neutrino magnetic moment (NMM) on the massive stars. We find that
the energy loss caused by the NMM does not have effects on the evolutionary
destiny of massive stars, and it does not significant change the compactness
at the time of Fe core collapse.

###### keywords:

neutrinos – stars: massive – stars: rotation – (stars:) Hertzsprung-Russell
and colour-magnitude diagrams – stars: evolution

††pubyear: 2023††pagerange: Evolutionary tracks of massive stars with
different rotation and metallicity in neutrino H-R diagram–Evolutionary tracks
of massive stars with different rotation and metallicity in neutrino H-R
diagram

##  1 Introduction

Massive stars, defined as those with a mass greater than 8 M⊙, undergo a
series of nuclear fusion reactions. These include hydrogen burning via the CNO
cycle (Bahcall, 1989), helium burning via the triple-alpha process, as well as
carbon, oxygen, and silicon burning via the alpha process. These reactions
lead to the production of heavier elements and eventually to iron core
collapse, which can result in a supernova explosion (Woosley et al., 2002).
During hydrogen and helium burning, energy is radiated away via photons
(Kippenhahn et al., 2013), while during carbon, oxygen, and silicon burning,
neutrino cooling becomes the dominant energy loss mechanism due to the rapid
increase in internal temperature and density (Raffelt, 2012). Gaining a
profound understanding of the critical role of neutrino losses in the
evolutionary path of massive stars is indispensable for the progression of our
comprehension of astrophysics (Sukhbold, 2016).

Neutrino losses are sensitive to temperature and therefore play an important
role in the late evolution of massive stars (Kutschera et al., 2009a; Kato et
al., 2020; Fuller & Haxton, 2022). The neutrino delayed heating mechanism is
also an important factor leading to supernova explosions (Couch, 2013; Mirizzi
et al., 2016; Janka, 2017). Due to their extremely small scattering cross
section, neutrinos can penetrate freely from the interior of the star and
provide valuable insights into the temperature, density, and fusion processes
in the core (Farmer et al., 2016; Raj et al., 2020). Pre-supernova neutrinos
can even serve as an early warning of an impending supernova explosion
(Nagakura & Vartanyan, 2022). Recent studies have also shown that neutrino
signals can be used to estimate the radius and mass of neutron stars (Nagakura
& Vartanyan, 2022).

Neutrino detection is rapidly developing as a new astronomical messenger, with
advanced neutrino detectors such as Super Kamiokande (SK) (Ikeda et al.,
2007), Borexino (Bellini et al., 2014), KamLAND (Eguchi et al., 2003), Irvine-
Michigan-Brookhaven(IMB) (Bionta et al., 1987) designed to detect the burst of
neutrinos emitted during supernova explosions. Borexino (Bellini et al.,
2014), SNO+ (Albanese et al., 2021), and SK-GD (Marti Magro & Super-Kamiokande
Collaboration, 2022) play an important role in solar neutrino observations.
While many new generation neutrino detectors, such as JNE (Beacom et al.,
2017), DUNE (Abi et al., 2020), JUNO (Abusleme et al., 2021), aim to explore
new neutrino events.

Usually, neutrino astronomy mainly focuses on the solar neutrino, supernova
neutrino, even or blazar neutrino. However, it is well known that stars
radiate their energy not only by releasing photons from the stellar surface,
but also by radiating neutrinos from the stellar interior. Very recently, Shi
et al. (2020) and Farag et al. (2020) calculated the neutrino luminosity for
stars during different evolutional stages. The former calculated the neutrino
luminosity of stars with initial masses ranging from 1 to 9 $\rm M_{\odot}$
during the main sequence (MS) to asymptotic giant branch (AGB) stages. The
latter gave the evolutional tracks in a neutrino Hertzsprung-Russell (H-R)
diagram for stars with initial mass lower than 40 M⊙.

In fact, the stellar structure and evolution depends on stellar rotation and
metallicity. Rotation can have a significant impact on the evolutionary tracks
of massive stars, by lowering the effective gravity $g_{\mathrm{eff}}$ and
enhancing the luminosity, increases the mass loss rates (Maeder & Meynet,
2000). More importantly, rotational mixing can change the size of stellar
core, which may influence neutrino luminosity. At lower metallicity, lower
mass and angular momentum loss can make rotational mixing more efficient so
that massive star with high rotation ($\omega/\omega_{\mathrm{crit}}=0.7$) can
experience quasi-chemically homogeneous evolution (Brott et al., 2011; Georgy
et al., 2013), it can cause a larger helium core mass and CO core mass
(Limongi & Chieffi, 2018). Therefore it is necessary to include rotation and
metallicity when taking about neutrino luminosity.

Meanwhile, many studies have shown that neutrino could have magnetic moment
because of its non-zero mass (Marciano & Sanda, 1977; Lee & Shrock, 1977;
Fujikawa & Shrock, 1980; Bell et al., 2005). When a non-zero neutrino magnetic
moment (NMM) is present, the neutrino energy loss can increase due to the
additional electromagnetic contribution to the neutrino emissivity (Mori et
al., 2020). This additional energy loss, induced by physics transcending the
Standard Model, can alter the rate of neutrino production, potentially
impacting the star’s structure (Heger et al., 2009). Based on GEMMA experiment
(Beda et al., 2013), the current best experimental limit $\mu_{\nu}<2.9\times
10^{-11}\mu_{B}$, where $\mu_{B}=e/2m_{\mathrm{e}}$ is the Bohr magneton,
$\mathrm{e}=1.602\times 10^{-19}\ \mathrm{C}$ is the electron charge and
$m_{\mathrm{e}}=9.109\times 10^{31}\ \mathrm{Kg}$ is the electron mass.

In this paper, we plan to investigate the effects of rotation and metallicity
on the neutrino luminosity of massive stars, along with the influence of the
NMM on stellar structure and evolution. We aim to explore the neutrino
luminosity of stars with masses ranging from 20 to 90 $\rm M_{\odot}$, from
the MS to Fe core collapse (FeCC) stages. Section 2 gives the details of
software instruments, input physics, reaction networks and model selection. In
section 3 we discuss the variation of our results in the luminosity of
neutrinos with different solar masses, different metallicity, rotations, and
the effect of NMM on the evolution of massive stars. The conclusions are given
in section 4.

##  2 MODEL

We employ the open-source stellar evolution code (MESA, version 10398, Paxton
et al. 2011; Paxton et al. 2013, 2015, 2018; Paxton 2019) to construct 20 to
90 $\rm M_{\odot}$ stellar models. Each stellar model is evolved from the Zero
Age Main Sequence (ZAMS) until FeCC. Following Farag et al. (2020), we assume
FeCC occurring when any location inside the stellar model reaches an infall
velocity of 300 $\rm km\cdot s^{-1}$. The neutrino normalization
$L_{\nu,\odot}=0.02398\cdot L_{\gamma,\odot}=9.1795\times 10^{31}\mathrm{erg\
s}^{-1}$ is adopted in Farag et al. (2020).

###  2.1 mass loss

The mass-loss rates employed in the present paper are consistent with those
used by Brott et al. (2011). For stars hotter than approximately 25 $\rm kK$
with a surface hydrogen mass fraction $X_{\rm{S}}$ > 0.7, we utilize the wind
recipe proposed by Vink et al. (2001). For stars with low surface hydrogen
mass fraction ($X_{\rm{S}}$ <0.4), we apply the mass-loss recipe for Wolf-
Rayet (WR) stars proposed by Hamann et al. (1995), scaled down by a factor of
ten. The mass-loss rate is obtained through a linear interpolation between the
mass-loss prescriptions provided by Vink et al. (2001) and Hamann et al.
(1995) for hot stars with $X_{\mathrm{S}}$ between 0.4 and 0.7.

###  2.2 convection

We use MLT++ to treat convection, this prescription can calculate models of
massive stars up to core collapse (Paxton et al., 2013). Convection boundaries
are calculated by the Ledoux criterion but not Schwarzschild criterion,
mixing-length parameter of 1.5 is adopted (Böhm-Vitense, 1958; Langer, 1991;
Lü et al., 2017; Cui et al., 2018). Semi-convection is a time-dependent
diffusion process, the diffusion coefficient is given by Langer et al. (1983)
as follows:

$\displaystyle
D_{\mathrm{sc}}=\alpha_{\mathrm{sc}}\left(\frac{\mathrm{K}}{6C_{\mathrm{P}}\rho}\right)\frac{\nabla_{\mathrm{T}}-\nabla_{\mathrm{ad}}}{\nabla_{\mathrm{L}}-\nabla_{\mathrm{T}}}$
(1)

where $\nabla_{\mathrm{T}}$ is temperature gradient, $\nabla_{\mathrm{ad}}$ is
adiabatic temperature gradient, K is radiative conductivity, $C_{\mathrm{P}}$
is specific heat at constant pressure, $\alpha_{\mathrm{sc}}$ is a
dimensionless constant. In our work, we set $\rm\alpha_{sc}=0.01$ (Wellstein
et al., 2001). The overshoot area is calculated by the step-function
overshooting which depends on the overshooting parameter $\alpha_{\rm ov}$.
Following Brott et al. (2011), $\alpha_{\rm ov}=0.335$ in our simulations.

###  2.3 neutrino emission

The production of neutrinos inside stars can be roughly divided into two
processes. One is weak nuclear processes, including $\beta^{\pm}$ decay and
$e^{\pm}$ capture (Kutschera et al., 2009b):

$\displaystyle\rm nuclei\ EC\ :\ (Z,A)+e^{-}\longrightarrow(Z-1,A)+\nu_{e}$
(2) $\displaystyle\rm freep\ EC\ :\ p+e^{-}\longrightarrow n+\nu_{e}$ (3)
$\displaystyle\rm PC\ :\ (Z,A)+e^{+}\longrightarrow(Z+1,A)+\bar{\nu}_{e}$ (4)
$\displaystyle\rm\beta^{+}\ :\ (Z,A)\longrightarrow(Z-1,A)+e^{+}+\nu_{e}$ (5)
$\displaystyle\rm\beta^{-}\ :\
(Z,A)\longrightarrow(Z+1,A)+e^{-}+\bar{\nu}_{e}$ (6)

such as $\rm{}^{13}N(,e^{+}\nu_{e})^{13}C$,
$\rm{}^{15}O(,e^{+}\nu_{e})^{15}N$, $\rm{}^{17}F(,e^{+}\nu_{e})^{17}O$ in CNO
cycle, $\rm{}^{18}F(,e^{+}\nu_{e})^{18}O$ during helium burning,
$\rm{}^{25}Al(,e^{+}\nu_{e})^{25}Mg$, $\rm{}^{26}Al(,e^{+}\nu_{e})^{26}Mg$
during C burning, especially isotopes with $\rm A\approx 50-60$ are major
contributor for neutrino luminosity during $\rm\beta$ processes, such as
$\rm{}^{51,52,53}V$, $\rm{}^{53,55,63}Cr$, $\rm{}^{54,55,56}Mn$,
$\rm{}^{53,54,55}Fe$, $\rm{}^{55,56,57}Co$ (Patton et al., 2017). For the
massive star model, the large nuclear reaction network can better track the
neutrinos produced by each nuclear reaction process, so we select mesa_204.net
as the nuclear reaction network.

The other is thermal processes, thermal neutrino energy losses are from (Itoh
et al., 1996), including:

pair annihilation$\rm\left(T>10^{9}\ K\right)$

$\displaystyle e^{+}+e^{-}\longrightarrow\nu+\bar{\nu}$ (7)

photon-neutrino$\rm\left(T<4\times 10^{8}\ K,\rho<10^{5}\ g\ cm^{-3}\right)$

$\displaystyle e^{-}+\gamma\longrightarrow e^{-}+\nu+\bar{\nu}$ (8)

plasma decay$\rm\left(10^{7}<T<10^{8}\ K,10^{4}<\rho<10^{7}\ g\
cm^{-3}\right)$

$\displaystyle\rm\gamma^{*}\longrightarrow\nu+\bar{\nu}$ (9)

bremsstrahlung$\rm\left(10^{8}<\rho<10^{10}\ g\ cm^{-3}\right)$

$\displaystyle e^{-}+(Z,A)\longrightarrow e^{-}+(Z,A)+\nu+\bar{\nu}$ (10)

The MESA nuclear reaction rates are a combination of rates from NACRE (Angulo
et al., 1999) and JINA REACLIB (Cyburt et al., 2010). Treatment of screening
corrections are from (Chugunov et al., 2007), which enhance nuclear reactions
rates in dense plasmas. All the weak reaction rates are based on (Langanke &
Martínez-Pinedo, 2000; Oda et al., 1994; Fuller et al., 1985).

###  2.4 Rotation and metallicity

Rotation and metallicity have great effects on the massive stars. Rapid
rotation may trigger some instabilities in the stellar internal, which can
result in chemical homogenous evolution (Maeder, 1980; Martins et al., 2009,
2013; Hainich et al., 2015; Schootemeijer & Langer, 2018). Following Heger et
al. (2000), we take into account several instabilities induced by rotation
that result in mixing: Eddington-Sweet circulation instability, dynamical and
secular shear instability, and the Goldreich-Schubert-Fricke instability.
Additionally, the transport of angular momentum facilitated by magnetic fields
via the Spruit-Tayler dynamo mechanism is also taken into consideration (Brott
et al., 2011). According to Heger et al. (2000), the efficiency of rotational
mixing, $f_{\mathrm{c}}$, is considered as 1/30. The inhibition parameter
$f_{\mathrm{\mu}}$ of chemical gradient on the efficiency of rotational mixing
process is set to 0.1 (Yoon et al., 2006). In order to discuss the effects of
rotation, we take $\omega/\omega_{\rm crit}$=0 and 0.7 for different
simulations, where $\omega/\omega_{\rm crit}$ is the ratio of rotational
angular velocity to critical rotational angular velocity. Simultaneously,
metallicity can affect the efficiency of chemical homogenous evolution.
Therefore, we select $Z=0.014$ for the high metallicity, 0.001 and 0.0001 for
metal-poor stars.

In addition, the mass-loss rate of rotating massive stars would be enhanced by

$\displaystyle\dot{M}=\left(\frac{1}{1-\omega/\omega_{\text{crit
}}}\right)^{\beta}\dot{M}_{v_{\mathrm{rot}}=0}$ (11)

where $\rm\beta=0.43$ (Langer, 1998).

###  2.5 neutrino magnetic moment

Following Mori et al. (2020), the additional energy loss induced by NMM has
two ways: plasmon decay and neutrino pair production. The energy loss rate due
to the plasmon decay is given by

$\displaystyle\epsilon_{\text{plas
}}^{\mu}=0.318\left(\frac{\omega_{\mathrm{pl}}}{10\mathrm{keV}}\right)^{-2}\left(\frac{\mu_{v}}{10^{-12}\mu_{\mathrm{B}}}\right)^{2}\epsilon_{\mathrm{plas}}$
(12)

where $\rm\epsilon_{\text{plas }}$ is the standard plasmon decay rate (Itoh et
al., 1996) and $\rm\omega_{\text{pl }}$ is the plasma frequency (Raffelt,
1996)

$\displaystyle\omega_{\mathrm{pl}}=28.7\mathrm{eV}\frac{\left(Y_{\mathrm{e}}\rho\right)^{\frac{1}{2}}}{\left(1+\left(1.019\times
10^{-6}Y_{\mathrm{e}}\rho\right)^{\frac{2}{3}}\right)^{\frac{1}{4}}}$ (13)

where $Y_{\mathrm{e}}$ represents the electron mole fraction, $\rm\rho$
denotes the density in units of $\rm g\ cm^{-3}$. The energy loss rate due to
the pair production is given by (Heger et al., 2009)

$\displaystyle\epsilon_{\text{pair }}^{\mu}=1.6\times
10^{11}\mathrm{erg}^{-1}\mathrm{~{}s}^{-1}\left(\frac{\mu_{v}}{10^{-10}\mu_{\mathrm{B}}}\right)^{2}\frac{e^{-\frac{118.5}{T_{8}}}}{\rho_{4}}$
(14)

where $\rm T_{8}$ is the temperature in units of $\rm 10^{8}\ K$ and
$\rm\rho_{4}=\rho/(10^{4}g\ cm^{-3})$. We added these two additional energy
loss calculation formulas to $\rm run\\_star\\_extras$.

##  3 Result

Using MESA, we simulate the structure and evolution of 8 massive star models
including different initial masses from 20 to 90 M⊙ with a mass interval of 10
M⊙, different initial rotational velocities ($\omega/\omega_{\rm crit}=0$ and
0.7), and different metallicities (Z=0.014, 0.001 and 0.0001). Following Farag
et al. (2020), we employ the neutrino H-R diagram to present the neutrino
luminosities throughout the entire life cycles of these massive stars, and
discuss how stellar rotational velocity and metallicity affect stellar
neutrino luminosity.

###  3.1 H-R diagram and neutrino H-R diagram of massive star

Figure 1 shows the 20 to 90 M⊙ stellar evolutionary tracks in both the H-R
diagram and the neutrino H-R diagram. It clearly illustrates that the neutrino
luminosity during the MS phase increases correspondingly as the stellar mass
increases. In massive stars, the primary source of neutrinos is the Carbon-
Nitrogen-Oxygen (CNO) cycle, which is highly temperature-sensitive
($\varepsilon_{\rm CNO}\propto T^{17}$) (Bahcall, 1989). As a result, nuclear
reactions are more intense in more massive stars, which leads to higher
neutrino luminosity.

As shown in Figure 1, the neutrino luminosity of the star decreases steadily
during the helium burning phase. During this phase, only one weak process,
namely $\rm{}^{18}F(,e^{+}\nu_{e})^{18}O$, generates neutrinos in the
3$\alpha$ reaction (Odrzywolek & Heger, 2010). Consequently, the neutrino
luminosity during the helium burning phase arises primarily from the burning
of the hydrogen shell. At this stage, the high mass loss rate leads to a
reduction in the hydrogen envelope, a decrease in the hydrogen shell
luminosity, and a corresponding decrease in the neutrino luminosity. For more
massive stars, this phenomenon is even more pronounced due to greater mass
loss rate.

Figure 1: The evolution of 20 $\rm M_{\odot}$ to 90 $\rm M_{\odot}$ star with
non-rotation and high metallicity ($Z=0.014$) in the H-R diagram (top panels)
and neutrino H-R diagram (bottom panels), the left is 20 $\rm M_{\odot}$ to 50
$\rm M_{\odot}$ star, right is 60 $\rm M_{\odot}$ to 90 $\rm M_{\odot}$ star.
The positions of the neutrino luminosity at the ignition of He, C, O, and Si
are given.

Then, the neutrino luminosity experiences a rapid increase following carbon
burning. Figure 1 also marks the neutrino luminosity of 20 $\rm M_{\odot}$ and
60 $\rm M_{\odot}$ stars when helium depletion and durning carbon, oxygen, and
silicon burning. When the Carbon is ignited, the internal temperature of the
massive star exceeds $8\times 10^{8}\rm\ K$, leading to the domination of
thermal process, particularly electron-positron pair annihilation, in neutrino
production. At this stage, massive stars are often called neutrino-cooled
stars (Odrzywolek et al., 2004). This is due to the trapping of photons as the
density inside the star increases, while neutrinos can freely escape with
their extremely small cross sections. In addition, neutrinos carry away a
significant amount of energy during C, Ne, O, and Si burning, which
drastically reduces the lifetime of nuclear burning. As the temperature inside
the star increases rapidly, the thermal neutrino energy loss becomes more
significant. Therefore, the sensitivity of neutrino energy loss to temperature
allows massive stars to evolve more rapidly towards core collapse.

Figure 2: The neutrino luminosity produced by the nuclear reaction and thermal
processes (the top two panels). The colorful solid lines are the ratio of the
nuclear reaction neutrino luminosity to the thermal neutrino luminosity
plotted along stellar evolution tracks in a neutrino H-R diagram. Gray curves
indicate where nuclear reaction neutrinos dominate, green curves where thermal
neutrinos dominate, and blue and red curves where the reaction and thermal
neutrino luminosities are within a factor of 10. The bottom two panels show
the ratio of the neutrino luminosity to the photon luminosity during different
evolutionary stages. The colorful solid lines give the ratio of the neutrino
luminosity to the photon luminosity in a neutrino H-R diagram. The left panels
are for 20 to 50 $\rm M_{\odot}$ models with non-rotation and high metallicity
(($Z=0.014$)), the right panels are for 60 to 90 $\rm M_{\odot}$ models with
non-rotation and solar metallicity. The positions of the neutrino luminosity
at the ignition of He, C, O, and Si are given.

As mentioned in the above section, both nuclear reaction and thermal process
can produce neutrinos. The upper two panels of Figure 2 depict the ratio of
nuclear reaction neutrino luminosity to thermal process neutrino luminosity,
traced along the stellar evolution trajectories in a neutrino H-R diagram. The
lower panels of Figure 2 illustrate the ratio of neutrino luminosity to photon
luminosity, also traced along the stellar evolution trajectories in a neutrino
H-R diagram. During hydrogen and helium burning, the majority of neutrinos are
generated via nuclear reactions, and massive stars primarily radiate their
energy through photons. Upon the commencement of carbon burning, the dominant
source of neutrinos transitions from nuclear reactions to thermal processes.
Consequently, neutrinos become the principal conduit for energy loss, and the
neutrino luminosity significantly surpasses the photon luminosity.

In particular, the neutrino luminosity produced by the nuclear reactions
becomes gradually important during Si burning. The main reason is given in
Figure 3, which illustrates the change of central electron mole number
$Y_{\mathrm{e}}$ in a 20 $\rm M_{\odot}$ star with non-rotation and high
metallicity ($Z=0.014$) after oxygen depletion. It shows that the number of
electrons in the star decreases during silicon burning. This reduction can be
attributed to an increase in the internal density and temperature, which leads
to an increase in the electron degeneracy pressure. As a result, the mean
energy of electrons becomes larger than that of many heavy nuclei, increasing
the rate of $\rm e^{-}$ capture processes in the stellar interior. Once the
core temperature of the massive star reaches $\rm 5\times 10^{9}\ K$, nuclear
statistical equilibrium (NSE) is reached, and the $\rm e^{-}$ capture process
is significantly enhanced (Kato et al., 2017). During the NSE phase, $\beta$
decay and electron capture (EC) on heavy nuclei dominate the production of
neutrinos. As demonstrated in the top rows of Figure 2, once the core begins
to collapse, a large number of $\rm e^{-}$ captures enable nuclear reaction
processes to produce far more neutrinos than thermal processes.

Figure 3: Evolution of $Y_{e}$ in the center of a 20 $\rm M_{\odot}$ star with
non-rotation and high metallicity ($Z=0.014$) after Oxygen depletion. The red
dot represents Si ignition, the blue dot represents Si depletion.

###  3.2 Effect of metallicity on neutrino luminosity

In order to discuss its effects on neutrino luminosity, we present Figure 4,
which shows the neutrino H-R diagram of massive stars with different
metallicities. In this section, we do not take into account the effect of
rotation on massive stars. For the low-metallicity models, the neutrino
luminosity is higher than that of models with high metallicity ($Z=0.014$)
during MS because they have a higher central temperature. Simultaneously, the
low-metallicity model has a weaker mass-loss rate, which results in much
hydrogen shell, and heavier helium core. Therefore, the massive stars with a
low metallicity have higher neutrino luminosity during helium burning.

Figure 4: The neutrino H-R diagram. The left panel is 20 $\rm M_{\odot}$ star
with non-rotation, right panel is 50 $\rm M_{\odot}$ star with non-rotation.
The red line is $\rm Z=0.014$ model, green line is $\rm Z=0.001$ model, blue
line is $\rm Z=0.0001$ model. The ZAMS is located on the dashed line. The
positions of the neutrino luminosity at the ignition of He, C, O, and Si are
given.

###  3.3 Effect of rotation on neutrino luminosity

Rotation is widely acknowledged to profoundly influence stellar structure and
evolution. Here, taking the models with 20, 50 and 90 M⊙ as examples, we
discuss its effects on the neutrino luminosity.

Figure 5: The neutrino H-R diagram. From left to right is 20, 50 and 90 $\rm
M_{\odot}$ star, respectively. The red line is high metallicity ($Z=0.014$)
without rotation model, green line is high metallicity with high rotation
($\rm\omega/\omega_{crit}=0.7$) model, blue line is $\rm Z=0.001$ with high
rotation model. The ZAMS is located on the dashed line. The positions of the
neutrino luminosity at the ignition of He, C, O, and Si are given. Figure 6:
The Kippenhahn diagrams of 20 $\rm M_{\odot}$ star. The green shaded areas
indicate convection zones, colors indicate the specific energy rate due to
neutrino emission $\rm\epsilon_{\nu}$. From left to right are $\rm Z=0.014$
and $\omega/\omega_{\mathrm{crit}}=0$, $\rm Z=0.014$ and
$\omega/\omega_{\mathrm{crit}}=0.7$, $\rm Z=0.001$ and
$\omega/\omega_{\mathrm{crit}}=0.7$ models.

Figure 5 displays 20, 50 and 90 $\rm M_{\odot}$ stellar evolutionary tracks
with different rotation and metallicity in the neutrino H-R diagram. At high
metallicity ($Z=0.014$), it is obvious that the neutrino luminosity increases
during MS for 20 $\rm M_{\odot}$ star with high rotation
($\omega/\omega_{\mathrm{crit}}=0.7$). As described by Georgy et al. (2013),
rotation induces mixing and brings light elements, like hydrogen, into the
core. This leads to larger convective and burning regions during MS, as shown
in middle panel of Figure 6, which increases the MS lifetime and enhances the
energy yield of the CNO cycle. For massive stars, neutrinos are mainly
produced by the CNO cycle during hydrogen burning. Therefore, the more
efficient the CNO cycle inside the star, the higher the neutrino luminosity
during hydrogen burning. The high rotation model
($\omega/\omega_{\mathrm{crit}}=0.7$) with high metallicity ($Z=0.014$) has a
0.4 dex larger neutrino luminosity than the non-rotating model during hydrogen
burning. At low metallicity, the surface angular momentum loss is reduced due
to diminished stellar wind losses. In our simulations, the 20 M⊙ model with
rapid rotation ($\omega/\omega_{\mathrm{crit}}=0.7$) can undergo a strong
mixing process and result in efficient chemically homogeneous evolution during
the MS phase. The right panel of Figure 6 illustrates that the convection and
burning regions of the models with low metallicity and rapid rotation are
larger than models with high metallicity and rapid rotation. This results in a
more efficient CNO cycle and, therefore, a higher neutrino luminosity.

Regarding the neutrino HR diagram of 50 and 90 $\rm M_{\odot}$ stars depicted
in Figure 5, at high metallicity ($Z=0.014$) and high rotation
($\omega/\omega_{\mathrm{crit}}=0.7$), these stars can undergo intense star
wind losses and angular momentum loss during its MS lifetime. Therefore,
despite the high rotation rate, mixing is not efficient, and there is no
significant increase in neutrino luminosity compared to a 20 $\rm M_{\odot}$
star. At low metallicity and high rotation, the evolution of 50 and 90 $\rm
M_{\odot}$ stars follow a trajectory similar to that of a 20 $\rm M_{\odot}$
star. The smaller mass loss and more massive core allow these stars also
undergo quasi-chemical homogeneous evolution. Therefore, the neutrino
luminosity of the 50 and 90 $\rm M_{\odot}$ stars increases significantly
during the MS. In fact, the rotation rate required for 50 and 90 $\rm
M_{\odot}$ star to achieve chemically homogeneous evolution is smaller than
that for a 20 $\rm M_{\odot}$ star (Yoon et al., 2006).

During the helium burning, neutrino H-R diagram of Figure 5 displays the high
rotation ($\omega/\omega_{\mathrm{crit}}=0.7$) models show more significant
drops in neutrino luminosity than non-rotation. The drop occurs because mixing
processes take away significant amounts of hydrogen from the stellar surface,
leaving a hydrogen shell that is too thin to produce more neutrinos. This
phenomenon is particularly pronounced in low-metallicity models, as the star
undergoes quasi-chemical homogeneous evolution and completely loses its
hydrogen shell, causing the neutrino luminosity to drop almost linearly during
the helium burning. Simultaneously, during the helium burning, masstive stars
with high rotation lose more hydrogen shell. For more massive stars, such as
our 50 and 90 $\rm M_{\odot}$ star models, their helium cores can even be
eroded by intense star wind losses.

Figure 7: The neutrino luminosity, center temperature($T_{c}$) and energy
yield of electron pair annihilation (from top to bottom) during C (left) and O
(right) burning. The solid line is $\rm Z=0.014$ and
$\omega/\omega_{\mathrm{crit}}=0$ star model, the dashdot line is $\rm
Z=0.014$ and $\omega/\omega_{\mathrm{crit}}=0.7$ star model, the dotted line
is $\rm Z=0.001$ and $\omega/\omega_{\mathrm{crit}}=0.7$ star model. The red
line is 20 $\rm M_{\odot}$, the green line is 50 $\rm M_{\odot}$, the blue
line is 90 $\rm M_{\odot}$.

Rotation and metallicity can also affect the neutrino luminosity at the
advanced burning stage of massive stars. The top panel of Figure 7 shows the
neutrino luminosity of 20 $\rm M_{\odot}$, 50 $\rm M_{\odot}$, and 90 $\rm
M_{\odot}$ stars with different rotation and metallicity during carbon and
oxygen burning. For the 20 $\rm M_{\odot}$ star, high rotation can increase
the neutrino luminosity, with this increase being more pronounced at low
metallicity. This is because the 20 $\rm M_{\odot}$ star with high rotation
($\omega/\omega_{\mathrm{crit}}=0.7$) can have a larger $\mathrm{He}$ and
$\mathrm{CO}$ core mass than those without rotation, and the lower the
metallicity, the larger the increase in mass. Additionally, massive stars with
low metallicity show a higher temperature after helium burning, as shown in
the middle panel of Figure 7. In section 3.1, it was mentioned that the
neutrino luminosity after helium burning to oxygen depletion is mainly
determined by the thermal process, particularly pair annihilation. It is
independent of density and sensitive to temperature, as indicated by Itoh et
al. (1996), so an increase in neutrino luminosity can occur significantly, as
shown in the bottom panel of Figure 7. These factors contribute to the
increase in neutrino luminosity at later stages of stellar evolution. Although
the $T_{\mathrm{c}}$ of a 20 $\rm M_{\odot}$ star with high rotation
($\omega/\omega_{\mathrm{crit}}=0.7$) is slightly lower than that of a non-
rotating star during carbon burning, the larger CO core mass leads to a higher
neutrino luminosity.

For the 50 $\rm M_{\odot}$ star with high rotation
($\omega/\omega_{\mathrm{crit}}=0.7$) and high metallicity ($Z=0.014$) in
Figure 7, the neutrino luminosity do not increase significantly during carbon
and oxygen burning, similar to the MS. This is because as the initial mass of
the massive star increases at high metallicity, the effect of rotation on
$M_{\mathrm{CO}}$ gradually decreases, as shown in Figure 8. For massive stars
with $\rm M>50\ M_{\odot}$, the global effect of mass loss surpasses the
influence of rotation, which even decreases the CO core mass (Limongi &
Chieffi, 2018). For example, as our 90 $\rm M_{\odot}$ star (which is shown by
blue lines in Figure 7) with high metallicity model, high rotation
($\omega/\omega_{crit}=0.7$) reduce the CO core mass by $32\%$ compared to
non-rotation. As a result, the average neutrino luminosity is always reduced
by about $60\%$ and $95\%$ during carbon and oxygen burning, respectively,
compared to the 90 $\rm M_{\odot}$ non-rotation and high metallicity model. At
low metallicity, the effect of mass loss can be neglected. Both 50 and 90 $\rm
M_{\odot}$ star models with high rotation
($\omega/\omega_{\mathrm{crit}}=0.7$) have a larger CO core mass than non-
rotating models after helium depletion, resulting in a significant increase in
neutrino luminosity during carbon and oxygen burning.

It is worth noting that our 90 $\rm M_{\odot}$ star model with high rotation
and low metallicity can enter the pulsation pair instability regime during O
burning. As described by Woosley et al. (2002), massive stars with a CO core
mass exceeding 35 $\rm M_{\odot}$ have a structural adiabatic index
$\Gamma_{1}$ below $4/3$, which can lead to instability. To avoid this
instability, we stop the calculation as soon as the instability sets in.

Figure 8: $\rm M_{CO}-M_{INI}$ relation. The dark line is initial mass. Figure
9: Similarly to Figure 2, but for models with different network. The left
panel is for approx21.net, while the right panel is for $mesa\\_204.net$,
respectively.

###  3.4 Effect of network on neutrino luminosity

The selection of nuclear network plays an important role in the nuclear
neutrino luminosity of massive stars. Figure 9 shows the ratio of thermal
process to nuclear process neutrino luminosity along the stellar evolution
track of 20 $\rm M_{\odot}$ star with non-rotation and high metallicity
($Z=0.014$) in a neutrino H-R diagram when adopted the basic network
approx21.net and large network $mesa\\_204.net$, respectively. Obviously, in a
small network, electron capture and $\beta$ decay contributes little to the
neutrino production due to the absence of a large number of isotopes, so that
nuclear neutrinos are still much smaller than thermal neutrinos after silicon
burning. In a larger network, however, it not only tracks a large number of
nuclear neutrinos after Si burning, but also shows a higher neutrino
luminosity.

###  3.5 The effects of neutrino magnetic moments on the evolution of massive
stars

Based on the above simulations, massive stars produce a huge number of
neutrinos during their lifetime. If neutrinos have a magnetic moment, their
effects on the evolution of massive stars may be significant. It is because
the NMM can cause additional energy loss when pair annihilation occurs inside
the star, as shown in Figure 10. The neutrino losses from pair annihilation
dominate the evolution of massive stars from He burning to oxygen depletion.
Therefore, we have evolved all the models in this paper with
$\rm\mu_{-12}=25$, where $\rm\mu_{-12}$ is the NMM in units of $\rm
10^{-12}\mu_{B}$.

Figure 10: The energy yield of electron pair annihilation during C (left) and
O (right) burning when the $\rm\mu_{-12}=0$ and $\rm\mu_{-12}=25$,
respectively. $\rm\mu_{-12}$ is the NMM in units of $\rm 10^{-12}\mu_{B}$, the
star mass is 20 M⊙.

When $\rm\mu_{-12}=25$, the effects of NMM on the evolution of massive stars
can be negligible. The energy loss caused by such a small NMM cannot change
stellar destiny of massive stars from ZAMS to FeCC. This result is consistent
with that of Heger et al. (2009), who examined the effect of NMM on stellar
evolution in the mass range of $7\ \mathrm{M_{\odot}}\lesssim M\lesssim 18\
\mathrm{M_{\odot}}$.

Figure 11: The change of $X_{{}^{12}\mathrm{C}}^{\mathrm{center}}$,
$M_{\mathrm{{CO}}}/\mathrm{M_{\odot}}$, and $\xi_{2.5}$ when the
$\rm\mu_{-12}=0$ and $\rm\mu_{-12}=25$, respectively. $\rm\mu_{-12}$ is the
NMM in units of $\rm 10^{-12}\mu_{B}$, FeCC is Fe Core Collapse,
$X_{{}^{12}\mathrm{C}}^{\mathrm{center}}$ is the faction of central
${}^{12}\mathrm{C}$ at the end of central helium depletion,
$M_{\mathrm{{CO}}}/\mathrm{M_{\odot}}$ is the CO core mass, $\xi_{2.5}$ is the
compactness parameter. According to O’Connor & Ott (2011), the dashed line at
$\xi_{2.5}=0.45$ serves as a demarcation, separating models that could
potentially undergo explosion (below the line) or implode (above the line).
The initial mass of stars ranges from 20 to 90 M⊙ with a mass interval of 10
M⊙.

For further study, we discuss the effect of NMM on the compactness of
presupernove stellar cores, which is usually expressed in $\xi_{\mathrm{M}}$.
It is defined by (O’Connor & Ott, 2011)

$\displaystyle\xi_{\mathrm{M}}=\frac{M/\mathrm{M}_{\odot}}{R\left(\mathrm{M}_{\text{bary
}}=M\right)/1000\mathrm{~{}km}},$ (15)

where $\rm M=2.5\ M_{\odot}$ is used for quantifing the mean density near the
iron core. In order to maintain substantial accuracy of compactness parameter,
we end the evolution of massive stars when the collapse speed first reaches
$\rm 1000\ km\ s^{-1}$ (Sukhbold & Woosley, 2014).

The extra cooling of massive stars can have an impact on the nuclear reactions
inside the star, ultimately affecting the mass of the final CO core, as stated
by (Heger et al., 2009). However, these changes are trivial for star with a
mass larger than 20 M⊙. Figure 11 shows the effects of NMM on
$X(^{12}\mathrm{C})$ after He depletion, CO core mass, and $\xi_{2.5}$. It
clearly explains that for stars with a mass larger than 20 M⊙, even after
considering additional cooling, the structure of massive stars does not
undergo significant changes, so the explodability of the massive stars remains
unchanged.

##  4 conclusions

Using MESA, we simulate the stellar structure and evolution for massive stars
with a mass from 20 to 90 M⊙, and calculate the neutrino luminosity and give
neutrino H-R diagram. Our results show: the greater the mass of the star, the
more intense the nuclear reactions inside it are, and the higher the neutrino
luminosity is. During helium burning, the neutrino luminosity decreases more
significantly as the initial mass increases due to mass loss. During the
hydrogen and helium burning stages, neutrinos are mainly produced by nuclear
processes. Compared with photons, neutrinos carry away less energy. After
carbon burning, neutrinos are mainly produced by thermal process (photo
neutrino, pair annihilation and plasmon decay), meanwhile, photons are
trapped, neutrinos dominate stellar energy loss. However, after Si burning,
the large number of weak processes allows nuclear reaction processes to
produce more neutrinos than thermal processes.

Low metallicity can facilitate a more condensed structure within the massive
star, leading to higher $T_{\mathrm{c}}$ and higher neutrino luminosities
during hydrogen burning. During helium burning, low mass loss due to low
metallicity can prevent a significant decrease in neutrino luminosity due to
neutrino release from the burning of a large number of hydrogen shells.

Rotation also affects the neutrino luminosity, which is determined by the
initial mass and metallicity. At the high metallicity ($Z=0.014$), for stars
with a mass lower than 40 M⊙, rotational mixing can increase neutrino
luminosity during hydrogen, C, and O burning. For stars with a mass higher
than 40 M⊙, the rotational mixing is less efficient due to mass loss, neutrino
luminosity does not increase much. For more massive stars, neutrino luminosity
even decreases during C and O burning, such as our 90 M⊙ star model. At low
metallicity, where mass loss is negligible, the neutrino luminosity of all
massive stars with high rotation can increase significantly during hydrogen,
C, and O burning.

The energy loss caused by the NMM does not have effects on the evolutionary
destiny of stars with a mass larger than 20 M⊙, and does not significant
change the compactness at FeCC. Therefore, the cooling via NMM has no effect
on the explodability of massive stars.

## Acknowledgements

This work received the generous support of the National Natural Science
Foundation of China, project No.U2031204, 12163005, 12288102, and 12263006,
the science research grants from the China Manned Space Project with NO.CMS-
CSST-2021-A10, the Natural Science Foundation of Xinjiang No. 2021D01C075,
2022D01D85, and 2022TSYCLJ0006.

## Data Availability

The data generated, analysed and presented in this study are available from
the corresponding authors on reasonable request.

## References

  * Abi et al. (2020) Abi B., et al., 2020, Journal of Instrumentation, 15, T08008
  * Abusleme et al. (2021) Abusleme A., et al., 2021, Chinese Physics C, 45, 023004
  * Albanese et al. (2021) Albanese V., et al., 2021, Journal of Instrumentation, 16, P08059
  * Angulo et al. (1999) Angulo C., et al., 1999, Nuclear Phys. A, 656, 3
  * Bahcall (1989) Bahcall J. N., 1989, Cambridge University Press,
  * Beacom et al. (2017) Beacom J. F., et al., 2017, Chinese Physics C, 41, 023002
  * Beda et al. (2013) Beda A. G., et al., 2013, Physics of Particles and Nuclei Letters, 10, 139
  * Bell et al. (2005) Bell N. F., Cirigliano V., Ramsey-Musolf M. J., Vogel P., Wise M. B., 2005, Phys. Rev. Lett., 95, 151802
  * Bellini et al. (2014) Bellini G., et al., 2014, Phys. Rev. D, 89, 112007
  * Bionta et al. (1987) Bionta R. M., et al., 1987, in Barger V., Halzen F., Marshak M., Olive K., eds, Neutrino Masses and Neutrino Astrophysics (Including Supernova 1987a). pp 63–77
  * Böhm-Vitense (1958) Böhm-Vitense E., 1958, Z. Astrophys., 46, 108
  * Brott et al. (2011) Brott I., et al., 2011, A&A, 530, A115
  * Chugunov et al. (2007) Chugunov A. I., Dewitt H. E., Yakovlev D. G., 2007, Phys. Rev. D, 76, 025028
  * Couch (2013) Couch S. M., 2013, ApJ, 765, 29
  * Cui et al. (2018) Cui Z., Wang Z., Zhu C., Lü G., Chen H., Han Z., 2018, PASP, 130, 084202
  * Cyburt et al. (2010) Cyburt R. H., et al., 2010, ApJS, 189, 240
  * Eguchi et al. (2003) Eguchi K., et al., 2003, Phys. Rev. Lett., 90, 021802
  * Farag et al. (2020) Farag E., Timmes F. X., Taylor M., Patton K. M., Farmer R., 2020, ApJ, 893, 133
  * Farmer et al. (2016) Farmer R., Fields C. E., Petermann I., Dessart L., Cantiello M., Paxton B., Timmes F. X., 2016, ApJS, 227, 22
  * Fujikawa & Shrock (1980) Fujikawa K., Shrock R. E., 1980, Phys. Rev. Lett., 45, 963
  * Fuller & Haxton (2022) Fuller G. M., Haxton W. C., 2022, arXiv e-prints, p. arXiv:2208.08050
  * Fuller et al. (1985) Fuller G. M., Fowler W. A., Newman M. J., 1985, ApJ, 293, 1
  * Georgy et al. (2013) Georgy C., Ekström S., Granada A., Meynet G., Mowlavi N., Eggenberger P., Maeder A., 2013, A&A, 553, A24
  * Hainich et al. (2015) Hainich R., Pasemann D., Todt H., Shenar T., Sander A., Hamann W. R., 2015, A&A, 581, A21
  * Hamann et al. (1995) Hamann W. R., Koesterke L., Wessolowski U., 1995, A&A, 299, 151
  * Heger et al. (2000) Heger A., Langer N., Woosley S. E., 2000, ApJ, 528, 368
  * Heger et al. (2009) Heger A., Friedland A., Giannotti M., Cirigliano V., 2009, ApJ, 696, 608
  * Ikeda et al. (2007) Ikeda M., et al., 2007, ApJ, 669, 519
  * Itoh et al. (1996) Itoh N., Hayashi H., Nishikawa A., Kohyama Y., 1996, ApJS, 102, 411
  * Janka (2017) Janka H.-T., 2017, Handbook of Supernovae, p. 1095
  * Kato et al. (2017) Kato C., Nagakura H., Furusawa S., Takahashi K., Umeda H., Yoshida T., Ishidoshiro K., Yamada S., 2017, ApJ, 848, 48
  * Kato et al. (2020) Kato C., Hirai R., Nagakura H., 2020, MNRAS, 496, 3961
  * Kippenhahn et al. (2013) Kippenhahn R., Weigert A., Weiss A., 2013, Peking University Press,
  * Kutschera et al. (2009a) Kutschera M., Odrzywołek A., Misiaszek M., 2009a, Acta Physica Polonica B, 40, 3063
  * Kutschera et al. (2009b) Kutschera M., Odrzywołek A., Misiaszek M., 2009b, Acta Physica Polonica B, 40, 3063
  * Langanke & Martínez-Pinedo (2000) Langanke K., Martínez-Pinedo G., 2000, Nuclear Phys. A, 673, 481
  * Langer (1991) Langer N., 1991, A&A, 252, 669
  * Langer (1998) Langer N., 1998, A&A, 329, 551
  * Langer et al. (1983) Langer N., Fricke K. J., Sugimoto D., 1983, A&A, 126, 207
  * Lee & Shrock (1977) Lee B. W., Shrock R. E., 1977, Phys. Rev. D, 16, 1444
  * Limongi & Chieffi (2018) Limongi M., Chieffi A., 2018, ApJS, 237, 13
  * Lü et al. (2017) Lü G., Zhu C., Wang Z., Iminniyaz H., 2017, ApJ, 847, 62
  * Maeder (1980) Maeder A., 1980, A&A, 92, 101
  * Maeder & Meynet (2000) Maeder A., Meynet G., 2000, A&A, 361, 159
  * Marciano & Sanda (1977) Marciano W. J., Sanda A. I., 1977, Physics Letters B, 67, 303
  * Marti Magro & Super-Kamiokande Collaboration (2022) Marti Magro L., Super-Kamiokande Collaboration 2022, in 37th International Cosmic Ray Conference. p. 1088, doi:10.22323/1.395.01088
  * Martins et al. (2009) Martins F., Hillier D. J., Bouret J. C., Depagne E., Foellmi C., Marchenko S., Moffat A. F., 2009, A&A, 495, 257
  * Martins et al. (2013) Martins F., Depagne E., Russeil D., Mahy L., 2013, A&A, 554, A23
  * Mirizzi et al. (2016) Mirizzi A., Tamborra I., Janka H. T., Saviano N., Scholberg K., Bollig R., Hüdepohl L., Chakraborty S., 2016, Nuovo Cimento Rivista Serie, 39, 1
  * Mori et al. (2020) Mori K., Balantekin A. B., Kajino T., Famiano M. A., 2020, ApJ, 901, 115
  * Nagakura & Vartanyan (2022) Nagakura H., Vartanyan D., 2022, MNRAS, 512, 2806
  * O’Connor & Ott (2011) O’Connor E., Ott C. D., 2011, ApJ, 730, 70
  * Oda et al. (1994) Oda T., Hino M., Muto K., Takahara M., Sato K., 1994, Atomic Data and Nuclear Data Tables, 56, 231
  * Odrzywolek & Heger (2010) Odrzywolek A., Heger A., 2010, Acta Physica Polonica B, 41
  * Odrzywolek et al. (2004) Odrzywolek A., Misiaszek M., Kutschera M., 2004, Astroparticle Physics, 21, 303
  * Patton et al. (2017) Patton K. M., Lunardini C., Farmer R. J., Timmes F. X., 2017, ApJ, 851, 6
  * Paxton (2019) Paxton B., 2019, ApJS,
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton B., et al., 2018, ApJS, 234, 34
  * Raffelt (1996) Raffelt G. G., 1996, Stars as laboratories for fundamental physics : the astrophysics of neutrinos, axions, and other weakly interacting particles
  * Raffelt (2012) Raffelt G., 2012, Physics
  * Raj et al. (2020) Raj N., Takhistov V., Witte S. J., 2020, Phys. Rev. D, 101, 043008
  * Schootemeijer & Langer (2018) Schootemeijer A., Langer N., 2018, A&A, 611, A75
  * Shi et al. (2020) Shi Y., Xue X., Zhu C.-H., Wang Z.-J., Liu H.-L., Li L., Lü G.-L., 2020, Research in Astronomy and Astrophysics, 20, 005
  * Sukhbold (2016) Sukhbold T., 2016, PhD thesis, University of California, Santa Cruz
  * Sukhbold & Woosley (2014) Sukhbold T., Woosley S. E., 2014, ApJ, 783, 10
  * Vink et al. (2001) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574
  * Wellstein et al. (2001) Wellstein S., Langer N., Braun H., 2001, A&A, 369, 939
  * Woosley et al. (2002) Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015
  * Yoon et al. (2006) Yoon S. C., Langer N., Norman C., 2006, A&A, 460, 199