# A comprehensive screening of plasma-facing materials for nuclear fusion

Andrea Fedrigucci1,2<EMAIL_ADDRESS>Nicola Marzari1,3 Paolo Ricci2
1Theory and Simulation of Materials (THEOS), and National Centre for
Computational Design and Discovery of Novel Materials (MARVEL), École
Polytechnique Fédérale de Lausanne (EPFL), Lausanne CH-1015, Switzerland
2Swiss Plasma Center (SPC), École Polytechnique Fédérale de Lausanne (EPFL),
Lausanne CH-1015, Switzerland
3Laboratory for Materials Simulations, Paul Scherrer Institut, Villigen PSI
5232, Switzerland

###### Abstract

Plasma-facing materials (PFMs) represent one of the most significant
challenges for the design of future nuclear fusion reactors. Inside the
reactor, the divertor will experience the harshest material environment:
intense bombardment of neutrons and plasma particles coupled with extremely
large and intermittent heat fluxes. The material designated to cover this role
in ITER is tungsten (W). While no other materials have shown the potential to
match the properties of W, many drawbacks associated with its application
remain, including: cracking and erosion induced by a low recrystallization
temperature combined with a high ductile–brittle transition temperature and
neutron-initiated embrittlement; surface morphology changes (He bubbles and
fuzz layer) due to plasma-W interaction with subsequent risk of spontaneous
material melting and delamination; low oxidation resistance in case of air
contamination. Exploring alternatives to W requires the design of a
multivariable optimization problem. This work aims to produce a structured and
comprehensive materials screening of PFMs candidates based on known inorganic
materials. The methodology applied in this study to identify the most
promising PFM candidates combines peer-reviewed data present in the Pauling
File database and first-principle DFT calculations of two key PFMs defects,
namely the surface binding energy and the formation energy of a hydrogen
interstitial. The crystal structures and their related properties, extracted
from the Pauling File, are ranked according to the heat-balance equation of a
PFM subject to the heat loads in the divertor region of an ITER-like tokamak.
The materials satisfying the requirements are critically compared with the
state-of-the-art literature, defining an optimal subset where to perform the
first-principles electronic structure calculations. The majority of previously
known and extensively studied PFMs, such as W, Mo, and carbon-based materials,
are captured by this screening process, confirming its reliability.
Additionally, less familiar refractory materials suggest performance that
calls for further investigations.

††preprint: APS/123-QED

##  I Introduction

One of the major challenges on the way to fusion energy - one of the few
options we have to provide mankind with a large-scale, sustainable and carbon-
free energy source [1, 2] \- is linked to the availability of in-vessel
plasma-facing materials (PFMs) capable to sustain the severe operational
condition of a fusion power plant [3] for an assumed lifetime of 1-2 full
power years (fpy) [4, 5]. In particular, in modern tokamaks, such as ITER [6],
the most external magnetic field lines are deviated toward the divertor, and
the thermal loads that the divertor must withstand are the highest of all the
components facing the plasma, making it one of the most challenging parts to
build [7]. The divertor is expected to exhaust the heat and control the level
of helium ash and other impurities, mainly resulting from erosion of the
plasma-exposed surface, to minimize plasma contamination. In addition, the
divertor protects the most delicate surrounding components from neutronic
loads. The heat load on the divertor primarily results from a combination of
five types of loads: ionic, neutral, radiation power, electrons, and surface
recombination energy of incoming ion-electron pairs [8]. The power deposited
on the divertor surface is dominated by radiative and particle fluxes, with
the latter mainly constituted by helium (He), hydrogen isotopes (deuterium and
tritium), and impurities. Their flux at the divertor is the highest of the
vacuum vessel, ${10}^{24}$ particles
${\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}$ with an average kinetic energy
of 10 $\mathrm{eV}$ (during steady-state) [9, 10]. Last, 14.1 $\mathrm{MeV}$
neutrons bombard the PFM with a flux of ${10}^{18}$ particles
${\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}$ [9].

In ITER the thermal load at the divertor in quasi-stationary conditions is
estimated to oscillate between 10 $\mathrm{MW}\text{\,}{\mathrm{m}}^{-2}$, in
steady-state conditions lasting around 400 $\mathrm{s}$, and 20
$\mathrm{MW}\text{\,}{\mathrm{m}}^{-2}$, during the slow-transient thermal
loads, where a transition from detached divertor operation to fully attached
operation is observed ($<$ 10 $\mathrm{s}$) [11, 9]. Moreover, intense and
short transient loads are predicted to occur due to edge-localized modes
(ELMs), that are quasi-periodic fluctuations of the plasma pressure at the
plasma edge [9]. The thermal energy released on the divertor target during
unmitigated ELMs can be as high as 10-15
$\mathrm{MJ}\text{\,}{\mathrm{m}}^{-2}$ with a frequency of 1-2 $\mathrm{Hz}$
and a duration of 0.1-1 $\mathrm{ms}$ [12, 13, 14, 15, 16]. In fact, ELM
mitigation is required in ITER, and it will be based on increasing ELM
frequency to 25-50 $\mathrm{Hz}$, consequently lowering their heat flux to
around 0.5 $\mathrm{MJ}\text{\,}{\mathrm{m}}^{-2}$ and reducing the heat peak
from tens to tenths of $\mathrm{GW}\text{\,}{\mathrm{m}}^{-2}$ [12, 17, 18].
Finally, off-normal events, which may not be excluded during ITER operations,
include vertical displacement events (VDEs), which result in energy pulses up
to 30 $\mathrm{MJ}\text{\,}{\mathrm{m}}^{-2}$ in 1-5 $\mathrm{ms}$, plasma
disruption, which can reach around 60 $\mathrm{MJ}\text{\,}{\mathrm{m}}^{-2}$
and last 100-300 $\mathrm{ms}$, and runaway electrons, which lead to 50
$\mathrm{MJ}\text{\,}{\mathrm{m}}^{-2}$ pulses for less than 50 $\mathrm{ms}$
[9, 15, 19]. The simultaneous thermal, plasma, and neutron wall loads lead to
complex dynamics and degradation process of the divertor armor. The initial
ITER divertor design included both carbon fiber composite (CFC) and tungsten
(W) as PFMs, but given the high-level of tritium retention and surface
erosion, the ITER Organization (IO) decided to switch to a fully W divertor
since 2013 [20, 21]. While W has the highest melting point of all the metals
and outstanding thermal properties, it suffers of (i) low recrystallization
temperature combined with high ductile–brittle transition temperature (DBTT);
(ii) neutron-induced embrittlement that leads to cracking and erosion of the W
surface at high thermal loads; (iii) surface morphology changes (He bubbles
and fuzz layer) due to plasma-W interaction with subsequent risk of
spontaneous material melting and delamination; (iv) neutron-induced damages
and transmutation affecting thermal and mechanical properties; and (v) low
oxidation resistance in case of air contamination [7, 22]. While the PFM
community is working to propose solutions for some of these issues, it is
still unclear if a W-based divertor will be able to sustain operating fusion
conditions in a fusion power plant. For instance, the design of the divertor
is still open for DEMO [23], the power plant that will demonstrate the
production of electric power from nuclear fusion. The primary PFM of DEMO
divertor will withstand the same thermal loads as the one in ITER, but higher
neutron loads (from 0.37–0.47 $\mathrm{MW}\text{\,}{\mathrm{m}}^{-2}$ to
1.8–2.4 $\mathrm{MW}\text{\,}{\mathrm{m}}^{-2}$) [24, 25]. The same risks
previously mentioned related to erosion, melting, and cracking of the W
surface will be faced. The more intense neutron bombardment will also lead to
stronger displacement damage and transmutation effects, resulting in swelling,
irradiation hardening, and radioactivity after exposure. These challenges pose
even greater difficulties in the material choice for DEMO and future
commercial fusion reactors.

Given the high cost and time required for both testing or simulating PFMs in
such a complex environment [26, 27, 28], only few classes of materials have
been evaluated as potential alternatives to W [29] and a high-throughput
protocol for PFMs screening is missing in the literature. This multivariable
optimization problem requires to define a set of design criteria that a
material need to satisfy in order to sustain divertor conditions. In this
study we define the first methodology to assess the potential of a large class
of materials as PFM. The heat loads that PFMs must endure are modeled using a
heat balance equation, which takes into account both steady-state and
transient loads. The thermal properties involved in this equation, such as
thermal conductivity, melting temperature, heat capacity, and other related
proxies are extracted from the inorganic crystal database MPDS [30], an online
platform providing access to the PAULING FILE [31], which contains the peer-
reviewed properties of, approximately, half a million experimentally-confirmed
structures. The best performing materials are ranked and then selected to
calculate, in the framework of density-functional theory calculations (DFT)
[32], surface binding energies (SBE), which are a proxy of the required energy
for physical sputtering and therefore the erosion rate, and the hydrogen
interstitial formation energies (H-IFE), as an indication of the hydrogen
solubility into the crystal [33] and therefore its tritium retention.

The impact of neutron bombardment is not explored here, due to the significant
computational burden associated with modeling its effects on candidate
materials. Studies in the literature, such as those in Ref. [34, 35], have
calculated the radioactivity of pure chemical elements following neutron
bombardment in DEMO-like reactors. This information could serve as an
additional design criterion, as threshold values have been proposed [36].

The paper is organized as follows. In Sec. II, we discuss the design criteria,
data processing of the data extracted from MPDS, and the ranking procedures
used to select the PMFs. In Sec. III, we show the best candidates emerging
from the ranking criteria and in Sec. IV we present a comprehensive literature
review on their applicability as PFMs. Only the materials that pass both
selection procedures undergo ab-initio calculations of SBE and H-IFE, and
their results are presented in Sec. V. In Sec. VI we offer a conclusive
evaluation of the top PFM candidates. Finally, Sec. VII provides a perspective
on the development of the field of PFMs.

##  II Methods

The set of design criteria needed for finding an exhaustive list of PFMs
candidates should include the effects produced by the three loads on the
exposed material: heat, neutrons, and plasma particles. In this work, a first
material screening, we consider known or easy-to-compute properties of
experimentally observed crystal, therefore excluding the analysis of complex
dynamics emerging at the atomic-scale level from degradation processes. The
design criteria include the constraint dictated by thermal loads (Sec. II.1),
and the plasma-wall interactions (Sec. II.2). For the heat load, the heat
balance equation of a solid material exposed to a base load of steady-state
heat flux and short-transient loads (ELMs) is being investigated as the
primary design criterion. In Sec. II.1.1, the maximum thickness that a given
material possesses to avoid surface’s melting is used as first ranking
procedure of the material retrieved from MPDS. If this cannot be estimated due
to a lack of available properties, the concepts of Pareto fronts and win-
fraction are defined in Sec. II.1.2 to determine how materials performs in
terms of thermal properties. To compare the consistency of the two ranking
procedures the idea of a comparative ranking is presented in Sec. II.1.3.
Plasma-wall interactions are then explored by calculating the formation energy
of two common PFM defects, the surface binding energy (SBE) and the hydrogen
interstitial formation energy (H-IFE). The relations between these crystal
defects, the sputtering phenomena, and the tritium retention process are
explained in Sec. II.2.1 and II.2.2, together with the Python workflows built
to automate the first-principle calculation of these properties starting from
a crystal structure.

###  II.1 Thermal loads

All materials-related data of this section are retrieved from the inorganic
crystal database MPDS [30]. In MPDS multiple values can be associated to a
single material property. Indeed, the database is composed of entries, and
each entry corresponds to the value of a peer-reviewed measured (or,
occasionally, calculated) property of a specified crystal. Multiple entries
for the same property and material can differ because of their literature
reference and/or experimental or simulated conditions, which are not always
reported. For this reason, in this work, the value of a property of a material
is defined through the use of the kernel density estimation method (KDE) [37].
This choice is based on the empirical observation that most of the reported
values of a given property are around standard conditions ($T$ = 298
$\mathrm{K}$, $P$ = 1 atm). The KDE creates a probability density function
estimation, $f$, from the finite data sample of the available entries,
$x_{i}$, that is:

$f(h,x)=\sum_{i=1}^{n}K\left(x-x_{i};h\right)$ (1)

where $h$ is the bandwidth parameters setting the level of smoothness of
$f(h,x)$ ($h$ = 0.05 in this work) and $K$ is a kernel function (here a
Gaussian $K=exp[-9(x-x_{i})^{2}/(2h^{2})]$). The value of a property of a
crystal used in this study is thus obtained as the value that maximises
$f(h,x)$, unless otherwise specified. We note that some entries are calculated
under high-pressure conditions; sometimes, MPDS identifies such entries with
the ”hp” flag in the crystal structure chemical formula. All entries with the
”hp” flag are excluded from the screening process. In addition, 17 more
materials, selected based on ranking criteria, are excluded from the dataset
because their reported properties are measured under high-pressure conditions,
as indicated in the literature references provided within the database.

The underlying assumption of the design criteria defined in this section is
that a candidate material cannot melt during the steady-state (SS) tokamak
operating conditions and ELM events, in order to keep its shielding functions
and to avoid plasma contamination. In SS conditions the heat flux $\Phi_{q}$
to the PFM equals the heat flux that is removed on the cooling side, by
assuming one dimensional thermal conduction through a material of thickness
$d$ [38]:

$\Phi_{q}^{SS}=\frac{1}{d}\int_{T_{bott}}^{T^{SS}_{surf}}k(T)\mathrm{d}T$ (2)

where $T^{SS}_{surf}$ and $T_{bott}$ are the temperature on the plasma-facing
and the cooling sides. Assuming constant thermal conductivity $k(T)$ =
$\langle k(T)\rangle$, Eq. (2) leads to:

$T^{SS}_{surf}=T_{bott}+\Phi_{q}^{SS}\frac{d}{\langle k(T)\rangle}$ (3)

In contrast, during short thermal transients, such as ELMs, the heat is
deposited on the surface as a power pulse. The temperature at the end of the
ELM event ($T^{ELMs}_{surf}$) can be estimated as:

$T^{ELMs}_{surf}=T^{SS}_{surf}+\Phi_{q}^{ELMs}\sqrt{t}\frac{2}{\sqrt{\pi\rho\langle
C_{p}(T)\rangle\langle k(T)\rangle}}$ (4)

if the penetration depth $d^{\ast}$ of the heat in the material surface is
smaller than the material thickness $d$ [39]:

$d\geq d^{\ast}=2\sqrt{t}\sqrt{\frac{\langle k(T)\rangle}{\rho\langle
C_{p}(T)\rangle}}$ (5)

where $t$ is the duration of an ELM event, $\Phi_{q}^{ELMs}$ is the power flux
density during ELMs, $\langle C_{p}(T)\rangle$ is average heat capacity of the
material and $\rho$ its density. Eq. (4) can be used to define a filter to
assess the potential of a material to withstand the divertor thermal loads
without melting ($T_{m}>T^{ELMs}_{surf}$), having excluded off-normal events.
The material-dependent properties appearing in Eq. (4) are retrieved from the
MPDS database. The thickness $d$ is a target parameter, that is partially
constrained by the erosion rate at the divertor, its assumed lifetime (e.g.
twice as much CFC is needed compared to W metal [40] for a comparable lifetime
of the armor material), and the material response to thermal stress.

Two ranking procedures are then pursued to select the most promising PFM
candidates. The first is the direct application of the heat balance equations
and is called thickness ranking; the second builds on the various material-
related properties extrapolated from the same equations and selects materials
based on a Pareto optimization procedure refined by a win-fraction ranking
equation, and is called Pareto ranking. A third ranking procedure is applied
with the main goal of comparing the accuracy of the two ranking methods and
establish a unique ranking for all the materials retrieved from MPDS, which we
denote as comparative ranking.

####  II.1.1 Thickness ranking

The first ranking method comes directly from the application of the heat
balance equation in Sec. II.1. The list of materials for which is possible to
collect all the variables in Eqs. (4) and (5) are ranked accordingly to the
maximum thickness allowed to avoid melting during type-I ELMs ($d^{ELMs}$). If
$C_{p}(T)$ is not available, then the maximum thickness allowed to avoid
melting during steady-state ($d^{SS}$) is used as ranking criterion. The input
values assumed to compute the material thickness are reported in Tab. 1.

Table 1: Input parameters used to compute the heat balance equations: Eq. (3), (4) and (5). Parameter | Unit | Value | Ref.
---|---|---|---
$T_{bott}$ | $\mathrm{K}$ | 415 | [41]
$Q^{SS}$ | $\mathrm{MW}\text{\,}{\mathrm{m}}^{-2}$ | 20 | [9, 17]
$t$ | $\mathrm{ms}$ | 0.5 | [12, 17]
$Q^{ELMs}$ | $\mathrm{GW}\text{\,}{\mathrm{m}}^{-2}$ | 1 | [12, 17]

Given the large uncertainty around the real value of each property
extrapolated from MPDS, the thicknesses calculated from Eqs. (3) and (4) are
reported with a standard deviation (SD) that is obtained by the propagation of
the standard deviation of every single property. Materials with $d^{ELMs}>$ 0
$\mathrm{mm}$ and $d^{ELMs}-d^{\ast}>$ 0 $\mathrm{mm}$ and materials with
$d^{SS}\geq 3$ $\mathrm{mm}$ are considered PFM candidates; materials not
passing this test are discarded. The threshold value of 3 $\mathrm{mm}$ for
steady-state conditions is chosen since no material with a positive $d^{ELMs}$
has a lower $d^{SS}$ thickness.

It is important to remember that these values are calculated using three
assumptions: (i) the value at the maximum of the KDE probability density
function captures well the real value of the material property, (ii) these
properties remain constant over large temperature ranges, (iii) degradation
processes and materials defects are neglected. Therefore, the calculated
$d^{ELMs}$ and $d^{SS}$ remain indicative of a trend and should not be
considered for the engineering design.

####  II.1.2 Pareto ranking

For materials where neither $d^{ELMs}$ nor $d^{SS}$ can be calculated a
ranking process that does not require the knowledge of an objective function
is used. This is performed in order to extract materials from MPDS that are
showing promising performance in a subset of the properties identified through
thermal balance analysis. In order to expand the number of available PFM
candidates, the properties retrievable from MPDS that showed a high Pearson
correlation coefficient ($\geq$ 0.7) with $T_{m}$, $k$, $C_{p}$ are added to
the list of properties. Of all the properties available on MPDS, $T_{m}$ shows
a Pearson correlation coefficient greater or equal than 0.7 along with 6 other
properties: shear modulus ($G$), Young’s modulus ($E$), adiabatic and
isothermal bulk modulus ($B_{S}$ and $B_{T}$), Mohs hardness ($MH$), Knoop
hardness ($KH$); $C_{p}$ with the heat capacity at constant volume ($C_{v}$).
At this stage, the total number of properties is 10 and they are categorized
as follows: (i) melting temperature + proxies, (ii) thermal conductivity,
(iii) heat capacity (at constant pressure and constant volume). In the
following, for the sake of clarity, these 3 macro-categories are identified
with the term categories ($c$), whereas the 10 features are called
interchangeably features or properties ($p$) (see Tab. 2).

Table 2: The 10 properties ($p$) retrieved from MPDS are classified into 3 categories ($c$). This classification is applied inside the Pareto and comparison ranking systems. Category ($c$) | Property ($p$)
---|---
Melting temperature + proxies | $T_{m}$
$G$
$E$
$B_{S}$
$B_{T}$
$MH$
$KH$
Thermal conductivity | $k$
Heat capacity | $C_{p}$
$C_{v}$

Note that the density is always available for every crystal structure reported
in MPDS, therefore instead of using the heat capacity in molar units, it is
converted in unit of volume as Eq. (4) suggests.

A rigorous way to treat an optimization problem for which the objective
function is unknown is through the so-called Pareto set (or Pareto front), a
collection of solutions (i.e. materials) that cannot be improved in one
objective (i.e. property) without worsening at least one of the others [42].
Given a set of materials $m$ of size $M$ that share the same set of properties
$p$ of size $P$, the Pareto set $\mathcal{P}$ of this system is defined as
follow:

$\displaystyle\begin{aligned} &m_{i}\in\mathcal{P}\Leftrightarrow\nexists
m_{j}:\left\\{\begin{array}[]{l}p_{j,k}\geq p_{i,k}\quad\forall
k=1,\ldots,P\\\ \exists k_{0}:p_{j,k_{0}}>p_{i,k_{0}}\end{array}\right.\\\
&\forall i,j=1,\ldots,M\end{aligned}$ (6)

First, upstream of the Pareto optimization process, the data are standardized
by subtracting its mean, so that the mean of the standardized values vanishes,
and divided by the standard deviation, ensuring that the standard deviation of
the standardized values is one. Then, different sets of materials that posses
the same set of properties to apply the Pareto optimization are defined. Since
the whole set of 10 $p$ is not available for all the materials, it is crucial
to define their subsets for which a Pareto set optimization is pursued. Here,
we apply Pareto for every material possessing at least 1 property out of the
10 previously defined. Therefore, three subsets of materials are defined as
follows: $3p$, which includes materials possessing at least one property for
each of the three category; $2p$, which encompasses materials having at least
one property in, at least, two different category; and $1p$, which includes
materials with only one property per category. Then, for every subset of
$3p$/$2p$/$1p$, all possible combinations of $n$ properties that belong to $n$
different categories are used to define the feature space. For instance, in
the $3p$ subset an allowed combination of features can be either
$T_{m}$-$C_{p}$-$k$ or $T_{m}$-$C_{v}$-$k$ but not $T_{m}$-$B_{S}$-$k$ or
$T_{m}$-$C_{v}$-$C_{p}$ since $T_{m}$-$B_{S}$ and $C_{v}$-$C_{p}$ are part of
the same category (respectively melting temperature + proxies and heat
capacity).

If the number of features is exhaustive for the description of a given
problem, any solution that is not a Pareto front can be discarded from the
list of interesting candidates, because it is worse in terms of performance.
In the context presented here, the list of properties used as design criteria
for an armor divertor material is far from being comprehensive; therefore,
multiple Pareto set are considered as optimal solutions. Once a Pareto set is
found, it is discarded from the list of candidates and those are re-evaluated.
Each Pareto set can be conceptualized as a Pareto layer ($l$), which, once
identified, is peeled off, and the subsequent layer is selected. This process
is repeated for a specified number of cycles.

To distinguish between any of the candidates thus selected, the concept of win
fraction, as proposed in Ref. [42], is revisited and applied.

$WF_{i}=\frac{\sum_{j}\sum_{k}\left|p_{i,k}-p_{j,k}\right|\times\mathcal{H}\left(p_{i,k}-p_{j,k}\right)}{\sum_{j}\sum_{k}\left|p_{i,k}-p_{j,k}\right|}$
(7)

where $WF_{i}$ is the win fraction of material $i$, $p_{i,k}$ is the value of
the property $k$ of material $i$ and $\mathcal{H}(x)$ is the Heaviside step
function that is equal to 0 for $x<0$ and to 1 for $x>0$. Contrary to the
original reference, the ranking equation is not defined as the smallest win
fraction that each materials $i$ posses when compared with the other
candidates $j$; instead, it corresponds to a summation extended over all
candidates of the selected Pareto fronts. In addition, it is implied that an
unitary weight is assigned to every property of the summation. Once the
ranking is pursued, the best candidates are selected. The number of Pareto
sets $l$ and the threshold value of $WF_{i}$ used to identify the best
candidates for each category of Pareto solutions is reported in Tab. 3.

Table 3: A Pareto ranking is utilized to identify the most promising candidates that do not meet the criteria for thickness ranking due to insufficient retrievable properties. The materials were categorized into three subsets based on the number of available properties. For each subset, the first $l$ Pareto fronts were selected, and only those with a $WF_{i}$ 7 value exceeding a specified threshold were considered as PFM candidates. Material class | $\bm{l}$ | $\bm{WF_{i}}$
---|---|---
3p |  $\leq$ 5 |  $\geq$ 50%
2p |  $\leq$ 5 |  $\geq$ 70%
1p |  $\leq$ 100 |  $\geq$ 90%

In order to avoid false positive, materials discarded from more stringent
criteria of selection should not be captured by weaker filters. Hence,
materials that are excluded by the thickness ranking scheme or by a Pareto
optimization with a higher number of properties are removed from the list of
materials subject to a less stringent procedure. Fig. 1 introduces a visual
representation that illustrates the materials space and the process of
exploring it.

Figure 1: Schematic representation of the material space in MPDS and its
exploration based on the number of available properties and the corresponding
ranking procedure, allowing for the identification of promising candidates
through the application of multiple ranking criteria.

####  II.1.3 Comparative ranking

The third ranking procedure is then created as a visualization tool to collect
all available materials in a single table, to rank them within a unique
framework, and as a benchmark to compare the results of the previous two
ranking scheme. For this procedure, the main issue is connected to the fact
that materials with different available properties must produce a comparable
ranking value. The win fraction concept might be re-used to define a ranking
value for every available material, after applying some changes to the
formulation presented in the previous section. Eq. (7) is intended to be
applied to the set of properties that a material possesses, and the absence of
certain features does not impact the final value of $WF_{i}$. If $WF_{i}$ is
used directly, materials with one single outperforming feature might rank at
the top of the list, obscuring the rest of the candidates. On the other hand,
if the win fraction is utilized for every material across the 10 retrieved
properties, even in the presence of missing data, it becomes essential to
specify how the missing values affect the win fraction. This depends on two
factors: the weight assigned to each property and the value to be assigned to
the missing ones. A unitary weight assigned to all the properties would reward
materials with a lot of available data at the expense of others that have a
smaller set of essential and high-performing features. A good compromise was
found by subdividing the properties in three different categories ($c$), as
proposed in Tab. 2. Then, the weights are built by multiplying two
coefficients:

  * •

The first term is a constant value that applies to all materials and is
defined as 1/$C$, where $C$ is the number of categories. In this work, we
always consider three categories, so the first term is equal to 1/3.

  * •

The second term is material-dependent and it is defined as 1/$P_{i,z}$, where
$P_{i,z}$ is the number of available properties of the $z$-th category of the
$i$-th material.

In this way, materials lacking in the key categories - melting temperature +
proxies, thermal conductivity, and heat capacity - are penalized, but the
number of properties retrievable within a category does not affect the weight.
Then, the weight of a given property is multiplied by the corresponding
single-property win fraction $WF_{i,k}$, and the summation over the all
features builds the weighted average win function $\overline{WF_{i}}$ as
presented in the following equation:

$\displaystyle\begin{aligned} \overline{WF_{i}}&=\sum_{z,k}w_{i,z}WF_{i,k}\\\
&=\sum_{z=1}^{3}\frac{1}{3}\sum_{k=1}^{P_{i,z}}\frac{1}{P_{i,z}}WF_{i,k}\\\
\end{aligned}$ (8)

where $k$ represent the $k$-th property of the $z$-th category of the $i$-th
material. If a category is not present for a given material, fictitious
$WF_{i,k}$ = 0 and $P_{i,z}$ = 1 are assigned. This methodology can be
generalized for any number of categories and any of their subsets. For
instance, the melting temperature and proxies category might be split between
actual $T_{m}$ and proxies, generating a subset with two sub-categories. The
results remain very similar whether or not this additional subcategory is
included. Therefore, for the sake of simplicity, this approach is not
discussed further.

###  II.2 Plasma-wall interactions

The interactions between hydrogen isotopes, alpha particles and the PFM have a
complex dynamics, which requires large-scale (105-106 atoms) and long-
timescale (102-103 $\mathrm{ns}$) molecular dynamics (MD) simulations [43,
44]. For instance, the entire formation process of a W fuzz layer under plasma
bombardment is out of the timescale explorable even with classical MD [45].
With the goal of finding indicators that reflect some of the undesired
phenomena PFMs may experience when interacting with a plasma particle flow,
while being computationally efficient, we use two key defects, as proxies,
surface binding energies (SBE) and hydrogen interstitial formation energies
(H-IFE), whose formation energies can be calculated within the DFT framework.

####  II.2.1 Surface binding energies

As mentioned above, the minimum thickness of a PFM is highly constrained by
the rate of erosion of the materials surface under the plasma impingement,
which in turn is measured through the sputtering flux of the surface atoms.
For physical sputtering, the energy distribution of the sputtered particle
flux, $\Phi(E_{e})$, can be described within the binary collision
approximation (BCA) by the Thompson energy spectrum [46]. The Thompson energy
spectrum highlights an important proportionality between $\Phi(E_{e})$ and the
surface binding energy (SBE or $E_{sb}$) of a material:

$\Phi(E_{e})\propto\frac{E_{e}}{(E_{e}+E_{sb})^{3}}$ (9)

This relation has a distribution that peaks at $E_{e}$=$E_{sb}$/2, where
$E_{e}$ is energy of the sputtered (emitted) particles leaving the surface. It
is notable that the tail of the Thompson energy spectrum deviates from
experimental values when the projectiles bombarding the surface have low-
energy ($<$ 500 $\mathrm{eV}$) and it is corrected by the Falcone energy
spectrum [47]. In both cases, SBEs highly affect the sputtered particle flux,
therefore the erosion rate and the minimum thickness required from a given
PFM. The SBE is defined as the energy that is necessary to remove an atom from
the top surface layer (in vacuum) during the physical sputtering process. We
calculate SBEs from first principles, and we use these as a proxy of the
material’s resistance to physical sputtering. SBEs are often estimated using
elemental enthalpies of sublimation in sputtering calculations. However,
recent literature emphasizes that first-principles SBEs are essential to
achieve more accurate outcomes [48, 49, 50]. In the context of density-
functional theory (DFT), these are calculated through the following relation:

$E_{sb}=E_{a}+E_{s+v}-E_{s}$ (10)

where $E_{a}$ is the DFT total energy of the single isolated (sputtered) atom
$a$, $E_{s}$ is the total energy of the crystal slab with a perfect surface,
and $E_{s+v}$ is the energy of the unrelaxed (for efficiency) slab with one
surface vacancy at the position that was atom $a$. $E_{sb}$ is obtained by
developing the Surface Binding Energy (SBE) workflow, a semi-automated Python
script created to make standardized SBE calculations (in-house, forthcoming).
Once the candidate is selected and its crystal structure is defined, 6 main
steps follow:

  1. 1.

Relaxation of atomic positions and cell shape of the primitive crystal
structure as obtained by MPDS, and calculation of its total energy.

  2. 2.

From the relaxed lattice constants, all the inequivalent slabs that can be
generated with Miller indices $hkl$ up to 1 are created and relaxed and their
surface energy is obtained through the following relation:

$E_{surf}^{hkl,\alpha}=\frac{E_{s}^{hkl,\alpha}-\frac{N_{s}^{hkl,\alpha}}{N_{b}}E_{b}}{2A_{s}^{hkl}}$
(11)

where $E$, $N$, $A$ are the total energy, the number of atoms, and the area of
the crystal slab, $s$ and $b$ are the subscripts indicating slab and bulk
structures, and $hkl$ and $\alpha$ indicate the Miller index and the
termination of a specific facet of the crystal structure. A minimum distance
of 10 $\mathrm{\SIUnitSymbolAngstrom}$ of vacuum between repeated slabs is
imposed as discussed in Refs. [51, 52].

  3. 3.

The slab with the lowest $E_{surf}^{hkl,\alpha}$ is selected together with all
possible terminations of the same Miller index, supercells are created and
their total energy ($E_{s}$) is obtained. A minimum distance of 9
$\mathrm{\SIUnitSymbolAngstrom}$ between periodic images is imposed to reduce
vacancy-vacancy interaction.

  4. 4.

One atom is removed from the surface of each supercell (per species and
inequivalent positions) and the total energy ($E_{s+v}$) is calculated. One
atom is considered a surface atom if it lies within 1
$\mathrm{\SIUnitSymbolAngstrom}$ from one of the two surfaces. Atomic
relaxation is neglected within this approximation.

  5. 5.

The total energy of a single atom ($E_{a}$) in a cubic cell of side
$\sqrt[3]{4000}$ $\mathrm{\SIUnitSymbolAngstrom}$ is calculated for each
element of the crystal.

  6. 6.

$E_{sb}$ is calculated following Eq. (10); if inequivalent sites of the same
chemical element are found on the surface of given terminations $\alpha$
and/or multiple termination $\alpha$ are available, $E_{sb}$ is calculated for
each element as the arithmetic mean between all inequivalent sites of each
termination.

All DFT simulations are performed by using the open-source computer code for
electronic-structure calculations Quantum ESPRESSO (QE) [53, 54] and all
parameters required to run the calculation are retrieved by Quantum Espresso
Input Generator website [55, 56, 57] set for spin-unpolarized calculations,
with maximum k-points distance of 0.2 ${\mathrm{\SIUnitSymbolAngstrom}}^{-1}$,
a smearing width for metals of 0.2 $\mathrm{eV}$ and using the
pseudopotentials library SSSP Efficiency PBEsol (version 1.2.0), where
exchange and correlation effects are approximated using the Perdew-Burke-
Ernzerhof revised for solids (PBEsol) functional. The kinetic energy cutoff
for both charge density and wavefunctions is consistent across all
calculations for a specific crystal and is set as the highest value suggested
by the Quantum Espresso Input Generator among the elements present in the
crystal.

Once $E_{sb}$ is known, the energy that a projectile must posses to sputter
one atom from a perfect-crystal surface in a head-on elastic collision has to
be greater or equal to the kinetic energy transferred to the recoil atom
$E_{T}$:

$E_{p}\geq E_{T}=E_{sb}\frac{\left(m_{p}+m_{a}\right)^{2}}{4m_{p}m_{a}}.$ (12)

where in the case of equality $E_{p}$ is the minimum energy that the
projectile (a plasma atom) of mass $m_{p}$ must posses to produced physical
sputtering on a target surface atom of mass $m_{a}$ that is bonded to the
surface with energy $E_{sb}$.

####  II.2.2 Formation energy for interstitial hydrogen and its isotopes

A similar design is used to provide an estimate of the interaction between a
hydrogen isotope (hydrogen, deuterium or tritium), denoted as $H$ for the sake
of simplicity, and the host crystal structure ($M$). The solubility of
hydrogen isotopes in a perfect crystal structure besides its porosity, is of
particular relevance inside the reactor to estimate hazardous tritium
retention. The formation energy of an interstitial hydrogen isotope $\Delta
E^{f}(H_{i})$ in a perfect crystal at 0 $\mathrm{K}$ can be calculated as
follows:

$M_{x}+H\rightarrow M_{x}H_{i}$

and

$\displaystyle\Delta E^{f}(H_{i})$
$\displaystyle=E^{DFT}(M_{x}H_{i})+E^{ZP}(M_{x}H_{i})$ (13)
$\displaystyle-E^{DFT}(M_{x})-E^{ZP}(M_{x})$
$\displaystyle-E^{DFT}(H)-E^{ZP}(H)$

where $E^{DFT}$ is the ground-state total energy calculated with a DFT code,
$E^{ZP}$ is the zero-point energy (ZPE). This expression can be approximated
by considering $E^{ZP}(M_{x}H_{i})\approx E^{ZP}(M_{x})+E^{ZP}(H_{i})$.
Therefore, Eq. (13) can be re-written as reported in Ref. [58]:

$\displaystyle\Delta E^{f}(H_{i})$ $\displaystyle\approx
E^{DFT}(M_{x}H_{i})+E^{ZP}(H_{i})$ (14)
$\displaystyle-E^{DFT}(M_{x})-E^{DFT}(H)-E^{ZP}(H)$

As demonstrated in Refs. [59, 60, 61], the solubility of a gas within a
crystal under low pressure conditions ($<$ 100 atm) can be derived for a
system in thermodynamic equilibrium between an ideal gas and the bulk of the
perfect crystal. It is expressed as a function $\Delta E^{f}(H_{i})$, along
with additional entropic corrections ($s_{corr}$). Specifically:

$\displaystyle\frac{N_{H_{i}}}{N_{i}}\approx exp-\left(\frac{\Delta
E^{f}(H_{i})-s_{corr}}{K_{b}T}\right)$ (15) $\displaystyle
s_{corr}=Ts^{vib}(H_{i})+\Delta\mu(H)(T,p_{0})+K_{b}T\ln\frac{p_{H}}{p_{0}}$

where $N_{H_{i}}$ and $N_{i}$ are the number H atoms and interstitial sites
respectively, $Ts^{vib}(H_{i})$ is the vibrational entropy of a single $H$
atom in an interstitial site, $\Delta\mu(H)(T,p_{0})$ and the last factor
account for the entropic contribution of $H$ at temperature $T$ and pressure
$p_{H}$. Eq. (15) is derived by the separating vibrational and conformational
entropy of the crystal with an interstitial hydrogen isotope ($WH_{i}$), as
the sum of the entropies of the perfect crystal ($W$) and the hydrogen
interstitial ($H_{i}$), while assuming $N_{i}\gg N_{H_{i}}$. To maintain
computational affordability for screening purposes, we calculate $\Delta
E^{f}(H_{i})$ neglecting zero-point energy (ZPE) and entropic contributions,
as outlined in Eq. (16). Consequently, under this formulation, the computed
interstitial formation energy remains equivalent for any hydrogen isotope.
However, it is worth noting that within the context of PFMs, a more accurate
approach would require to redefine the hydrogen reservoir as the hydrogen
stored in the crystal surface at equilibrium, resulting from the impingement
of high-energy plasma particles.

A second semi-automatized Python workflow, called the Hydrogen Interstitial
Formation Energy (H-IFE) workflow, is built to calculate $\Delta
E^{f}(H_{i})$. The main challenge of this calculation lies on the
identification of the interstitial hydrogen ($H$) sites with the lowest
formation energy. Below, are the steps used to accomplish this task:

  1. 1.

A supercell of the primitive crystal structure of a given material is created
with a minimum size of 9$\times$9$\times 9$ $\mathrm{\SIUnitSymbolAngstrom}$3.
The supercell is relaxed and its $E^{DFT}(M_{x})$ is calculated.

  2. 2.

$N$ supercells of the same size of step 1 but containing one interstitial H
atom are generated. The interstitial sites are randomly chosen from within the
primitive cell of the supercell, with the only constraint being that the
interstitial site of the $n+1$ supercell is at a minimum distance of 1
$\mathrm{\SIUnitSymbolAngstrom}$ from all crystal atoms and the interstitial
sites of all previous supercell. The number of interstitial sites, $N$, is
determined based on the size of the primitive cell, and the random search for
interstitial sites is terminated if more than 5,000 attempts are required.
This approach ensures that the primitive cell of each crystal structure is
uniformly explored in terms of interstitial site density. The interstitial
supercell are then relaxed and their total energy is calculated.

  3. 3.

$\Delta E^{f}(H_{i})$ is calculated from the total energy of the interstitial
supercell with lowest $E^{DFT}(M_{x}H_{i})$, using the following equation:

$\displaystyle\Delta E^{f}(H_{i})$ $\displaystyle=E^{DFT}(M_{x}H_{i})$ (16)
$\displaystyle-E^{DFT}(M_{x})-E^{DFT}(H)$

The same DFT code and set up of the input parameters used for the SBE workflow
are used for the H-IFE. For comparison with literature reference $\Delta
E^{f}(H_{i})$ is also computed for a hydrogen molecule ($H_{2}$) reservoir
($\frac{1}{2}E^{DFT}(H_{2})$).

##  III Ranking results

The MPDS database contains more than 500,000 crystal structures. When only the
materials possessing at least one out of the ten properties previously
selected are retained and radioactive and high pressure crystals are
discarded, the total number of available candidates is reduced to 7,056.

The $d^{ELMs}$ thickness ranking procedure is applied to the 177 crystal
structures having the required properties. These materials are ranked as a
function of their maximum allowable thickness to avoid melting during type-I
ELMs, and it is found that only 20 materials meet the set prerequisites. Out
of the 6879 remaining materials, 176 had the properties to be ranked on
$d^{SS}$. Of these, 13 candidates are retained due to their ability to
withstand steady-state conditions without melting. The remaining available
structures are, instead, classified using the Pareto front procedure. From the
latter, we identify 38 best-performing candidates, bringing the total number
of selected materials from MPDS to 71. Fig. 1 illustrates a schematic
representation of the material space in MPDS and its exploration, taking into
account the number of available properties and the corresponding ranking
procedure. All 7,056 candidates are then ranked according to the comparative
ranking scheme, to show consistency across different ranking procedures and
identify any candidates that exhibit discrepancies among the various rankings.

We then assess the potential of the top candidates thanks to the analysis of
literature on divertor materials. This serves as an additional benchmark to
validate our findings and identify the materials that have already been
explored and the corresponding advantages and disadvantages of their
application. New $d^{SS}$ and $d^{ELMs}$ values are calculated for the most
promising PFM candidates identified through the literature review. These
calculations are based on the literature values of the properties required
from heat balance equation (Eq. 4). Then, SBE and H-IFE workflows are applied
to highlight possible drawbacks of the selected structures. Ultimately, a
comprehensive ranking is provided to identify the most promising PFMs. Fig. 2
depicts the flowchart encompassing the entire process, beginning from the
retrieval of MPDS data to the generation of a list comprising the top
candidates for PFM. It should be noted that the chemical formulas of the
materials presented in the tables and figures of this chapter correspond
exactly to those documented in the MPDS database.

Figure 2: Flowchart illustrating the PFM selection process, which involves all
the steps described in Section II and the literature review. It starts from
the MPDS materials and narrows down to a set of selected PFM candidates. The
material selection is visually represented using a funnel chart.

###  III.1 Thickness ranking during ELMs

In Fig. 3 and Tab. 8 in the Appendix, we present 20 PFM candidates able to
sustain the divertor conditions during type-I ELMs without melting and
satisfying Eq. (5). These selected materials can be summarized into 3 main
categories.

  1. 1.

Carbon-based materials: diamond (C dia).

  2. 2.

Pure transition metals: W, Mo, Ir, Cu, Os, Ru, Rh, Re, Ta, Cr.

  3. 3.

Ceramics: (carbides) TaC, HfC; (nitrides) AlN, GaN, TiN; (borides) HfB2, TiB2,
ZrB2.

This list includes the most extensively studied PFMs proposed as divertor
materials in numerous experimental projects. In addition, other materials,
such as TiB2 and TiN, which have been employed as PFMs in tokamak limiters [9]
subject to comparable thermal loads as divertor materials, appear in the
selection.

It might result surprising that a few of the most well-known PFMs are not
included in this initial selection. For example, a few complex alloys, such as
stainless steel and Inconel [9], that have been utilized in past as PFMs do
not appear just because they are not included in the MPDS database. However,
these alloys do not withstand the high heat loads present in the divertor
region. On the other hand, while Ti and Be are listed in the MPDS, not all the
required properties are available to calculate $d^{ELMs}$. However, Be
exhibits a thermal conductivity that falls just short of satisfying divertor
conditions [62], while the thermal conductivity of pure Ti is significantly
inadequate for the intended purpose [63]. Another group of materials
previously employed as PFMs but excluded through thickness ranking are TiC,
pure B, Al, and boron carbide (B4C) [9], as they are designed to withstand
lower temperature regimes. Graphite and SiC-3C are excluded due to the
underestimate of their actual thermal property values in the MPDS database.

Figure 3: $d^{ELMs}$ logarithmic scale for the 20 selected PFM candidates
extracted from MPDS. $d^{ELMs}$ defines the maximum thickness for a material
to withstand type-I ELMs without melting. Only materials with $d^{ELMs}>0$
$\mathrm{mm}$ are selected.

###  III.2 Thickness ranking in steady-state conditions

Once the materials ranked with $d^{SS}$ are excluded, we identify 13 materials
with $d^{SS}\geq$ 3 $\mathrm{mm}$; these are reported in Fig. 4 (and Tab. 9).
A comparison of their characteristics with those of the $d^{ELMs}$ ranking
candidates reveals that most of these materials rank at the bottom of the list
in terms of performance, and they could potentially be excluded altogether
unless unreasonably high heat capacities are considered. Additional literature
investigations substantiate this hypothesis, except for three materials, W2C
and BN, which possess adequate melting temperatures, and VB2, which exhibits
sufficiently high thermal conductivity, warranting further scrutiny.

Figure 4: $d^{SS}$ defines the maximum thickness for a material to withstand
steady-state condition without melting. The 13 materials showing $d^{SS}\geq$
3 $\mathrm{mm}$ are selected from MPDS once materials ranked through
$d^{ELMs}$ criteria are excluded.

###  III.3 Pareto ranking

Pareto ranking is performed with materials whose properties cannot be assessed
by the two thickness rankings. Of the 6,703 crystal structures analyzed, 6 are
in the $3p$ category (Tab. 10), 6 in the $2p$ category (Tab. 11), and 26 in
the $1p$ category, including 5 polymorphs of BC2N and C3N4 which are excluded
from further analysis because they exhibit highly similar thermal properties
(see Table 12). We increase the portfolio of properties for the selected
materials with values found in the literature. Comparing these properties with
those of materials selected through $d^{ELMs}$ and $d^{SS}$ rankings, we find
that the number of remaining candidates qualifying as PFM is reduced to the
following:

  1. 1.

Metal alloys: Ta_0.5W_0.5.

  2. 2.

Ceramics: BC2N, C3N4, HfN, SiC-2H, Ta2C, TaN, VC, VN, WC, W2C.

  3. 3.

Carbon-based materials: graphite intercalated with CoCl2, FeCl3, SbCl5.

  4. 4.

Boron-based materials: Mg2B25C4, YB66.

  5. 5.

Perovskite nitrides: TiFe3N, VFe3N.

Among these materials, it is worth noting that the x-Fe3N perovskite nitrides
(TiFe3N and VFe3N) exhibit a high value for the isothermal bulk modulus.
However, the literature lacks a comprehensive analysis of their thermo-
mechanical properties, which is crucial for determining their applicability as
PFMs. For these reasons, they are not further considered as potential PFM
candidates.

###  III.4 Comparative ranking

The comparative ranking procedure is used as benchmark for the previous
ranking processes. From the weighted average win-fraction definition
($\overline{WF_{i}}$), materials belonging to the subsets $3p$, $2p$, and $1p$
naturally separate, with their scores limited to 100%, 66.7%, and 33.3%
respectively. The results of the comparative ranking procedure are presented
in Fig. 5, illustrating the materials belonging to category $3p$ ranked from
the 1st to the 100th position. Additionally, the Appendix includes materials
from the 101st to the 200th position of category $3p$, together with the top
200 materials from category $2p$ and $1p$ (Fig. 6 \- 10). $\overline{WF_{i}}$
is compared with the result coming from the thickness and Pareto ranking
procedures. For the latter, a binary value is assigned to every material in
order to show its positive or negative detection by the ranking process. In
order to emphasize the comparison between different ranking procedures, each
of these is applied to the entire set of available materials, regardless of
whether they are excluded by more stringent rankings.

The first 63 materials (out of 218) outperforming the $3p$ $\overline{WF_{i}}$
ranking scale (Fig. 5) contain the whole set of the PFM candidates able to
sustain divertor ELMs conditions, as indicated by the heat balance equation.
Of the top 15 materials, 12 of them are selected by the $d^{ELMs}$ ranking.
All of these candidates are also selected by the $d^{SS}$ ranking and one or
more Pareto rankings, indicating an excellent consistency between the three
ranking procedures. The importance of a comparative scheme is further
emphasized by its ability to avoid false negatives in screening studies of
this nature. Inconsistencies among different selection criteria can bring to
light previously dismissed candidates. Indeed, the two materials erroneously
discarded according to the $d^{ELMs}$ ranking (as indicated in Sec. III.1),
graphite and SiC-3C, are well-captured by the whole set of Pareto rankings
when included. Consistently, other known refractory materials exhibiting
excellent thermal properties are selected by lower level ranking procedures
when included: TiC, pure B, BeO, ZrC, B4C (as B13C2), and Li2TiO3. In
particular, the latter material has been proposed as a tritium breeding
material for the test blanket module of ITER [64]. Furthermore, another class
of materials recently proposed as W-alternative plasma-facing materials is
identified by examining the top 50 materials ranked by the 3p
$\overline{WF_{i}}$ and in Fig. 5; this is the class of MAX phase materials,
represented by Ti3SiC2 [65, 66]. Despite being these materials, and especially
Ti3SiC2, refractory and possessing a good thermal conductivity, their
decomposition temperature severely limits their suitability as divertor PFMs
[67]. In the $2p$ and $1p$ $\overline{WF_{i}}$ ranking scales, which consist
of 803 and 6,063 materials respectively, the top 200 candidates include the
complete set of $d^{SS}$ and $2p$ Pareto top candidates from the $2p$ set, and
20 out of 26 $1p$ Pareto candidates from the $1p$ set. The clustering of the
previously selected candidates near the top of the comparative ranking
provides additional evidence of the reliability of the ranking procedure we
have established.

Figure 5: Heatmap representation of the materials ranked from 1 (top left) to
100 (bottom right) according to the comparative ranking procedure. The
materials are classified into the $3p$ subset and arranged in order based on
their weighted average win-fraction ($\overline{WF_{i}}$) scores, with a
maximum potential score of 100%. The figure also includes the individual win-
fraction ($WF_{i,k}$) values for each property of the materials.

##  IV Literature review of the PFM candidates

We analyse the four main macro-categories of materials, identified using the
previous ranking criteria, based on previous literature results:

  1. 1.

carbon-based materials;

  2. 2.

boron-based materials;

  3. 3.

transition-metal-based materials;

  4. 4.

ceramics, which further divide into three main classes: borides, carbides, and
nitrides.

Among the most promising materials, we select specific prototypical crystal
structures as the foundation for calculating the formation energy of two
commonly observed defects, the surface binding energy and the hydrogen
interstitial formation energy, as described in Sec. V.

###  IV.1 Carbon-based materials

The ranking procedure highlights two distinct allotropes of pure carbon:
diamond and various types of intercalated graphite. Surprisingly, pure
graphite is discarded by the $d^{ELMs}$ ranking because of its low assigned
KDE thermal conductivity value (see Fig. 6). In fact, the majority of the
thermal conductivity values found in the MPDS database for graphite pertain to
the direction perpendicular to its hexagonal plane. However, in the case of
polycrystalline samples, such as those used for this applications, the primary
contribution to thermal conductivity arises from heat transport parallel to
the basal plane [68], which is thousands of times higher than the
perpendicular one [69, 63]. The resultant thermal conductivity of these
samples, accounting for the limiting factors inherent to the polycrystalline
averages, remains of the order of hundreds of
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ [68],
considerably higher than the value extrapolated from the database (5.7
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$) and
sufficient to avoid melting during type-I ELMs. Consequently, graphite can be
considered as a PFM candidate according to the $d^{ELMs}$ ranking.

Graphite and carbon fiber composites (CFCs) were identified as the primary
choices for (PFMs) for the ITER divertor [70]. Due to their remarkable thermal
conductivity and high sublimation temperature [71], these materials have
demonstrated the ability to withstand significant thermal loads with no
structural damage or, at most, minimal degradation of the polyacrylonitrile
fibers in the case of CFCs [72]. Furthermore, graphite-based materials pose
low concerns regarding plasma core contamination due to their low-Z number,
which minimizes radiative energy losses and allows for higher tolerance levels
within the plasma. For this reason, fusion experiments worldwide continue to
adopt graphite as the preferred PFM [7]. However, a material composed of low-Z
elements is subject to surface erosion by physical sputtering processes, with
the additional drawback of an enhanced chemical sputtering for graphite, due
to the high reactivity between carbon atoms and the hydrogen isotopes
contained in the plasma [73]. This leads to the formation of co-deposited
layers of tritium-containing hydrocarbon on all in-vessel components, which
increases the inventory of tritium [74]. In addition, the high thermal
conductivity of graphite deteriorates rapidly under the influence of energetic
neutrons [75, 76], an effect only partially mitigated by defects recombination
that occurs at high-temperature.

Graphite-based materials have been tested for many years [77] and, even
recently, carbon-based candidates such as diamond or diamond-like carbon
structures are re-considered and proposed as alternatives to W [78, 79].
However, effective suppression of the surface sputtering process, has not been
demonstrated. Despite many attempts on suppressing chemical sputtering [80,
81, 82] and radiation enhanced sublimation [83], physical sputtering process
plays a crucial role per se. Based on the literature data from the early 2000s
[40, 9], one can see that even when considering only physical sputtering, the
tritium co-deposition rate constant in graphite drops to approximately one-
fifth of its original value. This reduction is inadequate to offset the
tritium inventory limit that exceeds the threshold by two orders of magnitude
[84, 74]. For carbon too be the main constituent of a candidate material, this
constraint cannot be circumvented unless: (i) significantly higher SBE values
are assumed for the surface carbon atoms of the candidate material, (ii) a
different preferential target element is considered for the sputtering process
(e.g. low-Z carbides), (iii) the kinetic energy of plasma particles is reduced
below the threshold value of carbon sputtering, (iv) new tritium desorption
techniques are proposed.

Carbon-based materials are, therefore, not considered optimal as divertor PFMs
in ITER-like devices. However, due to their impressive ability to withstand
high thermal loads and their well-documented characterization in the
literature, they are included as reference materials in the list of top PFM
candidates. In particular, carbon diamond (C dia) and graphite (C gra) are
chosen as representatives of the carbon-based materials category for the
purpose of calculating defect formation energies.

###  IV.2 Boron-based materials

The boron-based materials reported in this study belong to the family of
compounds known as icosahedron-based borides. These materials consist of a
main structure of B12 icosahedra, which may be interconnected through chains
of carbon (or other nonmetallic) atoms. The voids within this framework can
also accommodate metal atoms. Boron-based materials are selected through the
$1p$ Pareto ranking procedure, due to their outstanding Knoop hardness, and
fall into two classes. The first class includes RB66, with R being a metallic
atom. RB66 materials combine the physical and chemical properties of
refractory crystals with the typical electrical, optical, and thermal
properties of amorphous semiconductors. Despite their interesting stability
under high temperature and abrasive conditions, they posses a really low
thermal conductivity value because of their complex crystal structure [85,
86]. An example of the second class is Mg2B25C4, where, beside the presence of
metal in the boron structure, carbon acts as a bonding element between the
various icosahedral B12 structures. The materials of this class have similar
thermo-mechanical properties to the RB66 crystal class [87]. Boron carbide B4C
has one of the highest thermal conductivity within the second class, reaching
a value of approximately 28
$\mathrm{W}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{m}}^{-1}$ at room
temperature. Under neutron irradiation, this value is significantly reduced
[88], but it can be largely restored through high-temperature annealing,
especially in the case of 11B enriched samples [89]. However, this value of
thermal conductivity is insufficient to avoid melting during ELM event, as
confrmed by $d^{ELMs}$ ranking procedure correctly discarding of this material
from the list of potential candidates. It is worth noting that B4C is
identified by a compositional variation of boron carbide, B13C2, in MPDS. Due
to its limited thermal conductivity, this class of materials is deemed
unattractive for withstanding high heat loads and frequent thermal shocks and,
therefore, it is not included in the list of top PFMs candidates. It is worth
noting that in combination with a highly conductive substrate, such as
graphite, B4C has been proposed as protective layer of the main divertor armor
material [90].

We finally remark that the presence of 10B (approximately 20% of natural
elemental boron) poses significant limitations on the use of boron-based PFMs.
In fact, 10B is known to undergo transmutation through the (n, $\alpha$)
reaction, resulting in the production of helium and lithium atoms and the
formation of cavities or bubbles in the matrix and grain boundaries [91, 92].
This process leads to void-swelling, high-temperature helium embrittlement,
and compromises the structural integrity of the crystal due to progressive
loss of 10B. The negative effects of 10B transmutation can be mitigated by
using isotopically enriched 11B, which has a considerably lower transmutation
cross-section to thermal neutrons [92, 93]

###  IV.3 Transition metal-based materials

##### Tungsten

The future and current generations of PFM candidates for the divertor region
are W-based. Despite the fact that W outperforms all other alternatives tested
so far, many downsides in its application remain, and are reported in
literature. W can sustain the divertor thermal loads in a very narrow
temperature window compared to the operating temperature interval of the PFM
surface. In steady-state conditions the surface of the metal is close to its
ductile–brittle transition temperature (DBTT), 200 to 600
$\mathrm{\SIUnitSymbolCelsius}$, depending on the manufacturing processes, and
below the DBTT the reduced thermal stress resistance leads to intense crack
networks [94]. While W does not encounter any surface changes at heat loads of
10 $\mathrm{MW}\text{\,}{\mathrm{m}}^{-2}$, even during typical slow-transient
of 20 $\mathrm{MW}\text{\,}{\mathrm{m}}^{-2}$, it shows severe surface damage
and deep crack formations [95]. In addition, when the low recrystallization
temperature of W (roughly 1,300 $\mathrm{\SIUnitSymbolCelsius}$) is surpassed,
grains growth and consequent increase of DBTT, brittleness and surface
roughening, and decrease of thermal shock resistance and mechanical strength
are observed [96, 97]. Finally, power pulses produced during ELMs can result
into cracking of the surface, partial melting of the tiles edges, and
depletion of the PFM surface [94]; in particular, alarming levels of W erosion
are reached when ELMs are combined with slow-transient loads [98, 99].

The heat load is associated with a flux of plasma particles and their
implantation in the W armor, for a distance of a few nanometers [100]. The
energy of the plasma particles (approximately 10 $\mathrm{eV}$) is not
sufficient to physically sputter the surface atoms of the W armor, as it is
much lower than the W sputtering threshold energy (209 $\mathrm{eV}$ for D,
102 $\mathrm{eV}$ for He, at 20 $\mathrm{\SIUnitSymbolCelsius}$ [101]).
Furthermore, the weak chemical reactivity between W bulk, hydrogen isotopes
and He atoms does not lead to any additional chemical sputtering [102, 103].
Nevertheless, the erosion of the W surface is a reason of concern due the
surface morphology changes produced by the interaction between He atoms and W
crystal defects in the presence of hydrogen, which leads to the formation of
helium bubbles and the so-called a fuzz layer, a process that is not yet
completely understood [104]. Despite the well-known decrease of the effective
thermal conductivity of the nano-fiber layer and consequent ejection of melted
droplets [29], the effects on the net material erosion rate is still under
debate [103, 105, 106]. Plasma particles can be retained by the PFM, which is
problematic as mentioned above for CFCs. ELMs and fuzz formations seem to play
an opposite role in terms of tritium retention in the W wall: the damaged
surface produced by ELMs enhances its retention compared to steady-state
conditions [107]; in contrast, the open porosity of the fuzz layer serves as
an additional release channel for the hydrogen isotopes [100]. However, in the
presence of He defects in W can affect its retention capacity; dislocation
loops punched by helium bubbles are responsible for tritium trapping sites
[108] and the total retention amount can increase drastically.

In this scenario, high-energy neutrons also play a role and irradiated-W
samples showed more than double tritium retention compared with undamaged W
[109]. Vacancy clusters of the neutron-damaged materials can act as additional
traps for hydrogen isotopes [110]. 14.1 $\mathrm{MeV}$ neutrons are also
responsible for the degradation of the wall armor under operational conditions
via reduction of the thermal conductivity and embrittlement, due to defect
accumulation and transmutation. For both properties, the worsening of
performance dictated by the contribution of dislocations and voids is in large
part recovered at higher temperature, thanks to the annihilation of the
irradiation-induced defects [111, 112]. W transmutations leads to the
production of rhenium (Re) and osmium (Os) in different percentages as a
function of the neutron’s flux [113]. In ITER, divertor replacement is planned
after a displacement per atom (dpa) of 0.7 with 0.15% of Re content and less
than half of that for Os content [114]; in DEMO this could be up to 20 times
higher [115]. Thermal diffusivity is negatively affected by the presence of Re
in the W bulk (around -25% of unirradiated W at 1,000 $\mathrm{K}$) and
temperature changes seem to be negligible at the relevant concentration [116,
117]. In terms of ductility (and recrystallization temperature) Re addition to
W is beneficial [29], but in presence of Os and dpa above 0.6-1 their
concentration is enough to produce an intermetallic second phase precipitate,
which leads towards strong hardening and embrittlement of the PFM, regardless
of temperature variations [111]. Transmutation products dominate over the
neutron-induced defects for the strengthening effect as long the irradiation
dose is higher than 0.6 dpa, while the opposite is true below 0.3 dpa.

To summarize the aforementioned points, W is presently considered the most
viable option to withstand the expected heat loads at the divertor of a
tokamak reactor, securing its position among the most promising PFM
candidates. However, its performance under extreme thermal stress conditions
and its long-term structural integrity when exposed to the combined effects of
plasma and neutrons remain unresolved issues and crucial areas of ongoing
research in the field.

##### Molybdenum

The most studied high-Z divertor PFM alternative to W is molybdenum (Mo). Mo
was used as divertor PFM in the Alcator C-Mod tokamak reactor and its behavior
extensively analyzed [118]. While Mo has a melting temperature and thermal
conductivity lower than W, it shows some advantages when compared with W, such
as a lower transmutation rate and therefore reduced formation of precipitates
[113], and slightly lower H retention [119]. In addition, as in the case of W,
Mo thermal conductivity is barely affected by high temperatures [120] and it
can be restored to its high values thanks to annihilation of defects on
neutron-irradiated samples [121]. However, neutron irradiation poses a clear
drawback on the application of Mo as plasma-facing material. Both the US and
UK waste disposal rating classify a Mo armor exposed to divertor neutron
irradiation conditions as high-level-waste [26, 35, 34]. Despite this, plasma-
Mo interactions under divertor conditions are still under investigation, and
shows many similarities with plasma-W interaction. The erosion in steady-state
conditions is mainly due to impurities in the plasma that have mass
sufficently large to sputter Mo (i.e. Mo itself and oxygen). Even if this
erosion is predicted to be slightly higher than the case of W, it remains
tolerable (1.5 $\mathrm{mm}$/burn-year) [26]. As in the case of W, the
interaction with He ions results in the creation of fuzz nanostructures on the
surface exposed to the plasma [122]. The temperature window for their
formation is narrower and lower than for the case of W [123], but it seems to
follow the same pattern [122]. Also, under unmitigated type-I ELMs-like pulses
the difference in erosion between the two metal with a pristine surface has
not been fully elucidated yet. A recent study reported that the volume of Mo
ablated from the surface by a single laser pulse is 10 times larger than W
[124], whereas another study reports that after 100 laser pulses (of similar
energy) the two materials do not show appreciable difference in the amount of
mass loss [125]. In both experiments the energy required from the laser pulse
for melting the two metals was found to be above the energy of mitigated
type-I ELMs, 1.0 $\mathrm{MJ}\text{\,}{\mathrm{m}}^{-2}$ for Mo and 1.4
$\mathrm{MJ}\text{\,}{\mathrm{m}}^{-2}$ for W [125]. Similar effects are
produced by pulsed lasers, simulating unmitigated ELMs on fuzz form Mo
surfaces [126]. When covered with a fuzz layer, the mass loss appears to be a
linear function of the laser pulse, instead of exponential as for the pristine
surface [125]. Under mitigated type-I ELMs, the fuzz layer gradually
disappears, without any melting, with a mass loss that increases as a function
of the surface temperature [103]. The extent to which this phenomenon, when
combined with arcing, could pose a risk of core plasma contamination remains
unclear. Additionally, it should be noted that a higher concentration
tolerance is predicted for Mo (${10}^{-3}$-${10}^{-4}$ [127, 128]) compared to
W ($5\text{\times}{10}^{-5}$ [129]) before plasma quenching. Overall, Mo is
still considered an attractive alternative to W and is therefore included in
the list of top PFM candidates for further investigations.

##### Copper

A high-strength copper alloy material is envisaged for the cooling tube system
running through the central region of the W monoblock of the ITER divertor
target. Divertors designed for DEMO envisage the use of pure copper (Cu) as
stress-relieving interlayer between W and the Cu-alloy [130, 131]. In fact,
despite its impressing thermal conductivity, Cu cannot be exposed directly to
the plasma flux, primarily because of its low hardness and softening
temperature [132] and the very intense erosion to which it is subject to under
hydrogen plasma flux [133]. In addition, even if Cu and Cu alloys are known to
maintain high thermal conductivities at elevated temperatures [134], the
exposure to neutron radiation can damage their crystal structure and cause
transmutations, which can negatively affect their thermal diffusivity [135,
136]. Furthermore, the induced radioactivity caused by the high-energy fusion
neutrons raises concerns about safety and waste disposal considerations [34,
35]. Due to these reasons, pure Cu is not directly employed as PFMs in fusion
tokamak devices and is not included in the list of promising PFMs. However,
application of Cu in the form of alloy as a thermal conductor material that is
not exposed to plasma has been extensively investigated, and it is important
to mention that, in a small number of studies, copper-alloys are also proposed
as primary divertor and first-wall materials, in particular W-Cu [133, 137]
and Cu-Li [138, 139, 140].

##### Platinum-group and noble metals

A set of transition metals, belonging to the so-called platinoids and noble
metals classes, are reported here for their valuable thermal properties. We
consider in descending order of performance (see Fig. 3): iridium (Ir), osmium
(Os), ruthenium (Ru), rhodium (Rh), rhenium (Re), and silver (Ag). The main
limit to their use as pure metals is dictated by their high price (they
typically range in cost from $1,000 to $100,000 USD per kilogram) and scarcity
[141] and their thermo-mechanical properties, typically slightly inferior than
more affordable alternatives (W and Mo). The main application of these
elements is thus as dopants or alloying elements. Re, being one of the
cheapest of this list, second only to Ag, has been proposed as alloying
element for W due to its enhanced thermo-mechanical performance: higher
ductility, high-temperature strength, plasticity, and corrosion resistance,
stabilization of the grain structure with increased recrystallization
temperature and weldability, and reduced degree of recrystallization
embrittlement [29]. However, the high cost and scarcity of Re is already
affecting other industrial sectors (platinum–rhenium catalysts and high-
temperature superalloy for jet engines) and strategies to reduce its
application, finding alternatives and recycling processes are under
investigation [142, 143]. Despite its cost-related limitations, Re is included
in the list of leading PFM candidates as the most favorable alternative within
this class of materials, in terms of cost and performance. The least expensive
Ag is discarded due to inferior thermal properties when compared to the
aforementioned transition metals, primarily dictated by a low melting
temperature.

##### Early transition metals

Together with the already discussed Mo, W, and Re, there are a few more
elements from the first 5 groups of the transition-metal class that exhibit
notable thermal capabilities, in particular Ta and Cr. Ta was suggested as a
potential PFM in the past: however two significant drawbacks are present: (i)
its inferior thermal properties when compared with W [144] (ii) the exothermic
hydrogen absorption process, which yields larger tritium retention than other
metals (such as W) with positive heat of solution for hydrogen [145]. In
contrast, Ta demonstrates a higher energy threshold for the formation of a
”fuzz” structure when exposed to divertor-like He ion flux, which suggests
that surface thermal and mechanical properties may not degrade as quickly as
for W in extreme fusion environments [146]. Chromium has been recently
proposed as PFM for the DEMO-divertor structural body of the monoblock.
Cr-10%W displays superior mechanical properties to Cr at low temperatures,
where W is brittle, and higher strength compared with pure Cr. Further
investigations of its behaviour under actual divertor conditions must thus be
carried out [147]. For this reasons, both Ta and Cr are added to the list of
promising PFM candidates.

##### Transition metal alloys

Only one material stands out in this category, exhibiting features comparable
to those of the positive thickness materials under type-I ELMs, this is the
Ta_0.5W_0.5 solid solution, possessing high level thermal conductivity, heat
capacity, and isothermal bulk modulus, with a melting temperature of 3,450
$\mathrm{K}$ is reported in the literature [148]. The two constituents of this
alloy exhibit complete solubility [149], but given the low mass contrast
between its components, the phonon contribution to its thermal conductivity is
not significantly diminished, in contrast to many other alloy systems.
Previous studies also investigate the thermo-mechanical properties of Ta-W
composites, wherein Ta fibers are dispersed within a nanostructured W matrix
to enhance fracture toughness [150]. However, the ability of these composites
to withstand plasma irradiation conditions remains a topic of ongoing debate
[150, 151, 152]. Our screening process did not identify any other promising
alloys; nevertheless, one of the most promising research areas in terms of
alloys is represented by the high-entropy alloys (HEAs). These materials are
currently not included in MPDS but, in a recent study [153], a W-Ta-based HEA
alloy (38% W, 36% Ta, 15% Cr, 11% V) was able to show outstanding neutron
radiation resistance properties, high hardness and small irradiation-induced
hardening [153]. The material capacity to sustain high-dose neutron
irradiation is attributed to the equal mobilities of vacancies and self-
interstitials, which maximizes recombination and reduces defect concentration.
This is likely due to the rough energy landscape induced by local lattice
distortion and the disparity in chemistry. While lattice distortions can
affect the thermal conductivity (a common drawback of this class of materials
[154]), very recent promising results for thermal conductivity have been shown
for the class of medium-entropy alloy [155]. As a prototype of this class, the
ordered cubic structure (B2-superstructure [156]) of TaW is added to the list
of top PFM candidates for which defect formation energies are then calculated
using the DFT workflows developed here.

###  IV.4 Ceramics

Most of the ceramics identified by the three ranking criteria can be
classified as ultra-high temperature ceramics (UHTCs). This family includes a
variety of materials, such as borides, carbides, and nitrides of early
transition metals, and are characterized by their extremely high melting
temperatures, typically above 3,000 $\mathrm{\SIUnitSymbolCelsius}$, and/or
their ability to withstand continuous use at temperatures exceeding 1,600
$\mathrm{\SIUnitSymbolCelsius}$ [157]. One of the major challenge for UHTCs in
thermostructural applications is their poor resistance to thermal shocks.
However, the manufacture of fibrous monolith ceramics limit the formation and
propagation of cracks to the cell boundaries and can partially address this
problem [158].

##### Borides

Refractory diborides, such as the here HfB2, TiB2, and ZrB2 reported here,
exhibit an unusual combination of high strength, moduli, melting temperature,
and thermal conductivity. The possibility of using diborides as PFMs is now
again under discussion, since the production of isotopically separated 11B
diborides is now demonstrated [159, 93]. In fact, neutron irradiation of
naturally occurring B is indicated as the primary causes of swelling and
catastrophic cracking of borides, due to the helium and lithium accumulation
deriving from the (n, $\alpha$) reaction of 10B, as showed in Sec. IV.2.
Neutron-irradiation experiments on Ti11B2 demonstrated reduced swelling
(especially at high temperature), suppressed macroscopic damage and moderate
micro-cracking formation [93]. Concerning the thermal properties, ZrB2 shows
the capability to maintain a level of thermal conductivity at high temperature
that is comparable to W. Furthermore, since the electronic contributions
mostly determine the total thermal conductivity value, neutron-induced defects
are expected to have a minimal effect [159]. Both TiB2 and ZrB2 show
encouraging results in terms of erosion resistance under high-energy He ions:
the TiB2 sputtering yield is reduced when compared with those of pure Ti [160,
161], while ZrB2 show comparable erosion depth with that of W [159]. ZrB2 also
shows an interesting oxidation mechanism that might lead to oxygen absorption
at higher temperatures [157, 158, 159]. However, experiments at divertor
conditions are still lacking. While showing the same oxidation mechanism as
ZrB2, HfB2 has raised some concerns regarding the production of 20
$\mathrm{MeV}$ gamma rays when irradiated with 14 $\mathrm{MeV}$ neutrons
[162].

Inspired by the growing field of high-entropy alloys (HEAs) where five or more
elements are mixed in similar proportions to form a new alloy, ZrB2-based
ceramics are also proposed as high- and medium-entropy ceramics (HECs) [163],
and their thermo-mechanical properties are investigated [164].
(ZrxTiyNbyTay)B2 systems containing Zr, Ti, Nb, and Ta at the highest mixing
entropy posses the highest hardness and moduli, but also the highest phonon
scattering with the lowest thermal conductivity. Similarly, carbon-boron-
titanium compounds with varying B/C/Ti ratios have been studied. It is
observed that the highest thermal shock resistance is achieved for
compositions where the constituents are mixed in similar proportions, which
also exhibited lower thermal conductivity compared to pure borides and
carbides [165]. ZrB2 is also proposed as alloying ceramic to produce W-ZrB2
alloys, with improved bending strength and hardness thanks to the
strengthening of the originally fragile W grain boundaries [166]. An
additional diboride, VB2, is also of interest, according to the $d^{ss}$
ranking, due to its remarkable values of high melting temperature and thermal
conductivity. Considering the combination of these properties, with the
noteworthy volumetric heat capacity reported in the literature [167, 168]
(exceeding 5 $\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{cm}}^{-3}$
at high temperatures), VB2 has been previously suggested as a protective layer
to mitigate the ingress of hydrogen isotopes into structural materials in
fusion devices [169, 170]. Consequently, the crystal structures of all the
presented diborides, namely HfB2, TiB2, ZrB2, and VB2, are included in the
compilation of the most favorable PFM candidates.

##### Carbides

Carbides of the same transition metal have usually larger melting temperature
and hardness than the corrisponding borides and nitrides [171], at the price
of a substantially reduced thermal conductivity when compared with borides
[172]. Despite obvious concerns associated with high carbon content, linked
the high level of tritium retention caused by carbon sputtering and co-
deposition, carbides have mainly been investigated as PFMs for doping or
alloying W. Three polymorphs of SiC are reported in MPDS, but only the
hexagonal polytype, SiC-2H, is selected according to the Pareto ranking thanks
to its excellent melting temperature, around 3,100 K. In fact, the cubic
polytype SiC-3C is discarded by the thickness (ELMs) ranking, due to a
reported very low thermal conductivity. This is quite surprising, since SiC-3C
is known for its outstanding thermal conductivity (SiC-3C possesses the
highest thermal conductivity value [173], along with a comparable melting
temperature and slightly inferior mechanical properties [174] than the other
polymorphs). Indeed, an erroneous value of thermal conductivity is reported in
the MPDS database, 0.49
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$, instead of
the correct value, 490
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$. Once this is
accounted for, SiC-3C appears at the top of the $d^{ELMs}$ ranking. SiC is
largely investigated as refractory and plasma-facing material due to its
excellent thermal properties and low activation under fusion or fission
neutron bombardment [175]. It is also considered as divertor PFM [176],
motivated by its favourable erosion properties when compared with graphite. In
particular, SiC-3C has been used as a coating PFM of graphite tiles during the
2017-19 DIII-D run campaigns [176], and the extrapolation of results to a
DEMO-type device shows a 70% decrease of the erosion yield of carbon is
expected with SiC as PFMs at the divertor with respect to graphite. The
reduced erosion was also confirmed in a later study [177]. This effect,
combined with the low solubility and diffusivity of tritium in SiC, [178]
associated with a lower predicted tritium inventory in the co-deposited layer
when compared with graphite (2/3 of the carbon material) [179], could motivate
further investigations of SiC as divertor PFM. However, the reduced erosion
has been extrapolated from very low concentration of C atoms in the plasma. In
addition, the value of thermal conductivity of 490
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ is not
realistic under operational conditions where high temperature and neutron
bombardment occur. The high temperature reduces significantly the thermal
conductivity of the most well-known SiC polytypes [180, 181]. In addition, the
presence of neutron-induced defects in SiC results in detrimental swelling and
a significant decrease of the thermal conductivity, as demonstrated in several
studies [182, 183, 184, 185]. Despite the fact that these defects are
partially recovered at high temperatures, the thermal conductivity values
obtained for $\beta$-SiC at dpa levels comparable to ITER/DEMO devices
decreases to around 1 to 10
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ (depending on
the irradiation temperature), making this material almost unsuitable for the
expected thermal loads. This issue is ameliorated in modern SiC fiber
composites that exhibit minimal degradation of their thermo-mechanical
properties even when subject to high levels of neutron damage [176, 186].
However, their thermal conductivity in the pristine state values are low, and
comparable to those of neutron-damaged monolithic CVD SiC. As a consequence,
based on the current state of the literature, SiC does not meet the criteria
to serve as a PFM in the divertor region.

With small percentages of TaC that fill up the grain boundaries, sintered W
shows increased microhardness as a function of the TaC content, no cracking
formation under ELMs-like thermal shock, and even increased thermal
conductivity at the low carbide percentage of 1% [187, 188]. Similar effects
on the mechanical properties, thermal stability, and thermal conductivity are
also reported for this W-HfC alloy [189]. However, the retention of hydrogen
isotopes is generally increased under divertor-like deuterium bombardment
[190] and can be one order of magnitude larger when compared with pure W in
the case of W-1.1%TiC and W-3.3%Tac at 800 $\mathrm{K}$, decreasing with
increasing temperature. This might arise concerns on the tritium inventory.
The direct use of TaC and HfC as PFMs has not been extensively investigated.
Despite concerns related to their high carbon content, their thermal
properties are at the threshold of meeting the requirements for withstanding
divertor conditions without undergoing melting. However, the pure alloy system
TaC-HfC, shows the highest melting point of the known materials at the right
ratio of Ta and Hf (4:1) [191, 192] and decent values of thermal conductivity
that, in contrast to pure metals or diborides, increase with rising
temperature [193]. Further studies, including tritium inventory and behaviour
under neutron irradiation, should be pursued, since at the experimentally
reported values of the melting temperature, heat capacity and thermal
conductivity [193], the 4TaC-1HfC alloy system would withstand type-I ELMs
heat loads at the divertor with an armor thickness of about 5 $\mathrm{mm}$.
Rather than for their individual potential, TaC and HfC are added to the list
of PFMs to calculate defects formation energies as a proxy of their ultra-high
temperature alloys.

Given the melting temperature, heat capacity, and the thermal conductivity
reported in literature [194], vanadium carbide (VC), can only withstand
steady-state conditions, contrary to its boride form. Furthermore, its high-
temperature stability is moderate, but its ability to act as a barrier for
hydrogen permeation is noteworthy, which has led to the exploration of various
type of vanadium carbides (V_xC_y) as a coating layer for structural fusion
materials [195, 196]. The small hydrogen diffusion is attributed to its high
hydrogen penetration barriers and large endothermic detrapping energy of H
atoms in C vacancies. Previous test of the V_xC_y behavior under neutron
irradiation reports hardening and increased electrical resistivity [197, 198].
Even if more investigation is required to understand degradation of the
thermal properties induced in the presence of a plasma flux, this material
does not qualify to sustain divertor heat loads and is discarded from the list
of top PFM candidates.

Ta2C, along with the other carbides emerging from the previous ranking
selection criteria - TaC, HfC, WC, and W2C - has also been proposed in the PMF
field, but solely as a dispersion agent to enhance mechanical properties,
suppress the embrittlement of W, and improve the migration of grain boundaries
[199]. Ta2C has been explored also as a bonding material for joining the
structural materials or as a cooling system to the W armor [200]. Although
Ta2C is not extensively investigated as an alternative to W, its promising
melting temperature raises interest. However, the current lack of available
literature on the heat capacity and thermal conductivity of this material
hinders the determination of its capability to withstand divertor conditions
and to compare it with the top PFMs. Therefore, we refrain from conducting
additional follow-up studies on Ta2C at this time.

The last class of ceramic materials that emerges from the Pareto set analysis
is that of tungsten carbides, in particular WC and W2C. As previously
mentioned, they might be applied as alloying compounds of a pure W armor [201,
202], but their thermo-mechanical properties and response under plasma and
neutron irradiation are considered good enough to be considered and studied as
an alternative PFMs at the divertor [203, 204] and therefore these are
retained in our ranking. At high temperatures, WC has a lower value of thermal
conductivity when compared with pure W (reaching around 2/3 of its value).
However, WC has a thermal shock coefficient that is 3 times higher, indicating
better resistance to thermal shocks, and a Vickers hardness that is higher by
a factor of 8 at room temperature to 20 at operating temperature. WC shows a
good tolerance to neutron irradiation when compared to other carbides in terms
of cracking resistance. Regarding the interaction with plasma, WC shows
slightly higher sputtering rate and hydrogen retention than W. W2C has been
mainly applied as composites compound of both W [202] and WC [205]. When W2C
is subject to He ion bombardment, contrary to W, He bubbles do not form inside
its crystal grains [202]. These bubbles are the early stages of the fuzz
layer. The improved mechanical properties and surface layer stability of
tungsten carbides are potentially offset by reduced thermal properties and
increased sputtering and hydrogen retention, which has led to ongoing debates
over their applicability as substitutes for W or as composites with the latter
metal.

##### Nitrides

The rhombohedral structure of boron nitride (BN) is selected according to the
$d^{SS}$ ranking, contrary to the cubic phase which is discarded due to the
$d^{ELMs}$ ranking. Other two polymorphs of boron nitride are well-known: the
hexagonal polytype, discarded during the data cleaning procedure due to
discrepancies in reported properties that matched another polymorph according
to references, and the wurtzite polytype, not documented in MPDS. The melting
temperatures and heat capacities of the polymorphs are in agreement with the
KDE value extracted from MPDS [206, 63], but significant discrepancies arise
in the thermal conductivities. Computed values of isotopically pure crystals
of BN are around 1,900
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ [207, 208],
while experimental measurements report values that fall within the range of
150-550 $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ [209,
210, 211, 212, 213]. Thermal conductivities of this order of magnitude would
justify the potential application of boron nitride as PFM, even at the
divertor region, contrary to what the KDE values suggest. Although the effect
of neutron bombardment on the thermal conductivity of BN is not extensively
investigated, thermal annealing shows limited impact on the neutron-induced
damage of an irradiated BN structure [214], indicating a potential
irreversible degradation of its thermal conductivity. Nonetheless, the
chemical and thermal stability of BN, its limited physical sputtering, high
thermal conductivity, and resistance to thermal shocks, has motivated its
application as first-wall protective layer in fusion devices in the past [215,
213]. The interaction of BN with the plasma requires further investigation. On
the one hand, cubic and hexagonal BN are highly promising for hydrogen barrier
applications within the refractory nitrides class, exhibiting a diffusion
barrier to hydrogen migration of approximately 3.0 $\mathrm{eV}$ [216]. On the
other hand, BN shows high hydrogen permeability and solubility, which may lead
to significant tritium inventories [217]. The ability of BN to efficiently
store hydrogen molecules is linked to specific structural characteristics,
such as nanotubes and hexagonal structures with high porosity and structural
defects [218]. While these features are intentionally pursued for certain
applications, they should be avoided when considering BN suitability as a PFM.
Given its intriguing thermal properties and the necessity to conduct further
investigations under operational conditions, BN is included in the list of
leading PFM candidates, specifically in its cubic polymorph (c-BN). Recently,
a material sharing the same crystal structure as c-BN, boron arsenide (BAs),
has emerged as the second-highest thermally conducting material among all
known isotropic crystals, after diamond [219, 220]. The highest experimentally
measured thermal conductivity for BAs single crystals is 1,300
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$. However, the
thermal stability of BAs is compromised as it decomposes into B13As2 at 1,200
$\mathrm{K}$, releasing arsenic gas under atmospheric pressure [221]. This
negates its application as PFM.

Another promising material, but quite unexplored for PFM application, is
represented by the hexagonal, cubic, rhombohedral, orthorhombic, and amorphous
boron-carbon nitride BC2N. This material posses thermal and mechanical
properties that lie between those of carbon diamond and boron nitride (c- and
h-BN). Cubic, hexagonal and rhombohedral BC2N show impressive Vickers hardness
value (greater than c-BN) [222, 223, 224, 225]. A thermal conductivity that is
almost comparable with carbon diamond is calculated for h-BC2N [224] and close
to that of c-BN is estimated for orthorhombic BC2N [226]. The Debye
temperature of both hexagonal and orthorhombic BC2N is estimated close or even
higher than diamond, with values between 1,700-2,000 $\mathrm{K}$ [227, 228].
The heat capacity of the orthorhombic phase is predicted to lie between that
of diamond and boron nitride [227, 229]. It can be assumed that the behavior
of this class of compounds under divertor conditions is similar to their well-
known precursors, as low-Z materials. However, differences are expected due to
their chemical reactivity, and specific studies are necessary to determine
their suitability as PFMs. For these reasons, the hexagonal phase h-BC2N is
added to the list of top PFM candidates.

The C3N4 compound has been proposed in the past as potential low-Z PFM
alternative to graphite and CFCs, and a possible reduction of hydrogen
diffusion and retention inside the material was hypothesized [230]. Given the
difficulties connected with its synthesis [231], its thermo-mechanical
properties are not well-known. Only estimates of its thermal conductivity are
reported in literature, which range from 150 to 500
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ as a function
of the phase considered [226, 232, 233, 234]. Bulk moduli rivaling that of
diamond are also predicted for this crystal [233, 235]. Given the interesting
thermo-mechanical features, the potentially reduced hydrogen inventory (a
bottleneck in the application of carbon-based materials), and that synthesis
paths to produce microcrystalline and noncrystalline C3N4 have been proposed
[231, 236], the hexagonal structure $\beta$-C3N4 is also added to the list of
top PFM candidates.

Vanadium nitride (VN) exhibits high heat capacity and moduli, but its inferior
thermal properties, in terms of melting temperature and thermal conductivity
[237], preclude it from any application as a PFM in a fusion tokamak reactor.

The case of aluminum nitride (AlN) can be compared to that of gallium nitride
(GaN). For both of these materials, scientific literature on their potential
use as divertor armor is limited (AlN) or absent (GaN). One of the biggest
limitation in terms of applicability as PFM is dictated by their decomposition
temperature in vacuum or at low N2 level. Group III nitrides start to
decompose in vacuum at temperatures that are in the order of one third of
their melting temperature, and for AlN and GaN are calculated to be 1,040 and
850 $\mathrm{K}$ respectively [238]. This sets a temperature limit for high-
temperature applications that is much lower than the estimate presented here;
therefore, they are discarded from the list of promising PFMs. Despite this,
AlN has also been proposed as a substrate for bolometric systems, which
provide information on the spatial distribution of radiated power from the
main plasma and the divertor region [239, 240]. It shows good impermeability
to hydrogen isotopes [241, 242] and low susceptibility of its composition
under low-energy H and He bombardment [243]. However, both neutron irradiation
and high temperature have a detrimental effect on its thermal conductivity
[244, 245, 246]. In contrast, GaN is mainly investigated for its electrical
and optical properties for semiconductor device applications [247] and its
thermo-mechanical behaviour under plasma and neutron bombardment is mostly
unknown.

HfN shows thermo-mechanical properties that are not dissimilar from those of
HfC [172], but a lower melting temperature which does not consent its
application as primary PFM. However, alloy or composites of this material
might be investigated in the future.

An ultra-high temperature ceramics that has only recently been identified as a
ultra-high thermal conductive material is the hexagonal $\theta$-phase of TaN
[248]. Its astonishing thermal conductivity, estimated between 800-1,000
$\mathrm{W}\text{\,}{\mathrm{m}}^{-2}\text{\,}{\mathrm{K}}^{-1}$ at room
temperature (above 200
$\mathrm{W}\text{\,}{\mathrm{m}}^{-2}\text{\,}{\mathrm{K}}^{-1}$ at 800
$\mathrm{K}$), is just below that of c-BN, and more than double than the
largest for a metal, silver. This might place $\theta$-phase TaN at the top of
the thickness ranking methods. To the best of our knowledge, its application
as PFM has not been discussed yet, and a careful and critical analysis of its
mechanical and thermal properties and their evaluation under plasma and
neutron irradiation is thus needed. As stated in the reference article [248],
phonons are the main heat carriers for this material, and the high thermal
conductivity of the $\theta$-phase TaN is caused by a particular combination
of factors: (i) weak electron-phonon scattering, (ii) low isotopic mass
disorder, (iii) large frequency gap between acoustic and optical phonon modes
[220] that, together with (iv) acoustic bunching, impedes three-phonon
processes (which is however dominated by four-phonon scattering). $\theta$-TaN
is also added to the list of top PFM candidates.

Thanks to its decent thermal properties, TiN is selected by the $d^{ELMs}$
ranking as one of the lowest performing material, but still able to sustain
divertor conditions in a thin layer. Consistently, TiN was used in a
stellarator fusion device as coating layer of multiple plasma-facing
components [249] and has been recently proposed for DEMO as tritium permeation
barrier of the first-wall in a sandwich-like structure [250]. It shows the
capability to withstand hydrogen plasma shots in IR–T1 tokamak with no
significant morphology defects [251] and its permeability to hydrogen is also
confirmed computationally [216]. Although it may not be considered a favorable
choice for a primary plasma-facing material, TiN has the potential to be used
as a coating layer on plasma-facing components. Therefore, it is included to
expand the list of the most promising candidates.

Finally, all the nitrides presented here have in common the material
activation issues associated with the presence of nitrogen. The transmutation
of nitrogen (14N) under neutron irradiation produces long-lived carbon
isotopes 14C [252, 253] that possess a time-to-low-level-waste (time-to-LLW)
in the order of thousands of years [35].

###  IV.5 Revisited thickness ranking

The collection of materials obtained through the two-step ranking procedure -
thickness (ELMs and SS) and Pareto ($3p$, $2p$, $1p$) rankings - has been
qualitatively evaluated with the indication of the comparative ranking, which
highlighted potential false negatives, and then with the review of the
literature for their application as PFMs. The finalized list of PFM candidates
is presented in Tab. 4. Their capacity to withstand heat loads is re-evaluated
by calculating $d^{SS}$ and $d^{ELMs}$ with heat properties extracted directly
from the literature references reported in the table. Both high and low
temperature properties are reported, in order to calculate the thickness
values at the two temperature extrema. Considering the two temperature ranges
reported in Tab. 4, high-temperature conditions are more representative of the
thermal demands that the designated PFMs in the divertor region must
withstand. Therefore, the candidates are ranked in descending order of maximum
$d^{ELMs}$ at high-temperature. The candidates with a larger maximum
thicknesses that can avoid melting are prioritized because their surface
temperature can be further decreased by reducing their thickness, while
keeping a minimum width that can sustain erosive processes.

Carbon-based materials are ranked at the top of the list, diamond and
graphite, followed closely by the ceramic materials C3N4 and BC2N, which are
novel in the field of PFM, and the well-known BN. Within the metal-based
materials, W shows the highest performance, followed by Mo and three
transition-metal borides: TiB2, ZrB2, and HfB2. The extensively researched WC
[203, 204] and the recently investigated $\theta$-TaN [248] exhibit adequate
thermal properties, with maximum thicknesses of approximately 4 $\mathrm{mm}$.
All the other materials possess $d^{ELMs}$ values $\leq$ 2.5 $\mathrm{mm}$,
significantly lower than the current minimum reported thickness design for a W
divertor of 6 $\mathrm{mm}$ [254]. Consequently, pure Ta, Re, VB2, TaC, Cr,
TaW, TiN are likely unsuitable for use as primary PFMs. However, they could
also be utilized as coating layers or dopants of the primary PFM.
Additionally, they might be incorporated as constituents within alloy systems
or composite materials. Only two materials within the selected set showed
$d^{ELMs}\leq$ 0 under both temperature conditions: W2C and HfC. Although
these materials do not exhibit exceptional performance, they are still
included in the list. While the former may be considered as a potential
contaminant of the high-performing WC, the latter represents one component of
the two prototype crystals that constitute the aforementioned 4TaC-1HfC alloy
system.

Table 4: List of the top 21 PFM candidates identified in this study. The four
properties required for the calculation of the heat balance equation (Eq. (4))
are retrieved directly from literature. Thermal conductivity $k$ and heat
capacity at constant pressure $C_{p}$ are reported at both low (around RT) and
high temperature. The actual values for low and high temperature may vary
depending on the conditions documented in the literature. The maximum
thickness to avoid melting in steady-state condition ($d^{SS}$) and during
ELMs ($d^{ELMs}$) is calculated at low ($T_{L}$) and high ($T_{H}$)
temperatures. Computational data are reported only when experimental measures
are not available.

Formula | $\bm{T_{m}}$ | Ref. | $\bm{k(T_{L})}$ | $\bm{T_{L}}$ | $\bm{k(T_{H})}$ | $\bm{T_{H}}$ | Ref. | $\bm{C_{p}(T_{L})}$ | $\bm{T_{L}}$ | $\bm{C_{p}(T_{H})}$ | $\bm{T_{H}}$ | Ref. | $\bm{\rho}$ | Ref. | MM | $\bm{d^{SS}(T_{L})}$ | $\bm{d^{SS}(T_{H})}$ | $\bm{d^{ELMs}(T_{L})}$ |  $\bm{d^{ELMs}(T_{H})}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
($\mathrm{K}$) | ($\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$) | ($\mathrm{K}$) | ($\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$) | ($\mathrm{K}$) | ($\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}\text{\,}{\mathrm{K}}^{-1}$) | ($\mathrm{K}$) | ($\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}\text{\,}{\mathrm{K}}^{-1}$) | ($\mathrm{K}$) | ($\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$) | ($\mathrm{g}\text{\,}{\mathrm{mol}}^{-1}$) | ($\mathrm{mm}$) | ($\mathrm{mm}$) | ($\mathrm{mm}$) | ($\mathrm{mm}$)
C (diamond) | 3,950(1)(2) | [255, 71] | 2,300 | 300 | 1,540 | 400 | [63] | 6.2 | 300 | 24 | 1,200 | [63] | 3.51 | [63] | 12.01 | 406.53 | 272.20 | 361.60 | 253.51
C (graphite) | 3,950(2) | [71] | 1,950 | 300 | 448 | 1,200 | [63] | 8.6 | 300 | 23 | 1,200 | [63] | 2.20 | [63] | 12.01 | 344.66 | 79.18 | 300.28 | 66.17
h-C3N4 | 3,500(3) | [235] | 520 | 300 | 520 | 300 | [232] | 30 | 300 | 80 | 1,200 | [256] | 3.53 | [235] | 92.06 | 80.21 | 80.21 | 53.39 | 63.78
h-BC2N | 3,500(3) | [229] | 2,000 | 300 | 400 | 1,000 | [224] | 28 | 300 | 84 | 900 | [227] | 2.90 | [222] | 48.84 | 308.50 | 61.70 | 264.74 | 50.40
c-BN | 3,240 | [63] | 550 | 300 | 350 | 650 | [210] | 20 | 300 | 42 | 1,200 | [257] | 3.49 | [210] | 24.82 | 77.69 | 49.44 | 60.04 | 39.73
W | 3,687 | [63] | 174 | 300 | 112 | 1,200 | [63] | 24 | 300 | 28 | 1,200 | [258] | 19.26 | [63] | 183.84 | 28.47 | 18.32 | 17.97 | 10.53
TiB2 | 3,193 | [63] | 65 | 300 | 70 | 1,300 | [259] | 63 | 300 | 83 | 1,200 | [260] | 4.95 | [63] | 69.49 | 9.03 | 9.72 | 4.23 | 5.38
Mo | 2,895 | [63] | 138 | 300 | 96 | 1,200 | [261, 120] | 22 | 300 | 30 | 1,200 | [120] | 10.23 | [63] | 95.94 | 17.11 | 11.90 | 7.44 | 4.99
ZrB2 | 3,323 | [63] | 58 | 300 | 65 | 1,300 | [259] | 45 | 300 | 75 | 1,200 | [262] | 6.17 | [63] | 112.85 | 8.43 | 9.45 | 2.31 | 4.43
HfB2 | 3,373 | [63] | 51 | 300 | 60 | 1,300 | [259] | 50 | 300 | 80 | 1,200 | [262] | 10.50 | [63] | 200.11 | 7.54 | 8.87 | 1.98 | 4.10
$\theta$-TaN | 3,363 | [63] | 70 | 300 | 70 | 300 | [263, 248, 264] | 41 | 300 | 41 | 300 | [265] | 13.70 | [63] | 194.95 | 10.32 | 10.32 | 4.10 | 4.10
WC | 3,022 | [63] | 120 | 300 | 70 | 1,000 | [204] | 36 | 300 | 52 | 1,000 | [266] | 15.60 | [63] | 195.85 | 15.64 | 9.12 | 7.48 | 3.94
Ta | 3,290 | [63] | 58 | 300 | 60 | 1,200 | [267] | 26 | 300 | 28 | 1,200 | [268] | 16.63 | [63] | 180.95 | 8.34 | 8.63 | 2.12 | 2.53
Re | 3,458 | [63] | 48 | 300 | 45 | 520 | [269] | 23 | 300 | 31 | 1,200 | [270] | 21.02 | [63] | 186.21 | 7.30 | 6.85 | 1.88 | 2.32
VB2 | 2,723 | [63] | 42 | 300 | 42 | 1,300 | [259] | 45 | 300 | 70 | 1,200 | [168] | 5.07 | [271] | 72.56 | 4.85 | 4.85 | 0.24 | 1.15
TaC | 4,153 | [63] | 22 | 300 | 22 | 300 | [272] | 30 | 300 | 52 | 1,200 | [273] | 14.65 | [63] | 192.96 | 4.11 | 4.11 | 0.19 | 1.13
Cr | 2,180 | [63] | 94 | 300 | 62 | 1,200 | [63] | 22 | 300 | 32 | 1,200 | [274] | 7.20 | [63] | 52.00 | 8.30 | 5.47 | 1.29 | 0.75
TaW | 3,450 | [148] | 54 | 673 | 60 | 1,400 | [275, 276] | 25 | 300 | 28 | 1,200 | SWA(4) | 17.78 | [156] | 364.79 | 8.19 | 9.11 | -0.20 | 0.74
TiN | 3,220 | [63] | 25 | 500 | 25 | 1,200 | [277] | 45 | 500 | 55 | 1,200 | [278] | 5.21 | [63] | 61.87 | 3.51 | 3.51 | 0.27 | 0.58
W2C | 3,073 | [63] | 30 | 300 | 30 | 300 | [279] | 59 | 300 | 77 | 1,000 | [280] | 14.80 | [63] | 379.69 | 3.99 | 3.99 | -0.57 | -0.001
HfC | 3,273 | [63] | 30 | 300 | 10 | 1,000 | [281, 278] | 37 | 300 | 49 | 1,000 | [278] | 12.20 | [63] | 190.50 | 4.29 | 1.43 | -0.20 | -0.82

(1) Diamond $T_{m}$ is the same as graphite due to graphitization process
happening above 1,900 $\mathrm{K}$. (2) Graphite sublimates at atmospheric
pressure. (3) $T_{m}$ of C3N4 and BC2N are defined between those of diamond
and boron nitride by looking at the predicted Debye temperatures. (4) SWA
indicates that the $C_{p}$ of TaW is calculated as a stoichiometric weighted
average of the $C_{p}$ of its metallic single-element constituents.

##  V SBE and H-IFE workflows results

The crystal structures of the 21 PFM candidates are retrieved from MPDS if the
selected polymorphs are available; otherwise, they are directly obtained from
literature crystallographic data. The crystal structures are then employed as
the only input required by the Surface Binding Energy (SBE) and the Hydrogen
Interstitial Formation Energy (H-IFE) workflows.

The results of the SBE workflow screening are reported in Tab. 5. $E_{b}$ is
the SBE of a surface atom, while $E_{p}(D)$ is threshold energy of physical
sputtering, in this case the minimum energy that an impinging deuterium atom
(D) must possess to remove the target atom from the surface (see Eq. 12). For
each material, a single value for $\overline{E_{b}}$ and for
$\overline{E_{p}}(D)$ is calculated by averaging the values of individual
elements using stoichiometric weighted average (SWA). The weights for this
averaging process are determined from the stoichiometric ratios of the
elements present in the material. The PFM candidates are then ranked in
descending order of $\overline{E_{p}}(D)$, which is the most crucial value in
assessing the material susceptibility to surface particle sputtering, thereby
influencing the erosion rate and the contamination of the plasma core. The
highest-performing materials are high-Z transition metals, in particular Re,
TaW, W, and Ta, followed by their carbide, nitride and boride forms. The
lowest $\overline{E_{p}}(D)$ values are obtained for low-Z transition metal
and ceramics. As shown in Eq. 12, the mass difference between plasma ions and
the target surface atom defines the amount of energy that can be transferred
during an elastic collision, which is maximum at equal masses. This implies
that, despite the chemistry at the surface, hence the SBE of a given atom, the
atomic mass governs its susceptibility to sputtering. Indeed, even if the
highest $\overline{E_{b}}$ are obtained for HfC, TaC, and graphite (C-194),
low $\overline{E_{p}}(D)$ are also shown, due to the low atomic mass of the
carbon atom. In particular, in the case of HfC, TaC, and WC, the transition
metals in the carbide surface display impressive $E_{b}$ values, higher than
in their pure metallic forms. It is important to emphasize that in this cases
preferential sputtering towards the carbon atom would occur, leaving a carbon-
depleted surface with modified thermal properties.

Table 5: Summary of the properties of the 21 PFM candidates, calculated
through the SBE workflow. The space group number ($SG$) is reported to define
the polymorph selected for each crystal. The surface orientation is chosen by
finding the lattice plane with minimum surface energy at a given termination,
min$\\{E^{hkl,\alpha}_{surf}\\}$. The surface binding energy, $E_{b}$, and
threshold energy of physical sputtering, $E_{p}(D)$, are calculated for each
inequivalent surface atoms of the 21 PFM candidates. The stoichiometric
weighted average values are calculated for each material, respectively
$\overline{E_{b}}$ and $\overline{E_{p}}(D)$. The PFM candidates are reported
in descending order of $\overline{E_{p}}(D)$.

Formula | $SG$ | min$\bm{\\{E^{hkl,\alpha}_{surf}\\}}$ | Miller index | element | $E_{b}$ | $E_{p}(D)$ | element | $E_{b}$ | $E_{p}(D)$ | element | $E_{b}$ | $E_{p}(D)$ | $\overline{E_{b}}$ |  $\overline{E_{p}}(D)$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\mathrm{eV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-2}$ | $\mathrm{eV}$ | $\mathrm{eV}$ | $\mathrm{eV}$ | $\mathrm{eV}$ | $\mathrm{eV}$ | $\mathrm{eV}$ | $\mathrm{eV}$ | $\mathrm{eV}$
Re | 194 | 0.19 |  $\langle$001$\rangle$ | Re | 11.72 | 276.81 |  |  |  |  |  |  | 11.72 | 276.81
TaW | 221 | 0.20 |  $\langle$110$\rangle$ | Ta | 12.34 | 283.44 | W | 10.90 | 254.28 |  |  |  | 11.62 | 268.86
W | 229 | 0.22 |  $\langle$110$\rangle$ | W | 11.16 | 260.38 |  |  |  |  |  |  | 11.16 | 260.38
Ta | 229 | 0.17 |  $\langle$110$\rangle$ | Ta | 10.92 | 250.67 |  |  |  |  |  |  | 10.92 | 250.67
HfC | 225 | 0.12 |  $\langle$100$\rangle$ | C | 9.83 | 19.98 | Hf | 16.68 | 378.03 |  |  |  | 13.26 | 199.00
TaC | 225 | 0.11 |  $\langle$100$\rangle$ | Ta | 15.12 | 347.25 | C | 9.39 | 19.10 |  |  |  | 12.26 | 183.17
W2C | 60 | 0.19 |  $\langle$111$\rangle$ | W | 11.10 | 258.83 | C | 9.80 | 19.93 |  |  |  | 10.67 | 179.20
WC | 187 | 0.25 |  $\langle$100$\rangle$ | C | 7.91 | 16.09 | W | 12.87 | 300.11 |  |  |  | 10.39 | 158.10
TaN | 187 | 0.19 |  $\langle$100$\rangle$ | Ta | 12.63 | 289.98 | N | 7.48 | 17.01 |  |  |  | 10.05 | 153.50
Mo | 229 | 0.19 |  $\langle$110$\rangle$ | Mo | 8.66 | 107.49 |  |  |  |  |  |  | 8.66 | 107.49
HfB2 | 191 | 0.21 |  $\langle$101$\rangle$ | Hf | 11.21 | 253.89 | B | 8.12 | 15.34 |  |  |  | 9.15 | 94.86
ZrB2 | 191 | 0.20 |  $\langle$101$\rangle$ | B | 8.03 | 15.17 | Zr | 11.13 | 131.63 |  |  |  | 9.06 | 53.99
TiN | 225 | 0.09 |  $\langle$100$\rangle$ | N | 8.90 | 20.25 | Ti | 11.32 | 73.06 |  |  |  | 10.11 | 46.65
Cr | 229 | 0.22 |  $\langle$110$\rangle$ | Cr | 6.21 | 43.26 |  |  |  |  |  |  | 6.21 | 43.26
C | 194 | 0.001 |  $\langle$001$\rangle$ | C | 16.71 | 33.96 |  |  |  |  |  |  | 16.71 | 33.96
TiB2 | 191 | 0.22 |  $\langle$101$\rangle$ | B | 8.16 | 15.40 | Ti | 10.31 | 66.51 |  |  |  | 8.87 | 32.44
VB2 | 191 | 0.15 |  $\langle$101$\rangle$ | B | 8.01 | 15.13 | V | 9.20 | 62.86 |  |  |  | 8.41 | 31.04
C | 227 | 0.35 |  $\langle$110$\rangle$ | C | 11.23 | 22.83 |  |  |  |  |  |  | 11.23 | 22.83
BC2N | 160 | 0.23 |  $\langle$001$\rangle$ | C | 12.58 | 25.57 | N | 7.68 | 17.46 | B | 10.33 | 19.52 | 10.79 | 22.03
BN | 216 | 0.04 |  $\langle$111$\rangle$ | B | 12.29 | 23.20 | N | 8.76 | 19.93 |  |  |  | 10.53 | 21.57
C3N4 | 176 | 0.04 |  $\langle$100$\rangle$ | N | 7.76 | 17.66 | C | 12.30 | 25.01 |  |  |  | 9.71 | 20.81

The results of the H-IFE workflow are presented in Tab. 6, where the hydrogen
interstitial formation energies $\Delta E^{f}(H_{i})$ of the 21 PFM candidates
are reported in descending order. Two values are reported, $\Delta
E_{a}^{f}(H_{i})$ and $\Delta E_{m}^{f}(H_{i})$; the first is calculated from
a reference hydrogen atom ($a$) as a reservoir, while the second from a
reference hydrogen molecule ($m$), which is reported for comparison with
literature. Their difference represent the $H_{2}$ binding energy, which in
our calculation (excluding ZPE) is equal to -2.312 $\mathrm{eV}$. A decrease
of $\Delta E^{f}(H_{i})$ leads to an exponential increase of the hydrogen
solubility in the crystal, as reported in Eq. 15; this is an undesirable
phenomenon for PFMs, which need to minimize tritium inventory for reasons of
safety and fuel efficiency. Only 4 materials exhibit a positive $\Delta
E_{a}^{f}(H_{i})$: BN, C graphite, BC2N, and C3N4, all of which possess high
atomic density. Hence, aside from chemical interactions, there appears to be a
correlation between smaller interstitial volumes, where atomic hydrogen must
fit, and higher $\Delta E^{f}(H_{i})$. The trend provided by $\Delta
E^{f}(H_{i})$ is only an initial indication that must be further refined by
selecting more realistic energy references for hydrogen reservoirs, by
including kinetic effects to define the hydrogen permeability inside the bulk
of the crystal, and by considering damage and diffusion along pores and
grains.

Table 6: Formation energy for hydrogen interstitial of the 21 PFM candidates. $\Delta E_{a}^{f}(H_{i})$ and $\Delta E_{m}^{f}(H_{i})$ are calculated using Eq. 16, hence not including ZPE effects, taking respectively the hydrogen atom and the hydrogen molecule as reference hydrogen reservoirs. The deviations from known literature values range around $\pm$ 0.1 $\mathrm{eV}$. This variance can be attributed to both differences in the size of the supercells (larger in this study compared to most of the cited works), the utilization of different functionals, and the choice of a different H2 chemical potential. Formula | $SG$ | $\Delta E_{a}^{f}(H_{i})$ |  $\Delta E_{m}^{f}(H_{i})$
---|---|---|---
($\mathrm{eV}$) | ($\mathrm{eV}$)
This work | This work | Ref.
BN | 216 | 2.81 | 5.12 | -
C | 227 | 2.71 | 5.02 | 4.9 [282]
BC2N | 160 | 1.45 | 3.76 | -
C3N4 | 176 | 0.27 | 2.58 | -
TiN | 225 | -0.13 | 2.18 | -
VB2 | 191 | -0.42 | 1.89 | -
TaC | 225 | -0.54 | 1.77 | -
HfB2 | 191 | -0.67 | 1.64 | -
C | 194 | -0.94 | 1.37 | -
ZrB2 | 191 | -1.01 | 1.30 | -
TiB2 | 191 | -1.09 | 1.22 | -
TaN | 187 | -1.50 | 0.81 | -
W | 229 | -1.53 | 0.78 | 0.87 [283]
HfC | 225 | -1.57 | 0.74 | -
WC | 187 | -1.77 | 0.54 | 0.7 [284]
Mo | 229 | -1.80 | 0.51 | 0.62 [285]
Cr | 229 | -1.81 | 0.51 | -
W2C | 60 | -1.88 | 0.43 | -
Re | 194 | -2.39 | -0.076 | -
TaW | 221 | -2.46 | -0.15 | -0.02 [286]
Ta | 229 | -2.82 | -0.51 | -

##  VI Comprehensive evaluation of PFM

The results obtained through the revisited thickness ranking,
$d^{ELMs}(T_{H})$, and the SBE and H-IFE workflows, $\overline{E_{p}}(D)$ and
$\Delta E_{a}^{f}(H_{i})$ respectively, are summarized in Tab. 7. The 21 PFM
candidates selected are ranked according to the $WF_{i}$ equation (see Sec.
II.1.2, Eq. 7).

The best PFM candidate according to this final ranking is carbon diamond
(C-227), which outperforms most of the alternatives in terms of thermal
properties and hydrogen interstitial formation energy. However, its low
resistance to physical sputtering and consequent hydrogen retention remains a
bottleneck for its usage as primary PFM. The energy threshold of physical
sputtering is even lower than that of graphite, which is notorious for its low
resistance to hydrogen sputtering. The application of diamond and diamond-like
carbon as coating layer might be more appealing, as it has been already
studied in the past [78, 79]; nevertheless, the current literature is missing
a realistic estimation of diamond erosion rates under the plasma flux of the
divertor region, a parameter that constrains the layer lifetime and defines
the tritium co-deposition rate and total tritium inventory. Diamond might also
be applied in composite materials, but further research on the effects of
degradation processes on the diamond thermal properties and intragranular
hydrogen solubility at divertor conditions is needed. BN shows similar
capabilities than carbon diamond, despite a lower $d^{ELMs}(T_{H})$; thus,
quantifying how BN thermal properties are affected by neutron-induced defects
at high temperatures is crucial to assess its applicability as primary PFM. BN
low sputtering resistance should be investigated to define its tritium co-
deposition rate, but less interaction with hydrogen isotopes is expected than
for pure carbon. Swelling and cracking of BN can be mostly avoided if
isotopically separated 11B is utilized to synthesize the crystal; however,
nitrogen transmutation into radioactive 14C can pose challenges in disposing
the material at the end of its lifecycle. As expected and reassuringly, W is
ranked at the top of the PFM candidates, showing the best thermal properties
among all metal-based materials and significant sputtering resistance, mainly
dictated by its high atomic mass. W shows a negative $\Delta
E_{a}^{f}(H_{i})$; however, the diffusion of hydrogen in the perfect crystal
is limited by the energy migration barrier. Indeed, it has been shown that
diffusion through vacancy hydrogen traps dominates over interstitial diffusion
below 1,500 $\mathrm{K}$ [287]. Both BC2N and C3N4 show properties that are
comparable with the high-performing BN. These materials are understudied for
their refractory properties due to difficulties with their synthesis, and
little is known about their interaction with the extreme environment of a
fusion reactor. Their low threshold energy of physical sputtering should be of
primary concern, given the presence of carbon atoms in the crystal structure.
Therefore, their response to plasma ion bombardment, coupled with tritium
retention assessment, is key in evaluating their potential applicability. TaC
and HfC have both good threshold energies of physical sputtering and
reasonable hydrogen solubility compared to the other PFM candidates. However,
a bottleneck is presented by their thermal properties: $d^{ELMs}(T_{H})$ is
slightly above 1 $\mathrm{mm}$ for TaC and even becomes negative for HfC.
Exploring their potential as primary PFMs should then involve their alloy
system 4TaC-1HfC. While availability of thermal properties for this alloy is
limited, it has an estimated $d^{ELMs}(T_{H})$ of about 5 $\mathrm{mm}$, as
shown in Sec. IV.4. Re and Ta exhibit comparable performance, a high
resistance towards physical sputtering due to their high atomic mass,
reasonable thermal properties and very high hydrogen solubility. Their use as
primary PFMs is limited by a $d^{ELMs}(T_{H})$ around 2 $\mathrm{mm}$, that
might severely restrict their lifetime under erosion conditions. Furthermore,
high tritium intake can be excepted given the strongly negative $\Delta
E_{a}^{f}(H_{i})$. TaW is the second-best candidate in terms of
$\overline{E_{p}}(D)$, surpassed only by Re. Its thermal properties seems to
limit its applicability as primary PFM; however, most of the literature
sources are based on experimental measures on alloys of various compositions,
and reliable data on the ordered crystal are still lacking. TaW shows some of
the highest potential in terms of hydrogen solubility. Graphite’s
$d^{ELMs}(T_{H})$ is second only to carbon diamond and its resistance to
physical sputtering is the highest among all carbon-based materials, with
hydrogen finding its lowest-energy interstitial site when adsorbed on top of a
carbon site between two hexagonal planes. The interlayer space can accommodate
interstitial hydrogen, in particular when the hexagonal structure is damaged.
Despite that, the primary processes governing hydrogen migration into the
graphite bulk involve intergranular diffusion on pore surfaces and penetration
through grain boundaries. In the latter section of the ranking, a limited set
of materials lacks outstanding performance, but their appealing thermal
characteristics render them worthy of attention. TaN, HfB2, WC, Mo, ZrB2, and
TiB2 have $d^{ELMs}(T_{H})$ values around 4 $\mathrm{mm}$, and
$\overline{E_{p}}(D)$ equal or higher than graphite. Evaluating their
properties under divertor conditions would be of great relevance, since most
of have not been tested as primary PFMs, except for WC and Mo, which are well-
studied materials considered as alternatives to pure W divertors. Finally, due
to the low $d^{ELMs}(T_{H})$, combined with mediocre sputtering and hydrogen
solubility resistances, W2C, TiN, VB2 and, Cr are unlikely to be employed as
primary PFMs, protective layers, or composite materials, given that superior
alternatives are available.

Table 7: The selected 21 PFM candidates are ranked accordingly to the $WF_{i}$ equation applied on the three properties calculated in this study: (i) maximum thickness to avoid melting during ELMs, $d^{ELMs}(T_{H})$; (ii) threshold energy of physical sputtering, $\overline{E_{p}}(D)$; (iii) hydrogen interstitial formation energy, $\Delta E_{a}^{f}(H_{i})$. The Pareto layer number $l$ of each material is also reported, showing consistency with the $WF_{i}$ ranking method. $WF_{i}$ | $l$ | Formula | $SG$ | $d^{ELMs}(T_{H})$ | $\overline{E_{p}}(D)$ |  $\Delta E_{a}^{f}(H_{i})$
---|---|---|---|---|---|---
(%) | ($\mathrm{mm}$) | ($\mathrm{eV}$) | ($\mathrm{eV}$)
86.3 | 1 |  C | 227 | 253.51 | 22.83 | 2.71
69.6 | 1 |  BN | 216 | 39.73 | 21.57 | 2.81
64.8 | 1 |  W | 229 | 10.53 | 260.38 | -1.53
63.1 | 2 |  BC2N | 160 | 50.40 | 22.03 | 1.45
60.3 | 1 |  TaC | 225 | 1.13 | 183.17 | -0.54
55.9 | 2 |  C3N4 | 176 | 63.78 | 20.81 | 0.27
54.4 | 1 |  Re | 194 | 2.32 | 276.81 | -2.39
51.7 | 2 |  TaW | 221 | 0.74 | 267.35 | -2.46
48.5 | 2 |  HfC | 225 | -0.82 | 199.00 | -1.57
45.8 | 2 |  Ta | 229 | 2.53 | 250.67 | -2.82
45.7 | 1 |  C | 194 | 66.17 | 33.96 | -0.94
39.8 | 1 |  TaN | 187 | 4.10 | 153.33 | -1.50
39.3 | 2 |  W2C | 60 | -0.001 | 177.19 | -1.88
37.6 | 1 |  HfB2 | 191 | 4.10 | 94.85 | -0.67
36.8 | 2 |  WC | 187 | 3.94 | 156.59 | -1.77
33.8 | 1 |  TiN | 225 | 0.58 | 46.66 | -0.13
25.6 | 1 |  VB2 | 191 | 1.15 | 31.04 | -0.42
23.7 | 2 |  Mo | 229 | 4.99 | 107.49 | -1.80
21.3 | 1 |  ZrB2 | 191 | 4.43 | 53.94 | -1.01
15.8 | 2 |  TiB2 | 191 | 5.38 | 32.44 | -1.09
6.3 | 3 |  Cr | 229 | 0.75 | 43.26 | -1.81

##  VII Conclusions and Outlook

The goal of finding divertor PFMs alternative to W has been assessed by
defining a target dataset of crystal structures and properties, the inorganic
crystal database MPDS [30, 31], and a set of physical constraints that every
candidate needs to satisfy. Three main loads affect the PFM of the divertor -
heat, plasma particles, and neutrons - and following for a few milliseconds
these interactions is computationally prohibitive even for a single material.
In this work we first focused on the heat balance equation that the PFM is
subject to during steady-state conditions and ELM events. Three experimentally
measurable thermal properties are used ($T_{melt}$, $k$, and $C_{p}$), and
when all of these are available on MPDS, the heat balance equation is used via
the thickness ranking procedure to discriminate between materials that melts
or not under divertor heat loads, discarding the latter. When the heat
equation could not be solved due to the lack of measure properties, materials
were ranked accordingly to the Pareto ranking method on the subset of the
available thermal proprieties. 71 PFM candidates have been selected at this
stage, and all of them have gone through a selective literature review of
their properties under relevant operational conditions. This has led to a
further screening of the candidates, which, following the use of the
comparative ranking criteria to test for false negatives (i.e., materials
erroneously discarded) results in a total of 21 potential PFMs. The crystal
structures of these materials are used as single input of the two first-
principles workflows developed to calculate the candidate surface biding
energy (SBE) and the hydrogen-interstitial formation energy (H-IFE). These two
properties are used as proxies for two undesired phenomena produced by the
plasma bombardment of exposed PFM: (i) the erosion of the surface due to
sputtering; (ii) the inventory of tritium inside the material due to hydrogen
isotopes diffusion. The materials selected via the screening process together
with the properties calculated in this study are summarized in Tab. 7. Well-
known and already proposed PFMs appear in this list, such as tungsten in its
pure metallic (W) and carbide forms (WC and W2C), diamond and graphite, boron
nitride, and other transition metals, such as Mo, Ta and Re. Other less
investigated refractory materials show interesting characteristic that are
worth further analysis, such as TaN in its $\theta$-phase, HfB2, ZrB2, and
TiB2. These do not excel in terms of threshold energy for physical sputtering
or hydrogen-interstitial formation energy, but they are still acceptable when
compared with those of graphite and tungsten. Especially, their good thermal
properties make these appealing as primary PFMs, in particular since the
recent possibility to synthesize borides with isotopically separated 11B,
which avoids the boron transmutation reaction. Novel PFMs are also reported,
in particular BC2N, C3N4, TaW in its ordered cubic structure, and 4TaC-1HfC,
which, for simplicity, is examined in this study in terms of its carbide
components. To the best of the authors’ knowledge, there is no literature
investigation of the thermal properties of these materials at high temperature
under plasma and neutron bombardment. Some limitation come from challenging
synthesizability, especially with the goal of obtaining homogeneous phases. In
this context, computational studies may offer valuable insights by predicting
their thermal properties across different conditions.

This work also highlights the importance of considering the engineering
constraints that a PFM in the divertor region needs to satisfy; the