\section*{Results}
\subsection*{Degenerate vs Non-Degenerate cores}
A dislocation is characterised by its Burgers vector \(\mathbf{b}\) and line direction \(\boldsymbol{\xi}\)~\cite{anderson_2017_disl}.
On a continuum level, \(\mathbf{b}\) is defined as \(\oint_C \left(\partial \mathbf{u}/\partial l \right) \text{d} l\), where \(\mathbf{u}\) is the displacement vector and the integral \(\text{d}l\) is along an arbitrary closed path \(C\) enclosing the dislocation core.
On a discrete level, the contour integral is replaced by a summation of discrete segments \(k\) along a closed path,
\begin{equation}
\label{eqn:dd_dfn}
\mathbf{b} = \sum_{k=1}^N  \left( \mathbf{u}^{(k)}_{ij} - \mathbf{U}^{(k)}_{ij} \right)=\sum_{k=1}^N \mathbf{d}_k,
\end{equation}
where \(\mathbf{u}_{ij}\) and \(\mathbf{U}_{ij}\) are  displacement vectors between atom \(i\) and \(j\) before and after the introduction of the dislocation and \(\mathbf{d}_k\) is the differential displacement (DD) between atoms \(i\) and \(j\).
Each \(\mathbf{d}_k\) contributes a fraction of the total Burgers vector \(\mathbf{b}\)~\cite{anderson_2017_disl}.

In the BCC structure, the  screw dislocation Burgers vector \(\mathbf{b}\) can split into three \(\mathbf{d}_k\) enclosing the core.
Figure~\ref{fig:gamma_disl_Li_Ta_W} shows the core structures and \(\gamma\)-surfaces (energy  associated with translating one part of a crystal with respect to another along a crystal plane) of BCC Li, Ta and W calculated using DFT~\cite{supp_mat}.
Lithium exhibits a D-core while TMs Ta and W show ND-cores.
The core structures are dictated by a core energy \(E_\text{c}\), which is primarily associated with the shear displacements \(\mathbf{d}_1\),  \(\mathbf{d}_2\), \(\mathbf{d}_3\) in the \(\langle 111\rangle\) direction between neighbouring atoms enclosing the core centre (orange arrows in Fig.~\ref{fig:gamma_disl_Li_Ta_W}a-c).
For Li, Ta and W, the current DFT calculations show equal splitting of \(\mathbf{d}_1 \cong \mathbf{d}_2 \cong \mathbf{d}_3 \cong \mathbf{b}/3\) on three equivalent \(\{110\}\) planes at the core centre, i.e., all core centres are identical (Fig.~\ref{fig:gamma_disl_Li_Ta_W}).
This is consistent with nearly all previous BCC TMs core structure simulation results,  irrespective of the description of atomic interactions~\cite{rodney_2017_am} (DFT or interatomic potentials).
This core-splitting implies that the core centre (screw components) fulfils the \(\langle 110 \rangle\)-dyad symmetry of the BCC lattice and is thus ND.
Nevertheless, examination of the \(\gamma\)-surface shows that between the two ground states, the minimum energy path (MEP, black line in Fig.~\ref{fig:gamma_disl_Li_Ta_W}d-f) deviates from \(\langle 111 \rangle\) (white line).
The core energy can thus be further reduced by core splitting along the MEP with some edge components in \(\mathbf{d}_k\).

The principal differences between the D-/ND-cores lie in the details of the \(\mathbf{d}_k\) outside the core centre.
Applying the topological requirement of Eq.~\eqref{eqn:dd_dfn}, each \(\mathbf{d}_k\) must be compensated by two additional DD vectors in the triangular loop formed by one orange and two blue arrows (i.e., \(\sum \mathbf{d}_k=0\)).
The D-core (Fig.~\ref{fig:gamma_disl_Li_Ta_W}a) resolves each \(\mathbf{d}_k\) by a DD vector \(\mathbf{d}^\prime_k\) (blue arrow) of \textasciitilde{}\(\mathbf{b}/3\) along a \(\{110\}\) plane and a small DD (not visible in Fig.~\ref{fig:gamma_disl_Li_Ta_W}a).
For the ND-cores (Fig.~\ref{fig:gamma_disl_Li_Ta_W}b,c), each \(\mathbf{d}_k\) is compensated by two identical \(\mathbf{d}^{\prime\prime}_k\) (blue arrows) of \textasciitilde{}\(\mathbf{b}/6\) on two equivalent \(\{110\}\) planes. The D-/ND-core preference thus appears to be determined by the relative energy costs of the structures outside the core centre: one with a single \(\mathbf{d}^\prime_k\approx\mathbf{b}/3\) and the other with a pair of \(\mathbf{d}^{\prime\prime}_k\approx\mathbf{b}/6\).  These energy costs can be approximated by the stacking fault energy \(\gamma\) at \(\mathbf{b}/3\) and \(\mathbf{b}/6\); the ND core structure is preferred when \(2\gamma(\mathbf{b}/6) <\gamma(\mathbf{b}/3)\), as a well-known \(\gamma\)-line criterion~\cite{duesbery_1998_am}.

\begin{figure*}[!htbp]
  \centering
  \includegraphics[width=0.7\textwidth]{fig_2.pdf}
  \caption{\label{fig:prop_core_type} \textbf{Dislocation core structure versus proposed governing criteria}. \textbf{a} \(\gamma\)-line criterion. \textbf{b} Material index \(\chi\)-criterion.
 Dislocation cores transform from non-degenerate (filled symbols) to degenerate (empty symbols) as \(\chi\) drops below \textasciitilde{}0.65. }
\end{figure*}

The above \(\gamma\)-line based criterion provides some rationalisation in pure BCC TMs~\cite{duesbery_1998_am,rodney_2017_am}.  However, the DFT-based result of Li violates this criterion (Supplementary Fig.~\ref{fig:gamma_line_li_ta_w}), so do other cases calculated by interatomic potentials~\cite{fellinger_2010_prb} or DFT~\cite{romaner_2010_prl}.  To determine the true mechanism governing the cores, we calculate the core structures of 7 BCC TMs (V, Nb, Ta, Cr, Mo, W and Fe) based upon 72 different interatomic potentials (Supplementary Table~\ref{tab:potential}).  While the accuracies of the potentials vary, we may view them as 72 different types of atoms/average alloys with properties perturbed from respective pure BCC TM elements.  These potentials yield both D- and ND-cores, thus providing a broad examination over the proposed \(\gamma\)-line and \(\Delta E\) criteria.  To consolidate data across all elements, we use the material index \(\chi\) by normalising \(\Delta E\) of interatomic potentials with the intrinsic \(\Delta E^\text{P}\) of the pure element calculated by DFT.

Figure~\ref{fig:prop_core_type}a shows that the $\gamma$-line criterion does not distinguish different core structures.  No correlation is seen between the core structure and the ratio \(2\gamma(\mathbf{b}/6)/\gamma(\mathbf{b}/3)\), or element type, or interatomic potential formalism.  Rather, Fig.~\ref{fig:prop_core_type}b shows that the material index \(\chi\) provides a prediction of core structure.
More specifically, D-cores are seen when \(\chi < 0.45\) and ND-cores are observed for \(\chi>0.65\). When \(\chi\) approaches 1, the intrinsic value of pure TMs, all potentials yield the ND core in agreement  with DFT calculations~\cite{rodney_2017_am}.
This simple \(\chi\)-criterion is highly predictive, failing in only 3 of 72 cases examined~\cite{supp_mat}.
Near the threshold \( 0.45 < \chi < 0.65\), both D- and ND-cores are seen, indicating that other subtle factors (e.g., elastic constants or the shapes of the \(\gamma\)-surfaces) may play a role.
The \(\chi\)-criterion thus appears to be broadly applicable (DFT validations are provided below).
While empirical, it is rooted in core structure geometry~\cite{kroupa_1964_cjpb}, as illustrated below.

\begin{figure*}[!htbp]
  \centering
  \includegraphics[width=0.8\textwidth]{fig_3.pdf}
  \caption{\label{fig:bcc_fcc} \textbf{Atomic plane stacking sequences in BCC, HCP and FCC structures and EADs of screw dislocations.}
\textbf{a-c} Atomic stacking of BCC \(\{110\}\), HCP \(\{0001\}\) and FCC \(\{111\}\) planes.
\textbf{d-e} Atomic structure of the screw dislocation core viewed from two directions.
Displacements of atoms on the second A-layer in the \([110]\) direction on the \((1\bar{1}0)\) plane (orange arrows) creates a quasi-ABC stacking, similar to that for  \(\{111\}\) planes in  FCC.
\textbf{f-i} EADs (arrows) in cores of Li-MEAM and W-EAM.  Atoms are coloured based on their positions in the \(\langle  111 \rangle\) direction.}
\end{figure*}

\subsection*{Screw Dislocation Core Geometry}
On BCC \(\{110\}\) planes,  the atomic stacking follows an ABAB stacking sequence, similar to the ABAB close-packing of  \(\{0001\}\) planes in hexagonal close-packed (HCP) structures (Fig.~\ref{fig:bcc_fcc}a-c).
This stacking can be converted to quasi close-packing if atoms in the B-layer are displaced to reside above the centroids of the triangles formed by the A-layer atom triplets (blue dashed triangle in Fig.~\ref{fig:bcc_fcc}a).
The required displacements are along the \(\langle 110 \rangle\)-direction on \(\{110\}\) planes (orange arrows in Fig.~\ref{fig:bcc_fcc}e) and can be provided by non-screw components of the dislocation displacement field.
Each fractional Burgers vector \(\mathbf{d}_k\) has an excess displacement relative to the elastic field~\cite{supp_mat}; the details vary between the D- and ND-cores.
We illustrate these excessive atomic displacements (EADs) (Fig.~\ref{fig:bcc_fcc}, f-i) in the D-core of Li computed by an MEAM potential~\cite{supp_mat} and in the ND-core of W described by an EAM potential~\cite{setyawan_2018_jap}.
In the D-core of Li (Fig.~\ref{fig:bcc_fcc}f), the EADs mainly shift 6 atoms (3 at the core centre and 3 neighbouring them) in \(\langle 110 \rangle\) directions on three \(\{110\}\) planes.
Taking any pair of atoms in the pink shaded box as an example (Fig.~\ref{fig:bcc_fcc}f, others have similar EADs by symmetry), both atoms have EADs along \(\langle 110 \rangle\). These displacements  (arrows in Fig.~\ref{fig:bcc_fcc}f,g) move the two atoms on the A-layer (yellow atoms in Fig.~\ref{fig:bcc_fcc}e) towards the centroids of triangles formed by the B-layer atom triplets above them (black triangle), thus forming a local close-packed structure similar to that in FCC/HCP (\emph{cf.} Fig.~\ref{fig:bcc_fcc}b,c,e).
The displacements also put these two atoms in the C-layer position with respect to the A and B layers beneath them, forming a local quasi-ABC stacking sequence.

Kinematically, the EADs lead to similar differential displacements (\(\mathbf{d}_k \approx \mathbf{d}^{\prime}_k \approx \mathbf{b}/3\)) for the two atoms on the same \(\{110\}\) plane and hence the D-core structure.
Energetically, the resulting close packed, quasi-ABC stacking is favourable at the core of Li since \(\Delta E = - 0.8\) meV/atom.
In contrast, the EADs of W-EAM-ND (Fig.~\ref{fig:bcc_fcc}h,i) are relatively small and not aligned on \(\{110\}\) planes.
In this case, displacements in the \(\langle 110 \rangle \) directions on \(\{110\}\) planes incur a large energy penalty since the ABC close-packed structure is highly unfavourable (\(\Delta E = 384\) meV/atom).
Therefore, the D-/ND-core competition is determined by the energy cost associated with displacements along \(\langle  110 \rangle\) on \(\{110\}\) planes that lead to locally close packed structures.
Li and W represent extreme scenarios of the D- and ND-cores; similar examination of the Fe-GAP-ND~\cite{dragoni_2018_prm} core (\(\Delta E= 159\) meV/atom) and Fe-MEAM-D~\cite{asadi_2015_prb} core (\(\Delta E= 44\) meV/atom) shows that their EADs fall between these two cases (Supplementary Fig.~\ref{fig:fe_core_eecd}).

The structure formed at the D-core is not precisely FCC, but only similar in close-packing sequence; the geometric description here should be viewed as suggestive.
Nevertheless, the model explains the  trend in Fig.~\ref{fig:prop_core_type}, where the transition from ND- to D-core starts as \(\chi\) drops below a critical value \textasciitilde{}0.65 (towards smaller \(\Delta E\)).
Smaller \(\Delta E\) implies that FCC structure becomes less unfavourable (relative to BCC) and forming the quasi-ABC structure is increasingly favourable.
Our extensive study of core structures using different interatomic potentials indicates that relating the energy of the quasi-ABC structure to that of the perfect FCC structure is of sufficient accuracy to predict D- and ND-cores.
The current model, based on \(\Delta E\), is also consistent with DFT calculations for all BCC TMs performed to date~\cite{rodney_2017_am}.
For all BCC TMs,  \(\Delta E\) is large (from 138 to 483 meV/atom, Supplementary Fig.~\ref{fig:deltaE}) and  FCC structures are highly unfavourable, suggesting limited EADs in \(\langle 110 \rangle\) directions on \(\{110\}\) planes.
Hence all screw dislocations prefer the ND-core structure in pure BCC TMs (Fig.~\ref{fig:gamma_disl_Li_Ta_W}b,c and other DFT calculations~\cite{dezerald_2014_prb}).
On the other hand, alkaline metals have nearly zero \(\Delta E\) and show D-cores (Fig.~\ref{fig:gamma_disl_Li_Ta_W}a for Li-DFT and Supplementary Fig.~\ref{fig:alkali_core} for Li, Na and K-interatomic potentials).

\begin{figure*}[!htbp]
  \centering
  \includegraphics[width=0.75\textwidth]{fig_4.pdf}
  \caption{\label{fig:vca} \textbf{Prediction of \(\chi\) as a function of solute concentrations for binary W-TM alloys from VCA DFT calculations.}
  The solubility of each element (El) in W is denoted by \(\widehat{\text{El}}\) at the top. The elements are colour-coded by their group numbers in the periodic table.
  }
\end{figure*}

\begin{figure*}[!htbp]
  \centering
  \includegraphics[width=0.90\textwidth]{fig_5.pdf}
  \caption{\label{fig:chi_prop} \textbf{Prediction of key material properties as a function of \(\chi\) for binary W-TM alloys. } \textbf{a} ND and D cores. \textbf{b} Unstable stacking fault energy \(\gamma_{\text{us}}\). \textbf{c} Peierls barrier \(\Delta E_{\text{PB}}\) of the screw dislocation. All properties are calculated by VCA DFT. Lower \(\chi\) favours the degenerate core, reduces \(\gamma_{\text{us}}\) and \(\Delta E_{\text{PB}}\), all independent of solute types.
  }
\end{figure*}

\subsection*{Validation and Prediction in Transition Metal Alloys}
The above analysis provides a geometric link between the material index \(\chi\) and the screw dislocation core structure;  \(\chi\) can be thought of as a surrogate for local bonding in the core.
It serves as an indicator of core structure and provides a physical/thermodynamic basis for understanding core structure in all BCC metals.  Since slip occurs along the Burgers vector direction on \(\{110\}\) planes, the slip process can be related to a transformation between the BCC structure and the quasi-FCC structure (Fig.~\ref{fig:bcc_fcc}).  Therefore, the energy barriers associated with dislocation nucleation \(\gamma_\text{us}\) and glide (Peierls barrier) \(\Delta E_\text{PB}\) should also be related to the transformation energy barrier or \(\Delta E\)/\(\chi\). Solid solution alloying can effectively modulate \(\Delta E\) and thus be used as a means of controlling the screw dislocation core structure and slip behaviour.  The material index \(\chi\) thus has practical, engineering application for BCC TM alloy design. Since \(\chi\) depends on the energy difference of two elementary structures, DFT calculations can be used to predict \(\chi\) as a function of solute concentration \(c\) and identify favourable solute types/concentrations. We demonstrate this strategy in binary BCC TM alloys, and extensively in W-TMs, a well-studied system in experiments~\cite{stephens_1970_mmtb,li_2012_am} and simulations~\cite{romaner_2010_prl}.

We employ the DFT-based virtual crystal approximation (VCA~\cite{nordheim_1931_ap}) to predict \(\chi\) and core properties for alloys.
VCA replaces all atoms in an alloy by virtual atoms with a composition-weighted pseudopotential and valence electron number based upon the alloy constituent elements~\cite{bellaiche_2000_prb}.
While VCA accuracy is limited to average, homogeneous alloy properties, it retains a general first-principles DFT framework~\cite{supp_mat}.
Figure~\ref{fig:vca} shows the predicted \(\chi\) calculated for 28 binary W-TM alloys (Supplementary Fig.~\ref{fig:vca_mo_vasp} for Mo-TM alloys).
At low solute concentrations ($<$10 \%),  \(\chi\) increases for elements to the left of the solvent TM  in the periodic table and decreases sharply (favourably) for elements to its right.
For W/Mo, increasing solute valence electron concentration (VEC~\cite{stephens_1972_jlcm,klopp_1975_jlcm,guo_2011_jap}) decreases \(\chi\) and its slope.

Figure~\ref{fig:chi_prop} shows the core structure, \(\gamma_\text{us}\) and \(\Delta E_\text{PB}\) as a function of \(\chi\)  for binary W-TM alloys.
For all alloys, ND- and D-cores are seen when \(\chi\) is above and below $0.8$, similar to predictions by interatomic potentials (Fig.~\ref{fig:prop_core_type}b).
Both \(\gamma_\text{us}\) and \(\Delta E_\text{PB}\) scale  linearly with \(\chi\); i.e.,
\begin{align}
  \label{eqn:chi_scaling}
  \gamma_\text{us} (\chi) &= \gamma_\text{us}^\text{P}\left[1+ k_\text{us} (\chi - 1)\right]  \\
  \Delta E_\text{PB} (\chi) &= \Delta E_\text{PB}^\text{P} \left[1 +  k_\text{PB} (\chi - 1)\right] .
\end{align}
where \(k_\text{us} = 0.66\) and \(k_\text{PB} = 2.75\) as fit from the data.
The above expressions can be extrapolated to Group VI-based TM alloys.
Solutes to the left/right of the solvent increase/decrease \(\gamma_\text{us}\) and \(\Delta E_\text{PB}\), which is broadly consistent with solute hardening/softening effects measured in Group VI binary BCC TM alloys (Supplementary Fig.~\ref{fig:hs_gp_vi_exp}) and DFT-based calculations in Mo/W-TMs~\cite{trinkle_2005_science,hu_2017_am}.

The VCA calculations for W-Re show that \(\chi\) decreases with increasing Re for all concentrations.
At 20\% Re (\(\chi=0.73\)), the ND-core transforms to the D-core (Supplementary Fig.~\ref{fig:wre_core}), in agreement with the transformation seen at 25\% Re in explicit solute-based DFT calculations~\cite{romaner_2010_prl} and the observation of a change to alternating \(\{110\}\) slips
in experiments~\cite{li_2012_am}.
The core transformation is also seen in Mo-20\% Re (Supplementary Fig.~\ref{fig:more_core}).
In addition, reducing \(\chi\) can favourably decrease the barrier for dislocation nucleation (\(\gamma_\text{us}\)) at crack tips and lattice friction (\(\Delta E_\text{PB}\)) of the screw dislocation, the two limiting factors governing DBT in BCC TMs~\cite{gumbsch_1998_science}.  At sufficiently low \(\chi\) (e.g., 0.66, or 25 at.\% Re), the gap between the stress intensity factor \(K_\text{Ie}\) for dislocation emission and \(K_\text{Ic}\) for cleavage may be reduced to enable ductile behaviour in W at higher temperatures (Supplementary Fig.~\ref{fig:kie_kic_wre_calc})~\cite{supp_mat}.  In fact, increasing the Re concentration to  $\geq25$\% is the basis of many commercial W-Re alloys which exhibit high ductility and a DBT below room temperature.
For W-Ta alloys, \(\chi\) is always greater than 1 and no core transformation is predicted, in agreement with DFT calculations~\cite{li_2012_am} and experiments where a switch in slip behaviour is not observed~\cite{li_2012_am}.
Ta also increases \(\gamma_\text{us}\) and \(\Delta E_\text{PB}\) (Fig.~\ref{fig:chi_prop}b,c), and thus has deleterious effects on W alloy ductility, again consistent with experiments~\cite{wurster_2011_jnm}.  While Fig.~\ref{fig:vca} suggests many solutes (Tc, Fe, Os, Ru, Rh, Co, Ir etc.) possessing the ductilizing/softening effects of Re~\cite{klopp_1975_jlcm,caillard_2020_am}, their practical application is limited by their low solubility~\cite{thaddeus_1986_bapd} and/or cost.
Re is thus one of the most practical solutes for ductilizing W.
The present predictions provide explicit recommendations of solutes and concentrations that may be explored through non-equilibrium techniques (e.g., mechanical alloying or plasma sintering) to overcome solubility limits.
Overall, the \(\chi\)-based approach identifies solutes capable of ductilizing/softening Group VI TMs that are broadly consistent with nearly all experiments performed to date~\cite{ren_2018_ijrmhm} (Supplementary Fig.~\ref{fig:hs_gp_vi_exp}).

The applicability of the \(\chi\)-approach is further examined in binary TM alloys with solvents from Group V.
For V-TM alloys, solutes to the left of V in the periodic table (decreasing VEC) reduce \(\chi\) (Supplementary Fig.~\ref{fig:vca_v_vasp}) and are expected to have softening effects, while solutes to its right increase \(\chi\) and lead to hardening.
The V-results are opposite to that in Group VI but fully consistent with hardening/softening measurements of Group V TMs (Supplementary Fig.~\ref{fig:hs_gp_v_exp}), as well as DFT calculations based upon inter-string potential~\cite{medvedeva_2007_prb} and a double-kink nucleation~\cite{zhang_2015_sm} models.
As another example, Ti was observed to ductilize V at 77 K~\cite{fraser_1962_cmq}, in agreement with  \(\chi\)-model predictions.
While solute/VEC effects vary among elements, the proposed material index \(\chi\) succinctly and quantitatively captures such effects of quantum-mechanics origin via inexpensive DFT calculations.

\section*{Discussion}
The importance of \(\chi\) for \(\gamma_\text{us}\) and \(\Delta E_\text{PB}\) and thus overall material plasticity and toughness is indisputable~\cite{gumbsch_1998_science}.
Effects of \(\chi\) on the screw dislocation core structures are more subtle but significant.
A change from the ND-core to D-core can have multiple effects on dislocation behaviour; the most important of which is the activation of additional slip systems.
The D-core slips on alternating \(\{110\}\) planes with net slip on multiple \(\{112\}\) planes, while the ND-core is highly restricted to a single \(\{110\}\) plane~\cite{vitek_2008_dis}.
A complete switching of slip plane only occurs at high solute concentrations (corresponding to \(\chi\)<0.8) and may be influenced by additional factors such as stress state and temperature~\cite{weinberger_2013_imr,caillard_2020_am}.
At moderate solute concentrations, both D- and ND-core dislocations can co-exist, depending on local solute concentrations.
Their co-existence activates 24 slip systems; i.e., twice the number of ND core-favoured \{110\} slip systems.
The importance of an increase in the number of active slip systems should not be underestimated.  The ductility of W depends on texture and grain size~\cite{zhang_2009_msea}; single crystals are generally ductile and polycrystals brittle.
Activation of additional slip systems can be critical; providing the flexibility required for strain compatibility at grain boundaries/junctions and reducing the texture-sensitivity of ductility in structural alloys.

At low temperature, pure BCC TMs tend to be brittle as slip is inhibited by high lattice friction and limited to \(\{110\}\) systems.
With increasing temperatures, the ductility of these pure TMs increases dramatically as \(\Delta E_\text{PB}\)/\(\gamma_\text{us}\) is reduced and additional slip planes (\(\{112\}\) and \(\{123\}\)) become thermally activated.
Such effects of entropy can be equivalently achieved by solute-addition changes to the core structure, \(\gamma_\text{us}\) and \(\Delta E_\text{PB}\), as captured by \(\chi\).
In other words, plastic deformation of pure TMs observed at high temperature can be achieved by alloying with ductilizing solutes (e.g., Re in W and Mo) at low temperature.
This implies that appropriate choice of solutes can increase ductility and lower the DBT temperature, as seen in W/Mo-Re alloys.

The \(\chi\)-based predictions are also fully consistent with general ductile/brittle behaviour in emerging compositionally complex alloys. Refractory high entropy alloy NbTaMoW can be considered as adding group V BCC TM to group VI TM, which raises \(\chi\). The ND core is expected to be prevalent, consistent with recent DFT calculations~\cite{yin_2020_npjcm}. NbTaMoW is also generally brittle at low temperatures.  In contrast, HfTiVNb can be considered as adding group IV elements to the group V TM, which favorably reduces \(\chi\).  Preliminary DFT calculations show that HfTiVNb adopts both the ND and D core, depending on the local solute environment. HfTiVNb possesses tensile ductility above 25\% at room temperature in experiments~\cite{an_2021_mh}.

Solutes also have local effects, which influence \(\chi\) and all core-related properties locally; e.g., dislocation lattice friction by modifying double-kink nucleation/propagation, dislocation cross-slip and interactions~\cite{trinkle_2005_science}.
These local effects vary between BCC TMs and should not be overlooked.  The \(\chi\)-approach does not replace established temperature-dependent solid-solution softening/hardening models~\cite{trinkle_2005_science,zhao_2018_msmse,ghafarollahi_2021_am}.
Instead, \(\chi\) offers a fundamental understanding of the physical basis for controlling dislocation behaviour and quantitatively predicts solute effects at low temperatures where softening/toughening is most needed.
Since \(\chi\) is built upon crystal geometry and bonding characteristics (captured through energy differences of two simple crystal structures), and is validated in extensive DFT calculations and across nearly all interatomic potentials available, it is expected to be general and can be used to rapidly identify favourable solutes in the entire family of BCC TM alloys (e.g., Nb/Fe-based alloys~\cite{romaner_2014_pml}).

In summary, we revealed the mechanistic origin of the critical screw dislocation properties and presented a general model with a fundamental material index \(\chi\) for predicting and controlling screw dislocation average properties in all BCC materials.  The \(\chi\)-model rationalises core structures in pure alkaline metals, BCC TMs and their alloys, and is quantitatively related to lattice friction and barrier for dislocation nucleation.  More importantly, \(\chi\) can be computed rapidly for any alloy composition using DFT, providing a practical computational approach for ductile and tough BCC TM alloy design.  This approach was tested in several binary BCC TM systems and, for example, correctly identifies Re as the most practical solute for ductilizing W and the appropriate Re concentration range.  The material index \(\chi\) does not predict ductility or toughness directly, but does, nevertheless, provide a fundamental constraint on the nature of plasticity in these materials. The proposed \(\chi\) approach is first-principles-based; hence it is applicable to a wide-range of complex alloys where ductility is controlled/limited by lattice friction, dislocation nucleation barrier or nature's choice of slip systems.

\newpage
\section*{Methods}
An extensive range of calculations are performed using various methods and models in the current work.  We also used two density functional theory (DFT) packages to overcome some limitations in determining some specific properties and to cross-validate the generality of the \(\chi\)-model.  We therefore first describe the general methods, followed by details of individual models.  The parameters in the general methods apply to all the calculations unless mentioned otherwise in the respective models.

\subsection*{Density-functional theory calculations using VASP}
\label{sec:dft_method}
First-principles calculations are performed within the DFT framework using the Vienna Ab initio Simulation Package (VASP~\cite{kresse_1996_prb,kresse_1999_prb}) and Quantum ESPRESSO (QE~\cite{giannozzi_2009_jpcm}).

In VASP, the exchange-correlation functional is described by the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE~\cite{perdew_1996_prl}) parameterization.  In the standard DFT calculations, the core electrons are replaced by the projector augmented wave (PAW~\cite{blochl_1994_prb}) pseudopotentials with the valence states shown in Supplementary Table~\ref{tab:dft_details}. A first-order Methfessel-Paxton method~\cite{methfessel_1989_prb} is used to smooth eigenstate occupancy.  The plane-wave cutoff energy, \(\Gamma\)-centered Monkhorst-Pack k-point sampling mesh~\cite{monkhorst_1976_prb} and the smearing parameters sigma are established through convergence tests of 2-atom BCC unit cells with a threshold energy variation \(\Delta E_\text{c} < 1\) meV/atom. The final converged parameters are shown in Supplementary Table~\ref{tab:dft_details}.

In the virtual crystal approximation (VCA~\cite{bellaiche_2000_prb}) DFT calculations, virtual atoms are created with effective cores and valence electrons constructed based on the alloy constituents. In particular, the virtual atoms have their core electrons replaced by pseudopotentials \(\phi_\text{virtual}\), which are obtained via mixing that of the solvent and solute atoms and weighted by their atomic fraction \(X\).  For example, for a binary \(\text{A}_{X} \text{B}_{1-X}\) alloy, a pseudopotential is created for the core electrons of the virtual atom representing the average alloy property, i.e.,
\begin{equation}\label{eqn:vca_weighted_core_ve}
  \phi_\text{virtual} =  X \phi_\text{A} + (1-X) \phi_\text{B}
\end{equation}
where \(\phi_{i}\) are the pseudopotentials of the respective atoms.   In addition, the valence electron numbers of virtual atoms are calculated according to Eqn.~\ref{eqn:vca_weighted_core_ve}. The plane-wave kinetic cutoff energy is set at 450 eV and a fine \(\Gamma\)-centered Monkhorst-Pack k-point  mesh~\cite{monkhorst_1976_prb} (k-spacing \(2\pi \times 0.02\)\AA\(^{-1}\)) is used for Brillouin zone integration.  A first order Methfessel-Paxton smearing with a width of 0.2 eV is used to smooth the partial electron occupancies. The valence states used for different elements are listed in Supplementary Table~\ref{tab:vs_state_vca_vasp_qe}.

\subsection*{Molecular dynamics/statics simulations using LAMMPS}
\label{sec:md_method}
Molecular dynamics and static relaxation simulations are performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS~\cite{plimpton_1995_jcp}). Structure optimization is performed using the conjugate gradient method.  Convergence is assumed when forces on all atoms drop below \(10^{-12}\) eV/\AA{} in calculations of \(\Delta E\) and \(\gamma\) surfaces,  and \(10^{-10}\) eV/\AA{} in calculations of dislocation core structures.  For the GAP potential for Fe~\cite{dragoni_2018_prm,bartok_2010_prl,bartok_2013_prb,bartok_2015_ijqc,gap_www}, convergence tolerance is relaxed to \(10^{-12}\) eV for energy variation and \(10^{-6}\) eV/\AA{} for atomic forces, due to its slow convergence and high computational cost.

\subsection*{Interatomic potential models for BCC metals}

The screw dislocation has been extensively studied in elementary BCC transition metals using DFT calculations. A rich set of interatomic potentials have also been developed using various formalisms over the past several decades.  In this work, we examined nearly all interatomic potentials (publically accessible) for BCC metals developed up to date.  These potentials are available from the National Institute of Standards and Technology (NIST) potential repository~\cite{nist_2013}, the Open Knowledgebase of Interatomic Models (OpenKIM) project~\cite{openkim_2011} or literature.  In total, we examined 72 potentials, as shown in Supplementary Table~\ref{tab:potential}.  These potentials use different functional forms in their descriptions of interatomic interactions. The choice made here is not biased toward any particular form, element or fitting procedure.

In addition, we developed two new MEAM potentials for BCC V (V1 and V2). The potential parameters are fitted with target properties shown in Supplementary Fig.~\ref{fig:meam_prop_v1_v2} using the particle swarm optimization algorithm.  The two fit-for-purpose potentials have nearly identical properties, except for \(\Delta E\), which is intentionally controlled at 91 meV/atom and 124 meV/atom, respectively.  With \(\Delta E\) from DFT at 243 meV/atom, V1 and V2 have \(\chi = 0.37 \) and 0.51, respectively. V1 adopts a D-core and V2 adopts a ND-core for the screw dislocation, consistent with the proposed \(\chi\)-criterion.  A new MEAM interatomic potential for Li is also developed with \(\Delta E\) = -0.8 meV/atom close to that from DFT.  All potential parameters for Li are shown in Supplementary Table~\ref{tab:meam_pot_param} and can be used in LAMMPS directly.  The MEAM potential for Li exhibits a D-core, in agreement with DFT and consistent with the geometric model where D-core is favored if the FCC structure is not hugely energetically unfavourable.

\subsection*{Calculation of energy difference \(\Delta E\) between FCC and BCC structures}
\label{sec:calc_chi_pure_ele}
The energy differences \(\Delta E\) per atom between the FCC and BCC structures are calculated using two-atom BCC and 4-atom FCC unit cells by DFT in VASP.   Convergence is assumed when energy variation drops below \(10^{-5}\) eV per electronic step and \(10^{-4}\) eV per ionic step.  The final \(\Delta E\) results are shown in Supplementary Table~\ref{tab:delta_E_bcc_fcc}.  Similarly, the energy difference \(\Delta E\) is also calculated using all the interatomic potentials in Supplementary Table~\ref{tab:potential} and convergence criterion described earlier. In all cases, the values of \(\Delta E\) do not depend on the employed supercell size within convergence tolerances.

\subsection*{Generalized stacking fault energy surface}
\label{sec:gamma_surf_calc_vasp}

The generalized stacking fault energy surfaces (\(\gamma\)-surface) of Li, Ta and W are calculated by DFT using VASP. A slab-supercell model is used, with its crystal orientation given in Supplementary Table~\ref{tab:supercell_dft}. The slab contains 12 \(\{110\}\)-plane atom layers and a 20 \AA{} vacuum layer.  \(\gamma\)-surfaces are calculated using the classical method by Vitek. For each specific stacking fault position, a homogeneous slip displacement \(\mathbf{s}\) is applied to atoms in the upper half block of the supercell to create the stacking fault between the sheared and un-sheared atom block.  Stacking fault energies are calculated with all atoms allowed to move in the direction perpendicular to the slip, \(\{110\}\) plane. The same convergence criteria as that for \(\Delta E\) are used for the \(\gamma\)-surface calculations.

In the calculations of \(\gamma\)-surface using interatomic potentials in LAMMPS, tilted supercells of dimensions \textasciitilde{}\(10 \times 10 \times 40\) (\AA{}\textsuperscript{3}) are used.  The stacking fault energies are calculated with atoms constrained to move in the direction perpendicular to the slip plane only (same as in DFT, see Ref.\cite{yin_2017_am} for details).

\subsection*{Dislocation core structure}
\label{sec:disl_core}

Dislocation core structures are calculated using both a dipole configuration (dipole method) and a single dislocation with fixed boundary conditions (cluster method). For the dipole method, a fully periodic prism supercell is chosen with the dislocations arranged in quadrupolar positions, as shown in Supplementary Fig.~\ref{fig:dis_sche}b.  The supercell vectors are
\begin{align}
  \mathbf{c}_1 &= \mathbf{e}_{1}, \\
  \mathbf{c}_2 &= 5\mathbf{e}_{2}, \\
  \mathbf{c}_3 &= \dfrac{1}{2}\mathbf{e}_{1} + \dfrac{5}{2}\mathbf{e}_{2}+9\mathbf{e}_{3},
\end{align}
where (\(\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3}) = (a/2[111],a[11\bar{2}],a[\bar{1}10])\) and  \(a\) is the lattice parameter of the BCC structure determined in unit cell calculations. In total, the supercell contains 135 atoms.

A dislocation dipole, i.e., a pair of dislocations with opposite Burgers vectors \(\pm \mathbf{b}\), is introduced at easy-core/hard-core positions (Supplementary Fig.~\ref{fig:dis_sche}c) using the Babel package developed by E. Clouet~\cite{clouet_2021_babel}.  In particular, the displacement field of arrays of dipoles (repeated 16 times in the \(\mathbf{c}_2\) and \(\mathbf{c}_3\) directions) is firstly applied to atoms in the supercell, followed by applying a homogeneous strain to accommodate the plastic strain introduced by the dipole (see Ref.~\cite{rodney_2017_am} for details).  During structure optimization, the cell vectors are fixed and all atoms are allowed to move until the ionic force is below \(5\times 10^{-3}\) eV/\AA. DFT calculations show that the easy-core is the minimum energy configuration and the hard-core is a local maximum. The DFT-computed core structures (ND vs D) are shown in Fig.~\ref{fig:gamma_disl_Li_Ta_W} for BCC Li, Ta and W.

The dipole approach allows modelling of dislocation cores with a fully periodic supercell of relatively small sizes.  It avoids the complexity of free surfaces, but introduces dislocation interactions and a non-negligible homogeneous strain scaling with supercell sizes.  Therefore, we also used a hexagonal prism supercell (Supplementary Fig.~\ref{fig:dis_sche}a) with a single dislocation and performed dislocation core structure calculations with interatomic potentials.  In this approach, supercells of radius \(R\) \textasciitilde{}92 \AA{} and length of one Burgers vector are first created with perfect BCC structures. The crystallographic \([11\bar{2}]\) direction is aligned with the \(x\)-axis, \([\bar{1}10]\) aligned with the \(y\)-axis and the screw dislocation line direction \([111]\) aligned with the \(z\)-axis, i.e., the hexagonal prism axis (Supplementary Fig.~\ref{fig:dis_sche}d).  The \(z\)-direction has periodic boundary conditions while the \(x\) and \(y\) directions are treated as surfaces.  This geometry keeps the \(\langle  111 \rangle\)-axis threefold symmetry of the BCC lattice.  Each supercell contains \textasciitilde{}3600 atoms (depending on element types).  A single dislocation is then introduced by displacing atoms according to the anisotropic elastic displacement field of the corresponding Volterra dislocation at the supercell center (easy-core position). Atoms within \textasciitilde{}12 \AA{} from the outer surface (grey area in Supplementary Fig \ref{fig:dis_sche}a) are always fixed to their elastic displacements.  The initial dislocation structure is first optimized using conjugate gradient algorithm at 0 K, then equilibrated for 100,000 time steps (100,000 femtoseconds) at finite temperatures (\textasciitilde{}50 K) and optimized again at 0 K.  This procedure is not always needed, but helps to obtain well-equilibrated core structures in some interatomic potentials. Only 0 K structure optimization is performed for the GAP potential for Fe~\cite{dragoni_2018_prm}.  The obtained core structures (ND vs D) are summarized in Fig.~\ref{fig:prop_core_type}.

\subsection*{Differential displacement plot and atomistic visualisation}
Dislocation core structures are visualised using the differential displacement (DD) map between neighboring atoms~\cite{vitek_1970_pma}.  For screw components (Fig.~\ref{fig:gamma_disl_Li_Ta_W}, Supplementary Figs.~\ref{fig:prop_core_type_v} and~\ref{fig:alkali_core}), the arrows are scaled so that the largest components (i.e., \textasciitilde{}\(\mathbf{b}/3\)) touch the neighboring atoms.

Excessive atomic displacements (EADs) (Fig.~\ref{fig:bcc_fcc}f-i and Supplementary Fig.~\ref{fig:fe_core_eecd}) are obtained by using the corresponding anisotropic Volterra dislocation field as references. The EADs of ND-/D-cores possess different symmetry; in the D-core, they have a \(\langle{111} \rangle\)-threefold screw axis symmetry while that in the ND-core have both the \(\langle{111} \rangle\)-threefold screw axis symmetry and the \(\langle {110} \rangle\)-diad axis symmetry. The EADs illustrate the key, subtle differences among different materials/models. They are generally very small, and are thus magnified by 6 times.

Atomic configurations (Fig.~\ref{fig:bcc_fcc}a-e) are visualised using the Open Visualisation Tool (OVITO~\cite{stukowski_2009_msmse}).

\subsection*{Calculation of \(\chi\) in binary alloys with DFT VCA in VASP }
\label{sec:chi_vca_vasp}

The values of \(\chi\) are calculated at a 2\%-increment for binary W-TM, Mo-TM and V-TM alloys using the VCA method in VASP.  The supercells and convergence criteria are the same as that used for calculating \(\Delta E\) in pure elements. The DFT settings are described in the Methods: Density-functional theory calculations using VASP, except that the plane-wave kinetic cutoff energy is increased to 550 eV for all considered elements.  We used the same cutoff energy for 84 alloys here to demonstrate the method of \(\chi\) to identify favourable solutes, as commonly employed in high throughput DFT calculations~\cite{lejaeghere_2016_science}.  Refined calculations can be performed once favourable solutes are identified. The final results are shown in Fig.~\ref{fig:vca}, Supplementary Figs.~\ref{fig:vca_mo_vasp} and ~\ref{fig:vca_v_vasp}.

\subsection*{Peierls barrier of binary W-TM alloys with DFT VCA in VASP}

Peierls barrier can be defined as the energy barrier per unit length for a straight dislocation to move from one energy valley (equilibrium position) to the next one at zero stress and zero temperature.  For the \(1/2\langle 111 \rangle\) screw dislocation in BCC structures, the two equilibrium positions correspond to the easy-core configuration in adjacent lattice positions (Supplementary Fig.~\ref{fig:dis_sche}c).  In the current work, we use the nudged elastic band (NEB~\cite{henkelman_2000_11_jcp}) method to calculate the dislocation migration path between these two adjacent easy-core configurations and thus obtain the height of the energy barrier, i.e., the Peierls barrier.  The transition path is calculated for pure W by DFT in VASP and for W-5\%Re, W-10\%Re and W-10\%Ta by DFT VCA in VASP.  The dipole supercell with 7 linearly interpolated replicas are used for all calculations (see Refs.~\cite{samolyuk_2013_jpcm,dezerald_2014_prb} for details).

Supplementary Fig.~\ref{fig:vca_wtms_pe} shows the energy variation along the transition path.  For all cases, the Peierls potential profiles have a near-sinusoidal profile with a peak in the middle replica corresponding to the saddle-core configuration.  Re and Ta respectively decreases and increases the Peierls barrier relative to that of pure W, in agreement with previous DFT calculations~\cite{samolyuk_2013_jpcm,dezerald_2014_prb}.

The NEB method requires simultaneous structure optimizations in all replicas and is computational expensive.  For calculations in a wide range of W-TM alloys, we employ the reaction coordinate method~\cite{xu_1998_cms} (i.e., the drag method) to compute their Peierls barriers as a function of solute concentrations and \(\chi\).  In the drag method, the 5 intermediate replicas are linearly interpolated between the initial and final configurations (easy-1 and easy-2 in Supplementary Fig.~\ref{fig:dis_sche}c). In the middle replica, a single atom near the core (the light grey atom in Supplementary Fig.~\ref{fig:dis_sche}c) is fixed along the dislocation line direction while all other degrees of freedom are fully optimized.  Under this constraint, the middle replica will evolve to the saddle-core configuration (peak energy in the transition path) found in the NEB method. Supplementary Table~\ref{tab:drag_vs_neb} shows the Peierls barriers computed by the NEB and the drag methods, respectively.  Their agreement shows that the drag method, as an alternative to NEB, can be used to compute the Peierls barrier in a reliable and computationally efficient manner.  The calculations for alloys are performed with the VCA method in VASP and convergence is assumed when the ionic forces fall below \(0.02\) eV/\AA .  The obtained Peierls barriers of W-TM alloys are shown in Fig.~\ref{fig:chi_prop}c.

\subsection*{Unstable stacking fault energy of binary W-TM alloys with DFT VCA in VASP}
\label{sec:usf_wre_vca_vasp}
The unstable stacking fault energy \(\gamma_\text{us}\) dictates the energy barrier to nucleate a dislocation.  In the current work, we calculate \(\gamma_\text{us}\) in the \(\langle 111 \rangle\) direction on the \(\{110\}\) plane as a function of solute concentration in binary BCC W-TM alloys by DFT VCA in VASP. In particular, \(\gamma_\text{us}\) of \{110\} and \{112\} planes are calculated for W-Re alloys and are used as material property inputs for analysis of crack tip behaviour using linear elastic fracture mechanics (LEFM). The DFT settings, supercell and convergence criterion are described in the Methods section: density-functional theory calculations using VASP, Generalized stacking fault energy surface and Supplementary Table~\ref{tab:supercell_dft}. The \(\gamma_{\text{us}}\) of \{112\} plane converges when the slab has 16 atomic layers or more. The converged \(\gamma_\text{us}\) values are shown in Fig.~\ref{fig:chi_prop}b and Supplementary Fig.~\ref{fig:usf_sf_wre_calc}a.

\subsection*{Surface energy of W-Re alloys with DFT VCA in VASP and QE}
\label{sec:wre_surf_e_vca_vasp}
The surface energies of binary W-Re alloys are calculated using VCA as implemented in both VASP and QE (see Supplementary Table~\ref{tab:supercell_dft} for supercells employed).  The pseudopotentials and parameters are described in the Methods section: Density-functional theory calculations using VASP and Cross-validation of \(\chi\) in binary W-TM alloys with DFT VCA in QE.  In addition, the kinetic cutoff energies for wavefunctions and charge densities are chosen as 952 eV (70 Ry) and 11424 eV (840 Ry) respectively in QE. For the ionic minimization, the convergence thresholds for total energy and forces are 1.36 $\times 10^{-3}$ eV (1 $\times 10^{-4}$ Ry) and 2.57 $\times 10^{-2}$ eV/\AA\ (1 $\times 10^{-3}$ Ry/bohr). The convergence criterion for electronic self-consistency is 1.36 $\times 10^{-5}$ eV (1 $\times 10^{-6}$ Ry).

For each surface, a fully periodic supercell is created with supercell vectors \((\mathbf{c}_1, \mathbf{c}_2, \mathbf{c}_3)\), where \(\mathbf{c}_1\) and \(\mathbf{c}_2\) are the in-plane supercell repeating vectors and \(\mathbf{c}_3\) is the out-of-plane vector.  We first obtain the energy per atom \(E_\text{c}\) in the bulk by optimizing the supercell and ionic positions.  The slab-vacuum supercell is then created by adding a vacuum layer in the out-of-plane direction. The slab-vacuum structure is optimized through ionic relaxation only in the out-of-plane direction while keeping all supercell vectors fixed.  In both calculations of the bulk supercell and slab-vacuum supercell, a consistent k-point sampling is used in the in-plane directions, leading to better convergence with respect to \(\mathbf{c}_3\), i.e., slab layers~\cite{sun_2013_ss}.  The surface energy \(\gamma\) is calculated as
\begin{align}
  \gamma = \dfrac{E_\text{sv} - n_\text{sv}E_\text{c}}{2A},
\end{align}
where \(E_\text{sv}\) and \(n_\text{sv}\) are the total energy and number of atoms of the slab-vacuum system and \(A\) is the area of the surface, i.e., \(A = \lvert  \mathbf{c}_1 \times \mathbf{c}_2 \rvert\).  The supercell vectors are shown in Supplementary Table~\ref{tab:supercell_dft} and a slab with 12 atomic layers gives the converged surface energies, as shown in Supplementary Fig.~\ref{fig:usf_sf_wre_calc}c,d.

For pure W in stardard DFT calculations, both VASP and QE give similar surface energies; \(\gamma_{\text{s} \{100\}}\) and \(\gamma_{\text{s}\{110\}}\) are 4.03 J/m\(^{2}\) and 3.27 J/m\(^{2}\) in VASP, and 4.15 J/m\(^{2}\) and 3.28 J/m\(^{2}\) in QE respectively.  In VCA QE, the surface energies of the \(\{100\}\) and \(\{110\}\) planes are predicted to gradually decrease with Re additions, consistent with previous calculations on W-Re alloys~\cite{yang_2018_prb,hu_2021_am}. However, in VCA VASP, the surface energies are unusual; they increase substantially with Re additions. For example, the surface energy is 11.12 J/m\(^{2}\) for the \(\{110\}\) plane of W-50\%Re, nearly 3.4 times of pure W (3.27 J/m\(^{2}\)). The VCA VASP-based results are also substantially different from that based on special quasi-random structure (SQS) supercell in standard VASP calculations~\cite{hu_2021_am}. In contrast, the results of VCA QE are close to those in SQS supercell calculations.  Therefore, we use the surface energies based on VCA QE as inputs in subsequent LEFM analysis.

\subsection*{Elastic constants of W-Re alloys with DFT VCA in VASP}
\label{sec:wre_elas_vca_vasp}
The elastic constants of binary W-Re alloys are calculated using VCA method with VASP.  The elastic contants are determined by fitting to the total energies versus strains within \(\pm\)1\%.  Supplementary Fig.~\ref{fig:usf_sf_wre_calc}b shows the elastic constants as a function of Re concentrations. Overall, Re has weak effects on the elastic constants of W.

\subsection*{Critical stress intensity factors for cleavage and dislocation emission}
\label{sec:lefm_wre}
Pure W exhibits brittle cleavage behaviour on \(\{100\}\) and \(\{110\}\) planes at low temperature~\cite{gumbsch_1998_science}.  The critical stress intensity factor \(K_\text{Ic}\) for Griffith cleavage is lower than \(K_\text{Ie}\) for dislocation emission~\cite{rice_1992_jmps} for sharp cracks on these two planes under mode-I loading. In binary W-Re alloys, \(K_\text{Ic}\) and \(K_\text{Ie}\) may change due to Re-effect. Here, we compute the respective \(K_\text{Ic}\) and \(K_\text{Ie}\) as a function of Re concentration for mode-I loaded sharp cracks.

Since W is nearly isotropic, we use isotropic LEFM theory. In particular, \(K_\text{Ic}\) for Griffith cleavage under plane strain Mode-I loading is computed as
\begin{equation}
  K_\text{Ic}=\sqrt{\frac{G}{D}},
\end{equation}
where
\begin{equation}
  G = 2 \gamma_\text{s}, \gamma_\text{s} \text{ is the surface energy}
\end{equation}
and
\begin{equation}
  D=\frac{1-{\nu}^{2}}{E}, \nu \text{ is the Poisson's ratio and } E \text{ is the Young's modulus}.
\end{equation}

The critical stress intensity for dislocation emission \(K_\text{Ie}\) is calculated according to the Rice criterion~\cite{rice_1992_jmps}
\begin{equation}
K_\text{Ie}=\frac{1}{f(\theta)} \sqrt{\frac{2\mu}{1-\nu} [1+(1-\nu)\tan^{2}\phi]\gamma_\text{us}}
\end{equation}
where \(f(\theta)=\cos^{2}(\theta/2)\sin(\theta/2)\), \(\theta\) is the angle between the crack plane and slip plane, \(\phi\) is the angle between the dislocation Burgers vector and crack front direction in the slip plane, and \(\mu\) is shear modulus, \(\nu\) is Possion's ratio, \(\gamma_\text{us}\) is the unstable stacking fault energy of the slip plane.

All material properties (\(\gamma_\text{s}, \gamma_\text{us}, E, \nu, \mu,\)) depend on Re concentrations in W and are calculated based on first-principle DFT VCA calculations. Supplementary Fig.~\ref{fig:kie_kic_wre_calc} shows the results for 4 crack systems (\(\{100\}\langle 001 \rangle\), \(\{100\}\langle 011 \rangle\), \(\{110\}\langle 001 \rangle\), \(\{110\}\langle 011 \rangle\)) susceptible to cleavage fracture in BCC W.  In pure W and for a \(\{100\}\) sharp crack (Supplementary Fig.~\ref{fig:kie_kic_wre_calc} c), \(K_\text{Ic} = 1.90 \text{ MPa}\cdot\text{m}^{1/2}\), which is \(\sim\)40\% and \(\sim\)42\% lower than \(K_\text{Ie}\) for dislocation emission in the \(\{100\}\langle 001 \rangle\) and \(\{100\}\langle 011 \rangle\) crack systems. Similar behaviours are seen in cracks on \(\{110\}\) planes (Supplementary Fig.~\ref{fig:kie_kic_wre_calc}g). In pure W, \(K_{\text{Ic}} = 1.72 \text{ Mpa}\cdot\text{m}^{1/2}\), which is \(\sim\)41\% and \(\sim\)31\% lower than \(K_\text{Ie}\) for dislocation emission in the \(\{110\}\langle 001 \rangle\) and \(\{110\}\langle 011 \rangle\) crack systems. The isotropic solution is nearly identical to that of the anisotropic solution~\cite{mak_2021_jmps} and both results suggest that \(\{100\}\) and \(\{110\}\) planes are intrinsically brittle in W.  With increasing Re concentrations, all \(K_{\text{Ie}}\) decrease sharply. However, \(K_{\text{Ic}}\) is insensitive to Re additions.  As a result, the gaps between \(K_{\text{Ie}}\) and \(K_{\text{Ic}}\) decrease sharply in all the four crack systems, which is potentially beneficial to make W ductile.  Nevertheless, no cross-over between \(K_{\text{Ie}}\) and \(K_{\text{Ic}}\) is predicted in the analysis based on material properties at 0 K.

To extend the prediction to room temperatures, we employ a finite-temperature ductility criterion \(D=K_{\text{Ie}}/K_{\text{Ic}}\) proposed recently~\cite{mak_2021_jmps} based on experimental observations. In BCC W, ductile behaviour can be achieved at room temperatures if \(D < D_\text{c}\); \(D_\text{c} \approx 1.56/1.26\) for \(\{100\}/\{110\}\) planes respectively. In W-20\%Re, \(D_{\{100\}\langle 001 \rangle} = 1.53\) and \(D_{\{110\}\langle 011 \rangle} = 1.30\), close to the respective threshold values \(D_\text{c}\) (Supplementary Fig.~\ref{fig:kie_kic_wre_calc}d,h).  With higher Re concentrations, \(D\) decreases further and satisfies the ductile criterion of both planes. At 20\% Re concentration (\(\chi=0.73\)), the ND core also transforms to the D core.  This concentration is also close to many commercial W-Re alloys exhibiting high ductility and a DBT below room temperatures.  The finite-temperature \(D_\text{c}\) is established empirically based on DBT experiments, it may include the effects of core transformation and associated change in plastic slip behaviour.  Since \(D > 1\) at 0 K, the mechanism of Re-ductilization effects is thus more likely related to core transformation.  However, we can not rule out thermal activation in dislocation nucleation at crack tips.  Nevertheless,  the above analysis is based on well-established linear elastic fracture mechanics and DFT-computed material properties. The primary assumption is linear elasticity.  Therefore, the results are robust and consistent with ductile behaviour in BCC W-Re alloys at room temperatures.  More importantly, the critical material properties (\(\gamma_\text{us}\) and ND/D core structures) and associated brittle vs ductile behaviour can be related to the material index \(\chi\).

\def\section*{}{\section*{}}
\section*{References}
\textbf{References 1-70 are for the main text.}
\vspace{-0.5cm}

\section{0pt}{4pt plus 4pt minus 2pt}{4pt plus 2pt minus 2pt}
\titlespacing\subsection{0pt}{4pt plus 4pt minus 2pt}{1pt plus 2pt minus 2pt}
\titlespacing\subsubsection{0pt}{4pt plus 4pt minus 2pt}{4pt plus 2pt minus 2pt}
\titlespacing\paragraph{0pt}{1pt plus 4pt minus 2pt}{2pt plus 2pt minus 2pt}

\setcitestyle{super}

\usepackage{graphicx}
\usepackage{natbib}
\usepackage{dcolumn}
\usepackage{amsmath, amsbsy, amssymb}
\usepackage{bm}
\usepackage{float,array}
\usepackage{textcomp}
\usepackage[version=3]{mhchem}
\usepackage{times}
\usepackage{helvet}
\usepackage{multirow}
\usepackage{booktabs}
\usepackage{mathtools}
\usepackage[]{hyperref}
\hypersetup{colorlinks,pdfborder={0,0,0},linkbordercolor={white},citecolor=MidnightBlue,filecolor=MidnightBlue,linkcolor=MidnightBlue,urlcolor=MidnightBlue}
\usepackage{pbox}
\usepackage{diagbox}
\usepackage{makecell}
\usepackage{tikz}
\usepackage{array}
\usepackage{booktabs,adjustbox}
\usepackage[UKenglish,USenglish,english,american]{babel}
\usepackage[]{upgreek}

\newcolumntype{L}[1]{>{\raggedright\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\newcolumntype{C}[1]{>{\centering\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\newcolumntype{R}[1]{>{\raggedleft\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}

\usepackage{forest}
\usetikzlibrary{arrows,arrows.meta,decorations.markings,shapes.arrows,shapes.geometric,calc,shadows}

\usepackage{verbatim
\makeatletter
\newcommand{\nextverbatimspread}[2]{%
  \def\verbatim@font{%
    \linespread{#1}\normalfont\ttfamily#
    \gdef\verbatim@font{\normalfont\ttfamily#2}
}
\newenvironment{myverbatim}[1][1.5]
  {\def\verbatim@font{%
     \linespread{#1}\normalfont\ttfamil
     \gdef\verbatim@font{\normalfont\ttfamily}
   \verbatim}
  {\endverbatim}
\renewcommand{\p@subsection}{}
\renewcommand{\p@subsubsection}{}
\makeatother

\usefont{T1}{cmr}{m}{n}
\renewcommand{\figurename}{{Figure}}

\newcommand*\mycommand[1]{\texttt{\emph{#1}}}

\newcommand{$\langle \bold{a} \rangle $ }{$\langle \bold{a} \rangle $ }
\newcommand{$\langle \bold{c} \rangle $ }{$\langle \bold{c} \rangle $ }
\newcommand{$\langle \bold{c}+\bold{a} \rangle $ }{$\langle \bold{c}+\bold{a} \rangle $ }
\newcommand{\langle \bold{a} \rangle  }{\langle \bold{a} \rangle  }
\newcommand{\langle \bold{c} \rangle  }{\langle \bold{c} \rangle  }
\newcommand{\langle \bold{p} \rangle  }{\langle \bold{p} \rangle  }
\newcommand{\langle \bold{q} \rangle  }{\langle \bold{q} \rangle  }
\newcommand{\langle \bold{c}+\bold{a} \rangle  }{\langle \bold{c}+\bold{a} \rangle  }

\newcommand{$\langle 111 \rangle /2 $ }{$\langle 111 \rangle /2 $ }

\newcommand{\mbox{\normalfont\AA}}{\mbox{\normalfont\AA}}
\newcommand{\mathbin{\!/\mkern-5mu/\!}}{\mathbin{\!/\mkern-5mu/\!}}

\newcommand{\Vect}[1]{\mathbf{#1}}
\newcommand{\Tens}[1]{\mathbf{#1}}

\def\Large Table of Contents{\Large Table of Contents}
\usepackage{nameref}

\hbadness=10000
\hfuzz=50pt