# High Chern number phase in topological insulator multilayer structures: a
Dirac cone model study

Yi-Xiang Wang<EMAIL_ADDRESS>School of Science, Jiangnan
University, Wuxi 214122, China  School of Physics and Electronics, Hunan
University, Changsha 410082, China  Fuxiang Li<EMAIL_ADDRESS>School of
Physics and Electronics, Hunan University, Changsha 410082, China

###### Abstract

We use the Dirac cone model to explore the high Chern number $(C)$ phases that
are realized in the magnetic-doped topological insulator (TI) multilayer
structures by Zhao et al. [Nature 588, 419 (2020)]. The Chern number is
calculated by capturing the evolution of the phase boundaries with the
parameters and then the Chern number phase diagrams of the TI multilayer
structures are obtained. The high-$C$ behavior is attributed to the band
inversion of the renormalized Dirac cones, along with which the spin
polarization at the $\Gamma$ point will get increased. Moreover, another two
TI multilayer structures as well as the TI superlattice structures are also
studied.

Keywords: topological insulator multilayer, Chern number, phase diagram

PACS: 05.30.Fk, 73.20.-r, 81.30.Bx, 68.65.Ac

##  I Introduction

The quantum anomalous Hall (QAH) insulator, or called the Chern insulator, has
aroused the ongoing interests in condensed matter physics, where the
dissipationless chiral currents can flow in a zero magnetic field [1, 2, 3].
It was first observed in magnetic topological insulator (TI) thin films [4]
that are formed by an intricate interplay between the ferromagnetism due to
the magnetic-doping and the intrinsic spin-orbit coupling [5]. The magnetic-
doping breaks the time-reversal symmetry (TRS) in the system and can open a
small Dirac mass gap in the topological surface states [6]. In these systems,
the unavoidable inhomogeneity of the magnetic-doping would lead to the complex
magnetic orders [7, 8]. As a consequence, the QAH effect can only be observed
at extremely low temperature, which is about several tens of milikelvin, one
to two orders of magnitude below the Curie temperature of the magnetic dopants
[4, 8, 9, 10, 11, 12]. When the crystalline structure and magnetic order are
self-organized into a well-ordered topological superlattice, the QAH insulator
can be observed at much higher temperature [13, 14].

Besides the $C=1$ Chern insulator phase, a recent experimental work by Zhao et
al. demonstrated the existence of the high Chern number $C>1$ phase in the TI
multilayer structures [15]. The authors explained the high-$C$ behavior
through the interface Dirac state mechanism, where a nontrivial interface
state appears at the interface between the magnetic-doped and undoped TI
layers and when such a state is occupied, it will contribute
$\frac{e^{2}}{2h}$ to the total Hall conductance $\sigma_{xy}$. In our
previous work [16], we instead attributed the inverted bands in the high-$C$
phase to the two-dimensional (2D) subbands, which may not appear at the
interface, but are splitted by the thickness confinement of the total stacking
layers. In both theoretical analysis [15, 16], the bulk $\bm{k}\cdot\bm{p}$
model of the three-dimensional (3D) TIs [17, 18] were used to describe the TI
multilayer structures and the plane-wave expansion method was adopted to
satisfy the open boundary condition at both ends. The plane-wave expansion
method can correctly describe the top and bottom surface states, so it works
well for the TI thin film [18, 19, 20]. For the TI multilayers, they were
fabricated by the MBE method in the experiment and thus the linear Dirac cones
should be present on each surface of the TI layers. However, this fact may not
be effectively captured by the plane-wave expansion method, as the electronic
states were found to even penetrate into the deep TI layer interior [15, 16].

On the other hand, as the bulk states of the TIs are gapped, a highly
simplified Dirac cone model was initially proposed to explain the QAH behavior
in the magnetic TI thin film [5]. The model was later developed to describe
the magnetic-doped TI and ordinary insulator multilayers that can host the 3D
Weyl semimetal (WSM) phase [21, 22]. In the Dirac cone model, only the Dirac
cone degrees of freedom on each TI surface were retained, with the intralayer
and interlayer Dirac cone hoppings being included; while the bulk electronic
dynamics were completely ignored. This was demonstrated to be valid by
comparing the model parameters with the DFT band calculations [23].
Furthermore, the Dirac cone model was used to explore the dependence of the
QAH behavior on the film thickness, the magnetic configuration as well as the
stacking sequences of the TI layers [23].

Motivated by these progresses, here we try to investigate the high-$C$
behavior in the magnetic-doped TI multilayer structures from the point of view
of the Dirac cone model. We focus on the Chern number phase diagram, where the
Chern number is calculated by capturing the evolution of the phase boundaries
with the parameters [16, 19]. For comparison, we also consider another two TI
multilayer structures as well as the 3D TI superlattice structures. Our main
findings are as follows: (i) Within the Dirac cone model, the high-$C$
behavior is attributed to the band inversion of the renormalized Dirac cones
that is driven by the exchange splitting, along with which the spin
polarization at the $\Gamma$ point gets increased. In the highest-$C$ phase,
the occupied states are fully downspin polarized at the $\Gamma$ point. (ii)
To explain the observation that the Chern number decreases with the decreasing
middle magnetic-doped layer thickness [15], we assume that the related
intralayer Dirac cone hopping changes from negative to positive. Based on this
assumption, it is interesting to find that at certain magnetic-doping, the
Chern number can also get increased when the layer thickness decreases. (iii)
For another two TI multilayer structures, the band inversion cannot be
achieved for all negative-mass Dirac cones even when the magnetic exchange
splittings are very strong. (iv) For the 3D TI superlattice structures,
besides the high-$C$ phases, the WSM phases can appear and span broad
parameter regions. Our work can provide more insights on the high-$C$ behavior
realized in the TI multilayer structures, which may pave the way for the
future topological electronic devices.

##  II Model and method

Figure 1: (Color online) Schematics of the TI multilayer structure, including
the alternating magnetic-doped ($N_{1}=3$) and undoped ($N_{2}=2$) TI layers.
The magnetization in the magnetic-doped TI layer is along the $z$ direction
and the interlayer distance is denoted by $d$. The Dirac cone hoppings
$\Delta_{S}$ and $\Delta_{D}$, and the exchange splittings $J_{S}$ and $J_{D}$
are indicated by the arrows. Without hoppings and exchange splittings, the
Dirac cones are localized on the top and the bottom surface of each TI layer.

The experimental configuration of the TI multilayers [15] is plotted in Fig.
1. It stacks along the $z$ direction and includes the alternating magnetic-
doped and undoped TI layers, with the layer number in the plot being $N_{1}=3$
and $N_{2}=2$, respectively, so the undoped TI layer is sandwiched between the
doped TI layers.

We use the Dirac cone model to describe the dynamics of the TI multilayer
structures. In the model, the linear Dirac cones are assumed to be present on
both the top and bottom surfaces of the magnetic-doped TI layer as well as the
undoped TI layer. The Hamiltonian is written as $(\hbar=1)$ [21, 22, 23]

$\displaystyle H(\bm{k})$
$\displaystyle=\sum_{i}\Big{[}(-1)^{i}v(-k_{x}\sigma_{y}+k_{y}\sigma_{x})+m_{i}\sigma_{z}\Big{]}c_{\bm{k}i}^{\dagger}c_{\bm{k}i}$
$\displaystyle+\sum_{i\in\text{odd},j}(\Delta_{S}\delta_{i,j-1}+\Delta_{D}\delta_{i+1,j}+\Delta_{D}\delta_{i-1,j})c_{\bm{k}i}^{\dagger}c_{\bm{k}j}.$
(1)

Here $c_{\bm{k}i}$ $(c_{\bm{k}i}^{\dagger})$ annihilates (creates) a surface
Dirac cone with the wavevector $\bm{k}=(k_{x},k_{y})$ and the surface index
$i$. For $i$, we use the odd/even number to represent the top/bottom surface
of each layer, respectively. $v$ is the velocity of the Dirac cone, $\sigma$
denote the Pauli matrice that act on the spin space, $\Delta_{S}$ is the
intralayer Dirac cone hopping, and $\Delta_{D}$ the interlayer Dirac cone
hopping across the van der Waals gap. When the mass term $m_{i}$ is absent,
the system owns the TRS as ${\cal T}H(-\bm{k}){\cal T}^{-1}=H(\bm{k})$, with
the time-reversal operator ${\cal T}=i\sigma_{y}K$ and $K$ being the complex
conjugate operator. The mass term $m_{i}$ of each Dirac cone that is induced
by the exchange interactions with the nearest local magnetic moments can break
the TRS of the system,

$\displaystyle m_{i}=\sum_{\alpha}J_{i\alpha}M_{\alpha},$ (2)

where $\alpha$ is the layer index. We assume that the magnetization is along
the $z$ direction, so that $M_{\alpha}=1$ in the magnetic-doped layer, and
$M_{\alpha}=0$ in the undoped layer. As shown in Fig. 1, each surface Dirac
cone has one nearby magnetic-doped TI layer, which may be the same layer with
the exchange splitting $J_{i\alpha}=J_{S}$ or the adjacent layer with the
exchange splitting $J_{i\alpha}=J_{D}$. The Dirac cones localized on different
surfaces can be coupled by the Dirac cone hoppings and exchange splittings and
thus get renormalized. Note that the hoppings beyond the nearest-neighbor
Dirac cones are neglected because they are minor [23] and will not change the
main conclusions in this paper. In our calculations, we take $-\Delta_{S}=1$
as the unit of energy and set the parameter $\delta=\frac{J_{D}}{J_{S}}$.

To calculate the Chern number in the TI multilayer structures, we follow the
previous works [16, 19] and try to judge the Chern number by exactly capturing
the evolution of the phase boundaries with the parameters, as the Chern number
change is closely connected to the gap closings. To help the Chern number
judgment, two limiting cases need to be considered: (i) $J_{S}=0$ of no
magnetic-doping case, where the TRS is preserved and the Chern number must be
vanishing, with the renormalized Dirac cones always appearing in pairs with
opposite masses; (ii) $\Delta_{D}=0$ of the limiting isolated-layer case,
where the mass of each renormalized Dirac cone $M_{i}$ can be calculated
analytically and then the Chern number is obtained. As the Dirac cones on all
surfaces own the positive chirality $\chi_{i}=1$ and each Dirac cone can
contribute a component $C_{i}$ to the Chern number, we have

$\displaystyle C=\sum_{i}C_{i},\quad\text{with}\quad
C_{i}=\frac{1}{2}\text{sgn}(M_{i}).$ (3)

Compared with the common Chern number calculations by using the Kubo’s formula
[1] or the Fukui’s algorithm [24], where the exact diagonalization is needed
to obtain the eigenenergy and eigenstate for each wavevector in the Brillouin
zone (BZ), the above method requires much less computational resources and
time, and is expected to be extended to the antiferromagnetic configuration as
well as other topological systems.

Figure 2: (Color online) The energy bands $\varepsilon\sim vk_{x}$ when
setting $k_{y}=0$, and the wavefunction distributions at the $\Gamma$ point on
each Dirac cone in the magnetic-doped TI thin film. We choose $J_{S}=0.4$ in
(a)$-$(b), and $J_{S}=1.6$ in (c)$-$(d). In the wavefunction distributions,
the blue and orange bars denote the upspin and downspin contributions,
respectively and the bar heights are proportional to the weights of the
distributions. For clarity, the neighboring distributions are shifted
vertically by 0.8.

##  III magnetic-doped TI thin film

First we consider the magnetic-doped TI thin film. Within the Dirac cone
model, a pair of Dirac cones with the same chirality are localized on the top
and bottom TI surfaces and there exist the intralayer hoppings between the two
Dirac cones. In the four-orbit basis
$\begin{pmatrix}|1\uparrow\rangle&|1\downarrow\rangle&|2\uparrow\rangle&|2\downarrow\rangle\end{pmatrix}^{T}$,
where $1/2$ denotes the top/bottom surface state and $\uparrow/\downarrow$ the
upspin/downspin state, respectively, the Hamiltonian is [5]

$\displaystyle H(\bm{k})=\begin{pmatrix}J_{S}&-ivk_{-}&\Delta_{S}&0\\\
ivk_{+}&-J_{S}&0&\Delta_{S}\\\ \Delta_{S}&0&J_{S}&ivk_{-}\\\
0&\Delta_{S}&-ivk_{+}&-J_{S}\end{pmatrix},$ (4)

here $k_{\pm}=k_{x}\pm ik_{y}$.

To see the Dirac cone renormalization, we use the unitary matrix [5, 25]

$\displaystyle U=\frac{1}{\sqrt{2}}\begin{pmatrix}1&0&1&0\\\ 0&1&0&-1\\\
1&0&-1&0\\\ 0&-1&0&-1\end{pmatrix}.$ (5)

Then the Hamiltonian is transformed as

$\displaystyle H_{T}(\bm{k})=UH(\bm{k})U^{-1}$
$\displaystyle=\begin{pmatrix}J_{S}-|\Delta_{S}|&-ivk_{-}&0&0\\\
ivk_{+}&-J_{S}+|\Delta_{S}|&0&0\\\ 0&0&J_{S}+|\Delta_{S}|&ivk_{-}\\\
0&0&-ivk_{+}&-J_{S}-|\Delta_{S}|\end{pmatrix},$ (6)

which becomes block diagonalized. It shows that the chirality of the
renormalized Dirac cones in the upper/lower block remains the same as in Eq.
(4), but the Dirac cone mass in the upper/lower block will get
reduced/increased by the intralayer Dirac cone hopping. So the Chern number of
the system is obtained as

$\displaystyle
C=\frac{1}{2}\text{sgn}\big{(}J_{S}-|\Delta_{S}|\big{)}+\frac{1}{2}\text{sgn}\big{(}J_{S}+|\Delta_{S}|\big{)}.$
(7)

We can see that (i) when $J_{S}<|\Delta_{S}|$, the Dirac cone mass is negative
in the upper block and positive in the lower block, and thus the Chern number
$C=0$; (ii) when $J_{S}>|\Delta_{S}|$, the Dirac cone mass in the upper block
becomes positive, thus the band inversion occurs and $C=1$. So in a magnetic-
doped TI thin film, the nontrivial $C=1$ phase can be observed only when the
magnetic-doping reaches a certain ratio, which is consistent with the
experimental observations in the Cr-doped or V-doped (Bi,Sb)2Te3 thin film [4,
11, 10, 9, 12]. Note that the Dirac cone model cannot support the high-$C$
phase in the TI thin film, which, although predicted in theory [26, 19, 16],
but to our knowledge, has not yet been reported in experiment.

Figure 3: (Color online) Phase diagrams of the TI multilayer structures in the
parameter space $(J_{S},\Delta_{D})$, with the Chern number $C$ and the spin
polarization $\eta$ at the $\Gamma$ point being labeled as $(C,\eta)$. The
contour scale represents the magnitude of log${}_{10}\Delta$, where $\Delta$
is the energy gap of the lowest bands, and the bright blue lines denote that
the gap is closed. We choose the parameter $\delta=0.6$, the layer number
$N_{2}=1-5$ in (a) to (e) and $N_{1}=N_{2}+1$.

In Fig. 2, we plot the energy bands and the wavefunction distributions at the
$\Gamma$ point in the $C=0$ and $C=1$ phases. For the Hamiltonian including
only the linear $k$ terms, the normal and inverted bands are identical [Figs.
2(a) and (c)], so there is no signature of the band inversions in the band
structures. If the second order $k$ term $(\propto k^{2})$ is included, the
band inversion can be found as an additional energy extremum at the nonzero
wave vector, whereas a single extremum at the $\Gamma$ point remains in the
normal band [5]. On the other hand, the band inversions can be manifested in
the wavefunction distributions at the $\Gamma$ point. In the $C=0$ phase, the
orbital contributions to the first valence band (VB) are
$\frac{1}{\sqrt{2}}(-|1\uparrow\rangle-|2\uparrow\rangle)$, and those to the
first conduction band (CB) are
$\frac{1}{\sqrt{2}}(-|1\downarrow\rangle+|2\downarrow\rangle)$ [Fig. 2(b)],
while in the $C=1$ phase, the contributions to the first VB and first CB are
interchanged [Fig. 2(d)], meaning that all VBs at the $\Gamma$ point are
downspin polarized.

##  IV TI multilayer structures

Next, for the TI multilayer structures, the Chern number phase diagrams are
plotted in the parametric space $(J_{S},\Delta_{D})$ in Fig. 3, with the layer
number $N_{2}$ increasing from 1 to 5 in (a) to (e) and $N_{1}=N_{2}+1$.
Experimentally, tuning the parameters $J_{S}$ and $\Delta_{D}$ is quite
feasible, as the exchange splitting $J_{S}$ is expected to increase with the
magnetic-doping or increase with the decreasing temperature, and the
interlayer hopping $\Delta_{D}$ is expected to increase when the van der Waals
gap is narrowed by the external pressure [23, 27]. In the phase diagrams, the
contour scale represents the magnitude of log${}_{10}\Delta$, where $\Delta$
is the energy gap of the lowest bands and the bright blue lines denote that
the gap is closed.

Fig. 3 shows that the Chern number $C$ can increase from zero up to its
highest value $C=N_{1}+N_{2}$. When $J_{S}=0$ of no magnetic-doping case, the
Chern number always vanishes due to the presence of the TRS. On the other
hand, when $\Delta_{D}=0$ of the isolated-layer case, similar to Eq. (7), the
Chern number of the TI multilayer structures is obtained as

$\displaystyle C=$
$\displaystyle\frac{N_{1}}{2}\Big{[}\text{sgn}\big{(}J_{S}-|\Delta_{S}|\big{)}+\text{sgn}\big{(}J_{S}+|\Delta_{S}|\big{)}\Big{]}$
$\displaystyle+\frac{N_{2}}{2}\Big{[}\text{sgn}\big{(}J_{D}-|\Delta_{S}|\big{)}+\text{sgn}\big{(}J_{D}+|\Delta_{S}|\big{)}\Big{]}.$
(8)

We can see that when the exchange splitting $J_{S}$ increases, the negative
mass $J_{S}-|\Delta_{S}|$ of the Dirac cone in the magnetic-doped TI layer
becomes positive at the critical point $J_{S}=|\Delta_{S}|$, which is
$N_{1}$-fold degenerate; whereas the negative mass $J_{D}-|\Delta_{S}|$ of the
Dirac cone in the magnetic-undoped TI layer becomes positive at the critical
point $J_{D}=|\Delta_{S}|$, which is $N_{2}$-fold degenerate. So the Chern
number evolutions with $J_{S}$ are as follows: (i) when $J_{S}<|\Delta_{S}|$,
the Chern number vanishes $C=0$; (ii) when
$|\Delta_{S}|<J_{S}<\frac{|\Delta_{S}|}{\delta}$, $C=N_{1}$; (iii) when
$J_{S}>\frac{|\Delta_{S}|}{\delta}$, the Chern number increases to its highest
value $C=N_{1}+N_{2}$, where the band inversions occur for all negative-mass
Dirac cones.

When the interlayer hopping $\Delta_{D}$ is nonvanishing, the Dirac cones in
different TI layers will get coupled and further renormalized, which thus
breaks the Dirac cone degeneracy. The renormalized Dirac cones may not be
localized in one TI layer, but extended to all layers [see Figs. 4(a) and
(b)]. As a result, with the increasing of $J_{S}$, the band inversion occurs
successively for the renormalized Dirac cones. Each time a phase boundary is
crossed, the band inversion occurs for one negative-mass Dirac cone and the
Chern number will be changed by one. When $\Delta_{D}$ increases, we can see
that all phase boundaries originating from $J_{S}=\frac{|\Delta_{S}|}{\delta}$
move to larger $J_{S}$; on the other hand, for the phase boundaries from
$J_{S}=|\Delta_{S}|$, only one will move to $J_{S}=0$, meaning that the region
spanned by the $C=0$ phase is shrinking so that a minor magnetic doping can
drive the transition from the $C=0$ phase to $C=1$, while the remaining
$N_{1}-1$ ones all move to larger $J_{S}$.

In the phase diagrams, if the magnetic-doping is small, the value of the
critical point $J_{S}=\frac{|\Delta_{S}|}{\delta}$ at $\Delta_{D}=0$ will
become larger; while if the magnetic-ordered septuple layer is formed, such as
in MnBi2Te4 [28], where the exchange splitting is greatly strengthened,
$\delta$ is expected be close to 1 and thus the critical point is close to
$J_{S}=|\Delta_{S}|$. In the experiment [15], the thicknesses of the magnetic-
doped and undoped layers are taken as $d=3$ nm and $d=4$ nm, respectively. The
different TI layer thicknesses and properties can lead to the unequal
$\Delta_{S}$ in the different layers, which, however, will not change the
basic physics and the structure of the phase diagrams (see Appendix).

Next, we study the wavefunction distribution and the spin polarization in the
TI multilayer structures. The spin polarization $\eta_{\bm{k}}$ is defined as

$\displaystyle\eta_{\bm{k}}=\frac{\sum_{n\in\text{occ}}\big{(}|\psi_{n,\downarrow}(\bm{k})|^{2}-|\psi_{n,\uparrow}(\bm{k})|^{2}\big{)}}{\sum_{n\in\text{occ}}\big{(}|\psi_{n,\downarrow}(\bm{k})|^{2}+|\psi_{n,\uparrow}(\bm{k})|^{2}\big{)}},$
(9)

where $\psi_{n,\uparrow/\downarrow}$ is the upspin/downspin component of the
wavefunction with the band index $n$, and $n$ is to be summed over all
occupied states.

Because the band inversion is accompanied by the change of the orbital
contribution at the $\Gamma$ point from upspin to downspin, so $\eta_{\Gamma}$
is directly related to the Chern number,

$\displaystyle\eta_{\Gamma}=\frac{C}{N_{1}+N_{2}}.$ (10)

The above relation can also be demonstrated by the numerical results in Fig.
3. For example, in Fig. 4(a) and (b), with the parameters chosen as the points
labeled by the crosses in Fig. 3(b), we plot the wavefunction distributions of
all VBs at the $\Gamma$ point on each Dirac cone. For Chern number phase
changing from $C=1$ to $C=2$, we observe that besides some band movements, the
orbital contribution to the first VB changes from upspin to downspin, and thus
the spin polarization $\eta_{\Gamma}$ will increase by
$\frac{1}{N_{1}+N_{2}}$. At the highest Chern number phase $C=N_{1}+N_{2}$,
the spin polarization can reach its saturation value, $\eta_{\Gamma}=1$.

Figure 4: (Color online) (a)-(b) The wavefunction distributions of all VBs at
the $\Gamma$ point on each Dirac cone in the TI multilayer structures. The
blue and orange bars denote the upspin and downspin contributions,
respectively, and the bar heights are proportional to the weights of the
distributions. For clarity, the neighboring distributions are shifted
vertically by 0.3. (c) The spin polarization versus $vk_{x}$, where $J_{S}$
increases from 0.12 to 3.6 and the corresponding Chern number is labeled in
the right. We choose $N_{2}=2$ and $N_{1}=N_{2}+1$, $\Delta_{D}=1$,
$J_{S}=0.84$ in (a), $J_{S}=1.08$ in (b), and $k_{y}=0$ in (c).

We also investigate the spin polarization at nonzero wave vectors. In Fig.
4(c), when setting $k_{y}=0$, $\eta$ is plotted as a function of $vk_{x}$,
where the different lines correspond to $J_{S}$ increasing from $0.12$ to
$3.6$. We can see that although $\eta_{\bm{k}}$ may increase with $|k_{x}|$
around the $\Gamma$ point, it decreases when $|k_{x}|$ is large enough. On the
other hand, $\eta_{\bm{k}}$ always increases with $J_{S}$. So at the
asymptotically strong Zeeman splitting, the system may be driven into the
ferromagnetic Chern insulator phase [29]. In the experiment, the spin
polarization can be detected by the spin-polarized angle-resolved
photoemission spectroscopy [30].

Figure 5: (Color online) Phase diagrams of the TI multilayer structures in the
parameter space of $(J_{S},\Delta_{Sm})$, with the Chern number $C$ being
labeled. The contour scale represents the magnitude of log${}_{10}\Delta$,
where $\Delta$ is the energy gap of the lowest bands. The bright blue lines
denote that the gap is closed. We choose the parameter $\delta=0.6$, the layer
number $N_{2}=2$, $N_{1}=N_{2}+1$, and $\Delta_{D}=0.5$ in (a) and
$\Delta_{D}=1.5$ in (b).

The recent experiment that implemented on the TI multilayers
[(Bi,Sb)2-xCrxTe${}_{3}-$(Bi,Sb)2Te3]${}_{N_{2}}-$(Bi,Sb)2-xCrxTe3 reported
the observations of the high-$C$ phase with $C=N_{2}$ when the magnetic-doping
ratio is $x=0.24$ [15]. According to the above analysis, this suggests that at
such a ratio, the band inversion only occurs in $N_{2}$ negative-mass Dirac
cones. When $\Delta_{D}=0$, all these Dirac cones are localized in the
magnetic-doped TI layers. In their work, the authors also reported that for
the TI multilayers $N_{1}=3$ and $N_{2}=2$, the Chern number can change from
$C=2$ to $C=1$ when the middle magnetic-doped layer thickness was tuned to
decrease from $d=4$ nm to zero [15]. They suggested that when the magnetic-
doped layer thickness decreases, a pair of nontrivial interface states will
disappear, so the Chern number is reduced by one. Here in the Dirac cone
model, this observation could be instead explained by assuming that the
intralayer hopping within the middle magnetic-doped TI layer $\Delta_{Sm}$
will change from negative to positive when the layer thickness decreases. This
means that when the layer is thick enough, the Dirac cones attract each other
so that $\Delta_{Sm}$ is negative, while when the layer thickness decreases
and becomes small, the Dirac cones become less attractive and even repel each
other when a critical thickness is acrossed, so $\Delta_{Sm}$ becomes
positive. In our calculations, the chosen model parameters $\Delta_{S}=-1$ and
$\Delta_{D}>0$ are consistent with the above assumptions. Moreover, in Ref.
[23], by comparing the model parameters with the DFT band calculations, the
authors obtained the opposite $\Delta_{S}$ and $\Delta_{D}$, which also
supports our assumptions.

To see the effect of $\Delta_{Sm}$ on the Chern number modulation, in Fig. 5,
the phase diagrams of the TI multilayer structures are plotted in the
parameter space $(J_{S},\Delta_{Sm})$, with $\Delta_{D}=0.5$ in (a) and
$\Delta_{D}=1.5$ in (b). We can see that in Fig. 5(a), there is no $C=2$ to
$C=1$ phase transition along the vertical direction. However, in Fig. 5(b),
the $C=2$ to $C=1$ phase transition can be captured along the left solid
arrow, which agrees with the experiment [15]. More interestingly, in Fig.
5(b), when the magnetic-doping increases and the exchange splitting lies
within the range $1.51<J_{S}<1.84$, there can exist the $C=2$ to $C=3$ phase
transition, e.g., along the right dashed arrow. Therefore, as long as the
layer thickness is not vanishing and the Dirac cones are still present on the
surfaces, the Chern number may get increased. The conclusion is quite
different from the interface Dirac state mechanism, but also roots in the band
inversion of the Dirac cone.

##  V other TI multilayer structures

To make comparisons, we study another two TI multilayer structures, with the
Chern number phase diagrams plotted in Fig. 6. In Fig. 6(a), we choose the
layer number $N_{1}=2$ and $N_{2}=3$ so that the magnetic-doped TI layers are
sandwiched between the undoped ones, whereas in Fig. 6(b), we choose
$N_{1}=N_{2}=2$ and the number of the magnetic-doped TI layer is equivalent to
the undoped layer.

Figure 6: (Color online) Phase diagrams of two alternative TI multilayer
structures in the parameter space of $(J_{S},\Delta_{D})$, with the Chern
number $C$ and the spin polarization $\eta$ at the $\Gamma$ point being
labeled as $(C,\eta)$. The contour scale represents the magnitude of
log${}_{10}\Delta$, where $\Delta$ is the energy gap of the lowest bands and
the bright blue lines denote that the gap is closed. We choose the parameter
$\delta=0.6$, the layer number $N_{1}=2$, $N_{2}=N_{1}+1$ in (a) and
$N_{1}=N_{2}=2$ in (b).

In the isolated layer case of $\Delta_{D}=0$, for the topmost undoped TI
layer, the masses of the two renormalized Dirac cones are calculated as

$\displaystyle M_{t,\pm}^{o}=$
$\displaystyle\frac{1}{2}J_{D}\pm\frac{1}{2}\sqrt{J_{D}^{2}+4\Delta_{S}^{2}}.$
(11)

Clearly, $M_{t,\pm}^{o}$ is always positive/negative and there is no sign
change. So the contributions of the two Dirac cones to the Chern number are
always zero. Cases are the same for the bottommost undoped TI layer, if
exists. This means that the Chern number is contributed by the magnetic-doped
TI layers as well as the sandwiched undoped TI layers, but not the topmost or
bottommost undoped TI layer. This conclusion also holds when $\Delta_{D}$ is
nonvanishing. Therefore, as shown in Figs. 6(a) and (b), the highest Chern
number can only reach $C=2N_{1}-1$.

For the spin polarization $\eta$ at the $\Gamma$ point, the results are shown
in Fig. 6. We can see that the relation in Eq. (10) of the spin polarization
to the Chern number still holds. Comparing Figs. 6(a) and (b) with Fig. 3(a),
we find that although the nontrivial Chern number phases span the similar
regions, the spin polarizations are quite different. As analyzed above, even
when the Chern number takes its highest value, the spin polarization at the
$\Gamma$ point cannot be saturated in the two TI multilayer structures.

We note that in Fig. 6, the gaps are well opened in the Chern insulator
phases. Because the Fermi energy needs to be tuned in the gap, so the large
gap favors the experimental observations. This conclusion is quite different
from our previous work [16], where the gaps of the Chern insulator phases were
found to be too small to hinder the experimental observations.

##  VI TI superlattice structure

In addition, we study the TI superlattice structure with the periodic boundary
conditions in the $z$ direction. Note that in a unit cell, the topmost and
bottommost Dirac cones may have two nearby magnetic-doped TI layers. As more
and more magnetic-ordered TIs have been successfully fabricated in experiments
[28, 31, 32, 33], the TI superlattice considered here may also be realized in
the future.

Figure 7: (Color online) Phase diagrams of the TI superlattice structures in
the parameter space of $(J_{S},\Delta_{D})$, and the corresponding $k_{z}$
value for the energy gap. In (a) and (b), the Chern number $C$ is labeled and
the contour scale represents the magnitude of log${}_{10}\Delta$, where
$\Delta$ is the energy gap of the lowest bands. The bright blue regions denote
that the gap is closed, and the system lies in the gapless WSM phase. We
choose the parameter $\delta=0.6$, the layer number $N_{2}=1$ in (a) and (c),
$N_{2}=2$ in (b) and (d), and $N_{1}=N_{2}+1$.

The phase diagrams for the TI superlattice are plotted in Figs. 7(a) and (b),
where the unit cell includes the same TI layer number as in Figs. 3(a) and
(b), respectively. The bright blue regions denote that the energy gaps are
closed so that the system lies in the gapless WSM phase. Figs. 7(a) and (b)
show that the WSM phases can span large parameter regions and separate the
distinct Chern insulator phases. This is consistent with the previous studies
[21, 23], but is quite different from the TI multilayer structures, where no
WSM phase exists in Fig. 3. This is explained that the wave vector $k_{z}$
needs to take a specific value to close the gap so that the WSM phase can be
accommodated in the TI superlattice [21, 22], but $k_{z}$ cannot be modulated
in the TI multilayer structures. We find that the wave vector $k_{z}$ for the
gap can also be used to distinguish the Chern insulator and WSM phases, where
in the former, $k_{z}$ takes the fixed value $0$ or
$\frac{\pi}{(N_{1}+N_{2})d}$ that situates in the BZ center or boundary, while
in latter, $k_{z}$ gradually changes from $0$ to $\frac{\pi}{(N_{1}+N_{2})d}$.
Note that in the superlattice structure, due to the enlarged unit cell, the BZ
in the $z$ direction has been folded into
$[-\frac{\pi}{(N_{1}+N_{2})d},\frac{\pi}{(N_{1}+N_{2})d}]$.

When $J_{S}=0$ of no magnetic-doping case, the masses of the $2(N_{1}+N_{2})$
renormalized Dirac cones can be calculated analytically as

$\displaystyle
M_{\pm}^{s}=\pm\sqrt{\Delta_{S}^{2}+\Delta_{D}^{2}+2\Delta_{S}\Delta_{D}\text{cos}\big{(}k_{z}d+\frac{2\pi\alpha}{N_{1}+N_{2}}\big{)}},$
(12)

with the layer index $\alpha=1,2,\cdots,N_{1}+N_{2}$. The above equation shows
that the masses always appear in pairs. One can check that the gap can get
closed only when $M_{N_{1}+N_{2},\pm}^{s}$ vanish, with the conditions of
$k_{z}=0$ and $\Delta_{D}=|\Delta_{S}|$, as shown in Fig. 7.

On the other hand, when $\Delta_{D}=0$ of isolated-layer case, the masses of
the two renormalized Dirac cones on the topmost magnetic-doped TI layer in a
unit cell are obtained as

$\displaystyle
M_{t,\pm}^{s}=J_{S}+\frac{1}{2}J_{D}\pm\frac{1}{2}\sqrt{J_{D}^{2}+4\Delta_{S}^{2}},$
(13)

which are the same as the bottommost layer in a unit cell. One can infer from
Eq. (13) that the mass $M^{s}_{t,-}$ changes from negative to positive at the
critical point $J_{S}=\frac{|\Delta_{S}|}{\sqrt{1+\delta}}$ and is twofold
degenerate. The other critical points at $\Delta_{D}=0$ are the same as those
in TI multilayer structures, suggesting that the Chern number behaviors are
similar.

##  VII Discussions and Conclusions

In summary, we have studied the Chern number phase diagrams in the TI
multilayer structures by using the Dirac cone model, where the Chern number is
effectively calculated by capturing the evolutions of the phase boundaries
with the parameters. We find that the magnetic doping can drive the band
inversion of the Dirac cones and thus induce the high-$C$ phase transitions.
In the highest Chern number phase, the occupied states are downspin polarized
at the $\Gamma$ point.

Compared with the previous theoretical works that used the plane wave
expansion method on the $\bm{k}\cdot\bm{p}$ model [15, 16], our results show
the following differences: (i) the inverted bands are related to the
renormalized Dirac cones; (ii) the Chern number may get increased even when
the thickness of the middle magnetic-doped TI layer decreases; (iii) the gaps
are well opened in the Chern insulator phase in another two TI multilayer
structures. The differences are attributed to the fact that the surface states
of each TI layer in the multilayer structures are well captured by the Dirac
cone model, but may not be correctly described by the plane-wave expansion
method. We hope that these conclusions, especially (ii) and (iii), would be
demonstrated in the future TI multilayer experiments.

Besides the magnetic-doped TI multilayers, the magnetic-ordered ones, as well
as the factors that cannot be included in the DFT calculations, such as
disorder and magnetic field, can also be described by the Dirac cone model. We
expect that the Dirac cone model to be widely used to explain more emergent
phenomena in the TI multilayer structures, such as the axion insulator [34,
35].

##  VIII Acknowledgments

This work was supported by the National Natural Science Foundation of China
(Grant No. 11804122 and No. 11905054), the China Postdoctoral Science
Foundation (Grant No. 2021M690970), and the Fundamental Research Funds for the
Central Universities of China.

##  IX Appendix: unequal intralayer Dirac cone hoppings

Figure 8: (Color online) Phase diagrams of the TI multilayer structures in the
parameter space $(J_{S},\Delta_{D})$ in (a) and $(J_{S},|\Delta_{S2}|)$ in
(b), with the Chern number $C$ and the spin polarization $\eta$ at the
$\Gamma$ point being labeled as $(C,\eta)$. The contour scale represents the
magnitude of log${}_{10}\Delta$, where $\Delta$ is the energy gap of the
lowest bands, and the bright blue lines denote that the gap is closed. We
choose the parameter $\delta=0.6$, the layer number $N_{2}=2$,
$N_{1}=N_{2}+1$, $\Delta_{S1}=-1$, and in (a) $\Delta_{S2}=-0.8$, in (b)
$\Delta_{D}=0.9$.

Here we discuss the effect of the unequal intralayer Dirac cone hoppings on
the phase diagrams. We use $\Delta_{S1/2}$ to represent the intralayer Dirac
cone hopping in the magnetic-doped/updoped TI layer.

When $\Delta_{D}=0$, the Chern number in Eq. (8) becomes

$\displaystyle C=$
$\displaystyle\frac{N_{1}}{2}\Big{[}\text{sgn}\big{(}J_{S}-|\Delta_{S1}|\big{)}+\text{sgn}\big{(}J_{S}+|\Delta_{S1}|\big{)}\Big{]}$
$\displaystyle+\frac{N_{2}}{2}\Big{[}\text{sgn}\big{(}J_{D}-|\Delta_{S2}|\big{)}+\text{sgn}\big{(}J_{D}+|\Delta_{S2}|\big{)}\Big{]}.$
(14)

The above Chern number expression shows that the critical points are located
at $J_{S}=|\Delta_{S1}|$ and $J_{S}=\frac{|\Delta_{S2}|}{\delta}$, as shown in
Fig. 8(a). More importantly, the structure of the phase diagram in Fig. 8(a)
is similar with Fig. 3(b). When setting $\Delta_{D}=0.9$, in Fig. 8(b), we
plot the phase diagram in the parametric space $(J_{S},|\Delta_{S2}|)$. We
observe that as $|\Delta_{S2}|$ increases, the left three phase boundaries
show minor changes, while the right two phase boundaries will move to higher
$J_{S}$. So for the unequal $\Delta_{S1}$ and $\Delta_{S2}$, the qualitative
results of the phase diagram remain unchanged.

## References

  * [1] Thouless D J, Kohmoto M, Nightingale M P, and Nijs M den, 1982 Phys. Rev. Lett. 49 405
  * [2] Haldane F D M, 1988 Phys. Rev. Lett. 61 2015
  * [3] Qi X L, Hughes T L, and Zhang S C, 2008 Phys. Rev. B 78 195424
  * [4] Chang C Z, Zhang J, Feng X, Shen J, Zhang Z, Guo M, Li K, Ou Y, Wei P, Wang L L, Ji Z Q, Feng Y, Ji S, Chen X, Jia J, Dai X, Fang Z, Zhang S C, He K, Wang Y, Lu L, Ma X C, and Xue Q K, 2013 Science 340 167
  * [5] Yu R, Zhang W, Zhang H J, Zhang S C, Dai X, Fang Z, 2010 Science 329 61
  * [6] Chen Y L, Chu J H, Analytis J G, Liu Z K, Igarashi K, Kuo H H, Qi X L, Mo S K, Moore R G, Lu D H, Hashimoto M, Sasagawa T, Zhang S C, Fisher I R, Hussain Z, and Shen Z X, 2010 Science 329 659
  * [7] Zhang J, Zhao B, Zhou T, and Yang Z, 2016 Chin. Phys. B 25 117308
  * [8] Kou X, Pan L, Wang J, Fan Y, Choi E S, Lee W L, Nie T, Murata K, Shao Q, Zhang S C, and Wang K L, 2015 Nat. Commun. 6 8474
  * [9] Checkelsky J G, Yoshimi R, Tsukazaki A, Takahashi K S, Kozuka Y, Falson J, Kawasaki M, and Tokura Y, 2014 Nat. Phys. 10 731
  * [10] Feng Y, Feng X, Ou Y, Wang J, Liu C, Zhang L, Zhao D, Jiang T, Zhang S C, He K, Ma X, Xue Q K, and Wang Y, 2015 Phys. Rev. Lett. 115 126801
  * [11] Kou X, Guo S T, Fan Y, Pan L, Lang M, Jiang Y, Shao Q, Nie T, Murata K, Tang J, Wang Y, He L, Lee T K, Lee W L, and Wang K L, 2014 Phys. Rev. Lett. 113 137201
  * [12] Chang C Z, Zhao W, Kim D Y, Zhang H, Assaf B A, Heiman D, Zhang S C, Liu C, Chan M H W, and Moodera J S, 2015 Nat. Mater. 14 473
  * [13] Rienks E D L, Wimmer S, Sánchez-Barriga J, Caha O, Mandal P S, Ruzicka J, Ney A, Steiner H, Volobuev V V, Groiss H, Albu M, Kothleitner G, Michalicka J, Khan S A, Minár J, Ebert H, Bauer G, Freyse F, Varykhalov A, Rader O and Springholz G, 2019 Nature 576 423
  * [14] Deng H, Chen Z, Woloś A, Konczykowski M, Sobczak K, Sitnicka J, Fedorchenko I V, Borysiuk J, Heider T, Pluciński L, Park K, Georgescu A B, Cano J and Krusin-Elbaum L, 2021 Nat. Phys. 17 36
  * [15] Zhao Y F, Zhang R, Mei R, Zhou L, Yi H, Zhang Y Q, Yu J, Xiao R, Wang K, Samarth N, Chan M H W, Liu C X, and Chang C Z, 2020 Nature (London) 588 419
  * [16] Wang Y X and Li F, 2021 Phys. Rev. B 104 035202
  * [17] Zhang H, Liu C X, Qi X L, Dai X, Fang Z, and Zhang S C, 2009 Nat. Phys. 5 438
  * [18] Liu C X, Qi X L, Zhang H J, Dai X, Fang Z, and Zhang S C, 2010 Phys. Rev. B 82 045122
  * [19] Wang J, Lian B, Zhang H, Xu Y, and Zhang S C, 2013 Phys. Rev. Lett. 111 136801
  * [20] Zhou B, Lu H Z, Chu R L, Shen S Q, and Niu Q, 2008 Phys. Rev. Lett. 101 246807
  * [21] Burkov A A and Balents L, 2011 Phys. Rev. Lett. 107 127205
  * [22] Zyuzin A A, Wu S, and Burkov A A, 2012 Phys. Rev. B 85 165110
  * [23] Lei C, Chen S, and MacDonald A H, 2020 Proc. Natl. Acad. Sci. USA 117 27224
  * [24] Fukui T, Hatsugai Y, and Suzuki H, 2005 J. Phys. Soc. Jpn. 74 1674
  * [25] Wang Y X and Li F, 2018 EPL 123 37001\.
  * [26] Jiang H, Qiao Z, Liu H, and Niu Q, 2012 Phys. Rev. B 85 045445
  * [27] Guo W T, Huang L, Yang Y, Huang Z and Zhang J M, 2021 New J. Phys. 23 083030
  * [28] Otrokov M M, Klimovskikh I I, Bentmann H, Estyunin D, Zeugner A, Aliev Z S, Ga$\beta$ S, Wolter A U B, Koroleva A V, Shikin A M, Blanco-Rey M, and et al. 2019 Nature 576 416
  * [29] Wang W, Ou Y, Liu C, Wang Y, He K, Xue Q K, and Wu W, 2018 Nat. Phys. 10 1038
  * [30] Sobota J A, He Y, and Shen Z X, 2021 Rev. Mod. Phys. 93 025006
  * [31] Wu J, Liu F, Sasase M, Ienaga K, Obata Y, Yukawa R, Horiba K, Kumigashira H, Okuma S, Inoshita T, and Hosono H, 2019 Sci. Adv. 5 eaax9989
  * [32] Vidal R C, Zeugner A, Facio J I, Ray R, Haghighi M H, Wolter A U B, Bohorquez L T C, Caglieris F, Moser S, Figgemeier T, Peixoto T R F, and et al. 2019 Phys. Rev. X 9 041065
  * [33] Klimovskikh I I, Otrokov M M, Estyunin D, Eremeev S V, Filnov S O, Koroleva A, Shevchenko E, Voroshnin V, Rybkin A G, Rusinov I P, and et al, 2020 npj Quantum Mater. 5 54
  * [34] Mogi M, Kawamura M, Tsukazaki A, Yoshimi R, Takahashi K S, Kawasaki M, Tokura Y, 2017 Sci. Adv. 3 eaao1669
  * [35] Xiao D, Jiang J, Shin J H, Wang W, Wang F, Zhao Y F, Liu C, Wu W, Chan M H W, Samarth N, and Chang C Z, 2018 Phys. Rev. Lett. 120 056801


# High throughput quantitative metallography for complex microstructures using
deep learning: A case study in ultrahigh carbon steel

Brian L. DeCost1 Bo Lei2 Toby Francis2 Elizabeth A. Holm2

( 1Material Measurement Laboratory, National Institute of Standards and
Technology, Gaithersburg MD, 20899, USA
2Materials Science and Engineering, Carnegie Mellon University, Pittsburgh PA,
15213, USA

Accepted for publication in Microscopy and Microanalysis Aug. 21, 2018
dx.doi.org/10.1017/S1431927618015635 )

###### Abstract

We apply a deep convolutional neural network segmentation model to enable
novel automated microstructure segmentation applications for complex
microstructures typically evaluated manually and subjectively. We explore two
microstructure segmentation tasks in an openly-available ultrahigh carbon
steel microstructure dataset [1, 2]: segmenting cementite particles in the
spheroidized matrix, and segmenting larger fields of view featuring grain
boundary carbide, spheroidized particle matrix, particle-free grain boundary
denuded zone, and Widmanstätten cementite. We also demonstrate how to combine
these data-driven microstructure segmentation models to obtain empirical
cementite particle size and denuded zone width distributions from more complex
micrographs containing multiple microconstituents. The full annotated dataset
is available on materialsdata.nist.gov [3].

##  1 Introduction

Quantitative microstructure analysis is central to materials engineering and
design. Traditionally this entails careful measurements of volume fractions,
size distributions, and shape descriptors for familiar microstructural
features such as grains and second-phase particles. These quantities are
connected to theoretical and/or empirical models for materials properties,
e.g. grain boundary [4] or particle [5] strengthening mechanisms. Contemporary
microstructure segmentation methods rely on specialized image processing
pipelines that often require expert tuning for application to a particular
microstructure system. Furthermore, the microstructures accessible to
quantitative analysis are limited by the use of segmentation algorithms that
rely on low-level image features (intensity and connectivity constraints). In
this work, we apply deep learning methods for image segmentation to complex
microstructure data, with the goal of extending the reach of quantitative
analysis to microstructure systems that are currently evaluated subjectively
or through laborious manual annotation.

Since 2012, deep learning methods [6] have dominated many computer vision
applications111See [7] for a comprehensive introduction to deep learning
methods, including architectural and training choices., including object
recognition and detection, scene summarization, semantic segmentation, and
depth map prediction. The success of deep learning is often attributed to the
ability of convolutional neural networks (CNNs) to learn to effectively
represent the hierarchical structure of visual data, composing low-level image
features (edges, color gradients) into higher level features corresponding to
abstract qualities of the image subject (e.g. object parts). Recently,
materials scientists have begun exploring a limited set of applications of
contemporary computer vision techniques for flexible and generic
microstructure representation. [8] and [9] explore these techniques in the
context of microstructure classification. [10] and [11] use pretrained CNN
representations to study relationships between processing conditions and
microstructure via dimensionality-reduction and visualization techniques. [12]
use a CNN segmentation model to identify constituent phases in steel
microstructures.

In this report, we train a pixelwise CNN [13] to segment microstructures at a
high level of abstraction, and investigate the potential for this technique to
enable quantitative microstructure analyses that conventionally would require
a large amount of hands-on image processing. We evaluate the feasibility of
this approach on a subset of the openly available Utrahigh Carbon Steel (UHCS)
microstructure dataset [1, 14, 15]. CNNs can distinguish between the four
principal microconstituents in this heat-treated UHCS: proeutectoid cementite
network, fields of spheroidite particles, the ferritic matrix in the particle-
free denuded zone near the network, and Widmanstätten laths. We also train a
network to segment individual spheroidite particles, and briefly explore
automated microstructure metrology techniques enabled by this kind of powerful
segmentation model. Our training data and annotations for both microstructure
segmentation tasks will be publicly available through the NIST materials
resource registry [3].

Our primary contributions are:

  * •

Establishing two novel microstructure segmentation benchmark datasets

  * •

Connecting microstructure science to the deep semantic segmentation literature

  * •

Exploring novel means of expanding contemporary quantitative microstructure
measurement techniques to more complex structures

For microstructure scientists, CNN-based microstructure segmentation tools
require an initial investment in annotation and training, but can enable
longer-term or larger-scale research and characterization efforts. This trade-
off is particularly attractive for its potential to enable microstructure-
based material qualification by making it easier/cheaper to obtain statistical
data on high-level microstructure features known to mediate critical
engineering properties of materials (e.g. particle size distributions; denuded
zone widths, and particle coarsening kinetics). In industrial settings where
reliance on semi-automated segmentation techniques is common, the barrier to
entry is even lower because the training data has already been collected. CNN-
based microstructure segmentation tools also offer a path forward to high-
throughput microstructure quantification techniques for accelerated alloy
design and processing optimization, where acquisition and analysis of high-
quality microstructure data is often a limiting factor.

##  2 Methods

###  2.1 Segmentation model

Recently a variety of deep CNN architectures have been developed for dense
pixel-level tasks [16], such as semantic segmentation [17], edge detection,
depth map, and surface normal prediction [18]. Conceptually, a modern deep CNN
computes a highly nonlinear function through a layerwise composition of
convolution, activation, and pooling (i.e. downsampling) functions, the
parameters of which are learned from large annotated datasets by some variant
of stochastic gradient descent [6, 7]. Classification CNNs reduce an input
image to a single latent feature vector, where CNNs designed for pixel-level
tasks produce a latent representation for every pixel of the input image. This
is typically accomplished by upsampling the intermediate feature maps via a
fixed bilinear interpolation [19, 13] or a learned deconvolution operation
[20]. In the latter class of networks, popular architectures include SegNet
[17], Bayesian SegNet [21], U-Net [22] with heavy data augmentation, and
fully-convolutional DenseNets [23]. In particular, U-Net [22] was designed for
application to medical image segmentation tasks with small dataset sizes,
relying on strong data augmentation to achieve good performance.

####  2.1.1 PixelNet architecture

The PixelNet [13] architecture is illustrated schematically in Figure 1.
PixelNet applies bilinear interpolation to intermediate feature maps to form
hypercolumn features $h(x)=[conv_{1}(x),\;conv_{2}(x),\;\ldots\;conv_{5}(x)]$,
which represent each pixel in the input image with information drawn from
multiple scales. A non-linear classifier implemented as a multi-layer
perceptron (MLP, i.e. a traditional artificial neural network (ANN)) maps the
hypercolumn features to the corresponding pixel-level target. Instead of
computing dense high-dimensional feature maps at the input resolution as in
other popular pixel prediction networks, at training time PixelNet performs a
sparse upsampling to efficiently obtain hypercolumn features only for a small
sample of the input pixels.222Our tensorflow implementation of PixelNet is
available at https://github.com/bdecost/pixelnet This is attractive for
quickly training segmentation networks from scratch with small training sets
because it reduces the memory footprint during training and makes training a
non-linear predictor with high-dimensional latent representations feasible
[13].

The feature extraction portion of our PixelNet variant uses the VGG-16
architecture [24] used by the original PixelNet [13]; this architecture
consists of 13 convolution layers and two fully-connected layers
${1_{1},1_{2},2_{1},2_{2},3_{1},3_{2},3_{3},4_{1},4_{2},4_{3},5_{1},5_{2},5_{3},6,7}$.
The MLP layers in our PixelNet variant consist of 1024 neurons with rectified
linear (ReLU) activations [25] ($ReLU(y_{i})=\max(0,y_{i})$ followed by batch
normalization [26]. Following the original PixelNet implementation, our
hypercolumn features consist of the highest convolution feature map within
each block of the VGG architecture
({1${}_{\text{2}}$,2${}_{\text{2}}$,3${}_{\text{3}}$,4${}_{\text{3}}$,5${}_{\text{3}}$,7}),
converting layer $7$ to a $7\times 7$ convolution filter as in [20] and [13].
We apply batch normalization [26] to each VGG-16 feature map before upsampling
via bilinear interpolation, immediately after the ReLU activations.

Figure 1: The PixelNet [13] image segmentation model. A pixel in the input
image (left) is represented by the concatenation of its representations in
each convolution layer (white dots). A multilayer perceptron (MLP) classifier
is trained to associate the pixel representation with membership in a
microstructure constituent (right).

####  2.1.2 Training details

We initialize the feature extraction portion of our networks with a pre-
trained VGG-16 [24] network trained on the ImageNet [27] classification
dataset. We train the pixel classification layers from scratch, randomly
sampling initial weights from Gaussian distributions with zero mean and
standard deviation $\sigma=\sqrt{2/c}$ [28], where $c$ is the dimensionality
of the input to the layer. To prevent overfitting, we use a combination of
batch normalization [26], Dropout regularization [29], weight decay
regularization [30] , and data augmentation. We set the weight decay strength
to 0.0005 and apply Dropout regularization with a rate of 10% after the final
MLP layer. Training images are subjected to local histogram equalization to
mitigate differences in overall brightness across different samples and
datasets. The training input and label images are augmented with random
rotations in the range $\mathopen{[}0,2\pi\mathclose{)}$, horizontal and
vertical mirror symmetry, scaling in the range $\mathopen{[}1,2\mathclose{]}$,
and a ± 5% random intensity shift. Rotated versions of the training input and
label images are computed with mirror boundary conditions, with bilinear
interpolation for the input images and nearest-neighbor interpolation for the
(discrete) label images. We train the networks with the AdamW optimizer [31,
30] with the recommended default parameters. First we fix the parameters in
the feature extraction portion of the network and train the pixel
classification layers with an initial learning rate of 10${}^{\text{-3}}$ for
20 epochs (125 gradient updates). Each gradient update is computed from a
random sample of 2048 pixels each from 4 augmented training images. We then
fine-tune the entire CNN for 125 additional gradient updates using AdamW with
an initial learning rate of 10${}^{\text{-5}}$.

Our dataset has a heavy class imbalance (e.g. Widmanstätten cementite only
accounts for $\sim$ 3% of pixels), so we compare a model trained using the
standard cross-entropy classification loss with another trained using the
focal loss [32], which is designed for unbalanced datasets. The focal loss
extends the cross-entropy loss function $CrossEntropy(p_{t})=-\log(p_{t})$,
where

$p_{t}(p,y)=\begin{cases}p&\mathrm{if}y=1\\\ 1-p&\mathrm{if}y=0\end{cases}$
(1)

with ground truth $y$ and predicted class probability $p=P(y=1)$. The focal
loss adds a modulating factor $(1-p_{t})^{\gamma}$ to emphasize examples about
which the classifier is less confident during training, and a scaling
parameter $\alpha$ to account for class imbalance:

$FocalLoss(p_{t})=-\alpha_{t}(1-p_{t})^{\gamma}\log(p_{t})$ (2)

We follow the recommendation of [32] in setting the focusing parameter
$\gamma=2$ and setting the class imbalance parameters $\alpha_{t}$
proportionally to the inverse frequency of each class.

###  2.2 Dataset

The semantic microstructure segmentation dataset consists of 24 manually
annotated333We used the medical image annotation system MITK [33]. micrographs
from the open UHCS dataset [1, 2]; examples are shown in Figure 2 and in the
online supplemental materials. These $645\times 484$ pixel micrographs focus
on the characteristic features of heat-treated UHCS: the proeutectoid
cementite network and the associated denuded zone, and spheroidized and
Widmanstätten cementite. Multiple heat treatment conditions and magnifications
are represented in the semantic microstructure segmentation dataset.

The particle segmentation dataset consists of 24 micrographs collected at a
single magnification in support of the particle coarsening analysis reported
in [15]. Particle annotations were obtained through a partially-automated
edge-based segmentation workflow [15]. A thresholded blur smooths contrast in
the matrix surrounding particles before application of the Canny edge detector
[34]. The particle outlines are filled in, and spurious edges (e.g. at grain
boundaries) are removed by a 2px median filter. The final particle
segmentations are verified and retouched manually where the contrast is
insufficient for the Canny detector to identify particle edges. Particles
intersecting the edge of the image are removed from the annotations to reduce
bias in the estimated particle size distributions.

###  2.3 Performance evaluation

Because our set of annotated images is small (24 annotated micrographs total),
we use cross-validation to estimate the generalization performance of the
PixelNet architecture on our two microstructure segmentation tasks. We use a
6-fold cross-validation scheme [35]: each dataset is split into six validation
sets of four micrographs each, and six PixelNet models are trained on each of
the complementary training sets. The quantitative performance metrics reported
in Tables 1 and 4 are averages over each validation image in the 6 validation
sets; uncertainties are standard errors computed over the six validation
images [35].

We report several standard evaluation metrics for semantic segmentation tasks:
pixel accuracy (AC), precision, recall, and region intersection over union
(IU) for individual microconstituents. For each of these metrics, a higher
score indicates better performance.

Precision is the fraction of instances predicted to have class $c$ that are
correct:

$Precision(c)=\frac{\sum_{i}\hat{y_{i}}=c\textrm{ and
}y_{i}=c}{\sum_{i}\hat{y_{i}}=c}$ (3)

where $\hat{y}_{i}$ indicates the predicted class label for each pixel $i$,
and $y_{i}$ indicates the corresponding ground truth class label.
Equivalently, precision is the ratio of true positives to total (true and
false) positives, which decreases when the model overpredicts the number of
member pixels in a class.

Recall is the fraction of instances with ground truth class $c$ that are
predicted to have class $c$:

$Recall(c)=\frac{\sum_{i}\hat{y_{i}}=c\textrm{ and }y_{i}=c}{\sum_{i}y_{i}=c}$
(4)

Equivalently, recall is the ratio of true positives to the total number of
pixels in a class, which decreases when the method underpredicts the member
pixels in a class. Since the overall accuracy is defined as the number of true
positives divided by the total number of pixels, it is straightforward to show
that the classwise average recall or precision equals the overall accuracy.

The intersection over union metric $IU(c)$ for class $c$ (also referred to as
the Jaccard metric) is the ratio of correctly predicted pixels of class $c$
(true positives) to the union of pixels with either ground truth or predicted
class $c$ (true and false positives plus false negatives):

$IU(c)=\frac{\sum_{i}\hat{y_{i}}=c\textrm{ and
}y_{i}=c}{\sum_{i}\hat{y_{i}}=c\textrm{ or }y_{i}=c}$ (5)

For the spheroidite particle segmentation task, we also report performance
metrics comparing particle size distributions (PSDs) obtained from the model
predictions with those obtained from the ground truth annotations (as reported
in [15]). We use the two-sample Kolmogorov-Smirnov (KS) test [36] to compare
each pair of predicted and ground truth PSDs. The KS score reported in Table 4
is the fraction of micrographs where the KS test indicates that the predicted
particle size distribution is consistent with the ground truth particle size
distribution (i.e. the fraction of micrographs where we fail to reject (at the
95% confidence level) the null hypothesis that the distributions are
equivalent).

###  2.4 Computing denuded zone widths

Given a microconstituent prediction map, we quantify the width of the denuded
zone by computing the minimum distance to the network phase for each pixel on
the matrix-particle interface. In practice, we compute a map of Euclidean
distance to the network phase, and select the measurements at the denuded zone
interface.

To obtain the denuded zone interface, we apply a series of image processing
techniques to clean up the microconstituent prediction map, so that only the
matrix predictions associated with the diffusion-limited denuded zone adjacent
to the proeutectoid cementite network remain. A morphological filling
operation removes any matrix pixels within the network. Matrix regions that
are not connected to the network are identified by application of a
morphological closing to matrix phase: any matrix segments that do not
intersect the network phase after the morphological operation are removed.
Finally, we remove any matrix predictions that are closer to a widmanstatten
region than to a network region, and subsequently remove the widmanstatten
regions. The region boundaries on the cleaned up label image (shown in Figure
5) include only the interface of the proeutectoid cementite network phase
(indicated in blue) and the diffuse interface of the denuded zone (indicated
in yellow).

##  3 Results and Discussion

###  3.1 Semantic microconstituent segmentation

Figure 2 shows microconstituent annotations (e-h) and predictions (i-p) for
the four validation set micrographs (a-d) in one cross-validation iteration;
results for all six validation sets are included in the online supplemental
materials. Microconstituent predictions using the focal loss function and the
cross-entropy loss function are compared in Figure 2 (i-l) and (m-p)
respectively. The predictions show reasonable correspondence with the
annotations, despite nontrivial differences in features, such as particle size
and appearance that arise from differences in heat treatment and
magnification. Intensity variations and polishing damage evident in the input
images have little impact on the predictive capability of the models. One
notable exception is the cluster of spurious network predictions associated
with the damaged areas in the lower left of Figure 2 c. Both models do a good
job respecting the edges of the network carbide phase, with a few exceptions
where the network is very fine or the contrast between network carbide and
metal matrix is poor (see supplemental Figures S1.1 d and S1.5 d). Predicted
boundaries between spheroidite particles and the denuded zone have little
noise and tend to be smoother than in the annotations. The Widmanstätten
predictions show the highest amount of noise, especially where the
Widmanstätten lath are fine or are beginning to break up. The focal loss also
tends to surround Widmanstätten cementite with wider swaths of the metallic
matrix compared to the annotations. In addition to the low area fraction of
Widmanstätten cementite, one potential contributing factor for these failure
modes is labeling bias where the microstructure is ambiguous even to the human
expert. For example, some areas with a low density of spheroidite particles
are labeled by the model as metallic matrix where the annotation has made no
such distinction. This phenomenon is evident in the lower half of Figure 2 i,
where the model correctly identifies large patches of bare metal in the
neighborhood of some large grain boundary cementite particles (refer to
supplementary Figure S1.13 a for a more detail).

The cross-entropy model segmentation maps (Figure 2 (m-p)) tend to be more
consistent with the annotations for the majority microconstituents. However,
each model errs from the annotations in distinct ways. In general, the focal
loss model seems to emphasize constituent contiguity, while the cross-entropy
model tries to resolve fine features. For example, compared with the
annotations in Figure 2 f, the focal loss model (Figure 2 j) more liberally
identifies the very fine lath structures in the bottom right corner of the
frame, consolidating them into a single patch, while the cross-entropy model
(Figure 2 n) produces a noisier map that attempts to track finer-grained
details of the lath structure but loses some area fraction. One conclusion
from these results is that while both models give acceptable results, neither
is necessarily the best that can be achieved. For a given segmentation
problem, the user must select model parameters and loss functions to achieve
the desired quantitative and qualitative performance.

Figure 2: (a-d) Validation set micrographs, (e-h) microconstituent
annotations, (i-l) PixelNet predictions using the focal loss segmentation
model, and (m-p) PixelNet predictions using the standard crossentropy
classification loss. Microstructural constituents include proeutectoid grain
boundary cementite (light blue), ferritic matrix (dark blue), spheroidite
particles (yellow), and Widmanstätten cementite (green). Scale bars indicate
$10\mu m$.

Table 1 shows the average validation set performance with standard errors for
the focal loss model. The model obtains 86.5 $\pm$ 1.6 % overall accuracy (AC,
equivalent to the average of the classwise recall or precision) in reproducing
the pixel-level annotations. The model is consistently good at identifying
spheroidite and network regions. The less prevalent microconstituents (matrix
and Widmanstätten) are not as well captured, and show higher variation between
images. For these microconstituents, the recall score is better than the
precision score, meaning that the CNN tends to mistake other classes for
matrix and Widmanstätten more than it tends to miss genuine matrix and
Widmanstätten pixels. This effect is demonstrated on the fine Widmanstätten
lath in the lower right portion of Figure 2 j, where the model fills in the
fine spacing between Widmanstätten lath in its prediction. The low proportion
of Widmanstätten pixels in the dataset enhances this effect. In the case of
the matrix class, the difference in recall and precision scores is partly due
to the overprediction of metallic matrix in areas containing a low density of
spheroidite particles, as discussed in reference to Figure 2 i.

Table 1: Focal model semantic segmentation performance averaged over validation images. Uncertainties are standard errors calculated across validation images. | IU | precision | recall
---|---|---|---
matrix | 49.1 $\pm$ 3.4 | 60.3 $\pm$ 4.4 | 72.3 $\pm$ 3.7
network | 72.9 $\pm$ 5.3 | 85.5 $\pm$ 4.0 | 80.7 $\pm$ 5.9
spheroidite | 85.7 $\pm$ 1.8 | 95.1 $\pm$ 1.2 | 89.8 $\pm$ 1.7
widmanstatten | 42.7 $\pm$ 2.9 | 50.2 $\pm$ 3.6 | 73.5 $\pm$ 3.9
overall | 62.6 $\pm$ 2.5 | 86.5 $\pm$ 1.6 | 86.5 $\pm$ 1.6

In contrast, the spheroidite and network classes have slightly higher
precision compared with their recall scores. The standard error for the
network scores is large, and is therefore likely accounted for by the small
number of gross errors discussed in supplemental Figures S1.1 d and S1.5 d.
Finally, the small difference in precision and recall score for the
spheroidite class is likely also due to the overprediction of the metal matrix
in regions with low particle density.

Table 2 shows the same performance metrics for the cross-entropy model trained
with revised training hyperparameters. The crossentropy shows clearly superior
overall numerical performance, including a nearly ten point bump in overall IU
score. While most of the per-microconstituent scores are higher for the
crossentropy model, the recall score for Widmanstätten cementite is
consistently depressed due to underprediction. We also briefly experimented
with the popular U-Net [22] architecture, but found that this architecture
slightly underperformed compared with Pixelnet. These results indicate that
various CNN architectures and training schemes can achieve reasonable results,
so the user can select an approach based on desired outcomes. Because of its
excellent performance in the spheroidite segmentation task (reported in the
next section), we present results only from the Pixelnet model trained with
the focal loss function throughout the rest of this manuscript.

Table 2:  Crossentropy model semantic segmentation performance averaged over validation images. Uncertainties are standard errors calculated across validation images. | IU | precision | recall
---|---|---|---
matrix | 67.3 $\pm$ 6.2 | 81.8 $\pm$ 8.0 | 79.4 $\pm$ 4.0
network | 89.4 $\pm$ 3.6 | 96.1 $\pm$ 0.6 | 92.8 $\pm$ 3.5
spheroidite | 91.9 $\pm$ 3.5 | 94.9 $\pm$ 2.6 | 96.6 $\pm$ 1.7
widmanstatten | 52.9 $\pm$ 6.4 | 72.2 $\pm$ 7.1 | 66.8 $\pm$ 8.1
overall | 75.4 $\pm$ 3.7 | 92.6 $\pm$ 2.5 | 92.6 $\pm$ 2.5

These quantitative metrics are useful for interpreting the strengths and
weaknesses of a particular CNN model, but they do not necessarily directly
quantify the quality of the predicted segmentation maps due to inherent
subjectivity and bias in the labeling process. Even a single human annotator
will not be able to consistently label an entire dataset, especially for
ambiguous higher-level microconstituents such as the spheroidite class. For
example, the annotator must decide how closely to track cementite particles
when tracing out the edge of the denuded zone. In some cases, it is unclear
whether a carbide should be labeled as grain boundary cementite or as a piece
of Widmanstätten lath.

Furthermore, the low resolution of the input images relative to some of the
finer features of interest also places a practical upper bound on these
numerical performance scores, especially for microconstituents with large
interfacial areas like the Widmanstätten lath. Many of the Widmanstätten lath
in this dataset are just a few pixels wide, which can lead large shifts in
numerical scores for what a human might consider a minor difference in
labeling (e.g. dilating or eroding the Widmanstätten lath by one pixel).

###  3.2 Spheroidite particle segmentation

Figure 3 shows some validation results for the individual particle
segmentation task, with numerical performance reported in Tables 3 and 4;
additional examples are included in the online supplemental materials.
Particle predictions are overlaid in red on the input micrographs (a-d). The
second row (e-h) shows the empirical particle size distributions for both
particle predictions and annotations, as well as the results of the two-sample
Kolmogorov-Smirnov hypothesis test for distribution equivalence. Predictions
for larger particles relative to the image frame (Figures 3 b and c) are
consistently good, even where contrast gradients across particles and non-
trivial background structure challenge thresholding and edge-based
segmentation methods. The primary failure mode of the particle segmentation
model is underprediction of very small particles, particularly in Figure 3 a
and d. The vast majority of the fine particles in Figure 3 are missing
entirely, and many are only partially labeled by the CNN with just one or two
foreground pixels. These particles are typically one to five pixels in size,
suggesting that higher- or multi-resolution inputs are necessary for general
microstructure segmentation CNNs. However, the CNN does avoid spuriously
labeling the small segments of Widmanstätten in Figure 3 as particles.

Figure 3: (a-d) Validation set predictions for the spheroidite particle
segmentation task, along with (e-h) corresponding derived particle size
distributions for the particle predictions (blue) and annotations (green).
Scale bars indicate $5\mu m$.

The PixelNet model performs slightly better than Otsu’s thresholding method
[37] on all metrics. One source of bias in these performance measurements are
missing particles in the annotations, either from the removal of particles
intersecting the image border, or from failure of the semi-automated
annotation method itself. An additional source of bias stems from the
application of the watershed algorithm [38] to split conjoined particles in
the annotations; watershed segmentation is not presently applied to the
particle predictions, increasing the relative rate of larger particles.

Table 3: Particle segmentation performance averaged over validation images. Uncertainties are standard errors calculated across validation images. | IU | precision | recall
---|---|---|---
matrix | 90.0 $\pm$ 1.0 | 95.0 $\pm$ 0.6 | 94.5 $\pm$ 1.1
spheroidite | 54.8 $\pm$ 3.4 | 74.6 $\pm$ 2.8 | 70.3 $\pm$ 4.3
overall | 72.4 $\pm$ 3.1 | 91.1 $\pm$ 0.9 | 91.1 $\pm$ 0.9
Table 4: Particle segmentation performance metrics. Uncertainties are standard errors calculated across validation images. model | IU${}_{\text{matrix}}$ | IU${}_{\text{spheroidite}}$ | IU${}_{\text{avg}}$ | AC | PSD KS
---|---|---|---|---|---
otsu | 86.2 $\pm$ 7.2 | 53.7 $\pm$ 12.1 | 69.9 $\pm$ 9.3 | 88.1 $\pm$ 6.1 | -
pixelnet | 90.0 $\pm$ 1.0 | 54.8 $\pm$ 3.4 | 72.4 $\pm$ 3.1 | 91.1 $\pm$ 0.9 | 0.042

Despite good numerical performance on the particle segmentation task, the KS
test suggests we reject the null hypothesis that the predicted and ground
truth particle size distributions are equivalent for all but one of the 24
validation micrographs (shown in 3 b). The difficulty in detecting small
particles explains the discrepancies between empirical particle size
distributions that contribute to the KS score. For the two validation
micrographs in Figure 3 containing fine particles, the particle size
histograms and prediction maps show that the model often entirely misses
particles with radii smaller than 5px. Many of these missing ~5px particles
are partially labeled in the CNN predictions, leading to a severe
overrepresentation of single-pixel particles, especially in Figure 3 h.

###  3.3 Quantitative analysis of higher-order features

High-quality automated segmentation techniques for complex microstructure
constituents expand the scope of conventional quantitative microstructure
analysis by reducing the manual labor required to obtain statistically
meaningful amounts of data. In our UHCS case study, the CNN segmentation model
allows us to collect volume and shape statistics for the proeutectoid carbide
network, spheroidite particles, and Widmanstätten lath directly from SEM
micrographs with no manual intervention. Additionally, the microconstituent
prediction maps enable automated acquisition of interesting microstructural
statistics that were previously intractable, such as particle size
distributions conditioned on spatial relationships with other microstructure
features, or denuded zone widths [15].

Combining the two microstructure segmentation models allows us to filter out
irrelevant microstructure features in order to estimate particle size
distributions. Figure 4 shows combined microstructure predictions from both
the abstract microstructure model and the particle model, using the same color
scheme as Figures 2 and 3. We run the input image through separately-trained
particle segmentation CNN and microconstituent CNN, suppressing particle
predictions (red) outside of the predicted spheroidite regions (yellow). With
an appropriate number of images, one could also compute particle size
distributions spatially conditioned on other microstructure features (e.g.
distance from the network phase), which could help lead to insights into
operative microstructure evolution mechanisms (particle coarsening vs
precipitation). The resolution of these input micrographs is insufficient to
yield quantitatively accurate particle size distributions, especially with the
underprediction of small particles discussed in Section 3.2, as evident in
Figures 4 b and c. However, higher quality input and training micrographs will
mitigate this effect.

Figure 4: (a-d) Micrographs with (e-h) validation set microconstituent
predictions and (i-l) derived particle size distributions obtained by applying
the particle segmentation CNN to the semantic microstructure segmentation
dataset. Scale bars indicate $10\mu m$.

Figure 5 shows the predicted network and denuded zone boundaries for four
validation images with corresponding computed denuded zone width
distributions. The denuded zone width distributions are calculated by
aggregating the minimum distance to the network interface for each pixel on
the denuded zone boundary, as described in detail in Section 2.4. Generally,
these empirical denuded zone widths are reasonable, but some care is required
to interpret them. Specifically, the denuded zone width distributions in
Figures 5 b and d have high frequencies at small spacings that result from
spurious cementite network predictions. Figures 5 a and d also exhibit some
overprediction of the denuded zone width where the particles are very fine,
particularly in the upper portion of Figure 5 a.

Figure 5: (a-d) Validation set microconstituent predictions with (e-h)
corresponding denuded zone width distributions. The network interface is shown
in blue and the particle matrix interface is shown in yellow. Scale bars
indicate $10\mu m$

The initial investment of micrograph annotation and training a CNN makes sense
where a statistical number of samples must be characterized in the context of
alloy and processing optimization studies, and in the context of
microstructure and process validation or verification. Microconstituent
annotation accounted for a substantial portion of the time we spent on this
project. After an initial learning curve, a typical micrograph in our dataset
cost between 20 and 30 minutes to annotate. In contrast, a commodity GPU
performs Pixelnet-based segmentation at a rate of approximately one micrograph
per second, after an initial training period of two to four hours depending on
training hyperparameters. Success in a practical microstructure science
setting will depend on establishing higher-quality training data and deeper
understanding of the biases and variance of the labeling process.

The CNN predictions provide some useful feedback on these subjective labeling
decisions: consider the micrograph, annotation, and predictions in
supplemental Figure S1.6 a, e, and i. In the bottom half of this micrograph
(and in the other micrographs in this validation set), the annotator neglected
to label the metal matrix surrounding the Widmanstätten lath as such, while
the CNN consistently includes some matrix predictions associated with
Widmanstätten predictions. This subjective labeling decision can be mitigated
with higher-fidelity labeling of individual carbide particles – at much
greater labeling expense. A high quality dataset might be obtained via crowd-
sourcing (e.g. students in a microstructure analytics course), generation of
realistic synthetic datasets through e.g. phase field modeling, or through the
substantial expense of high-resolution elemental mapping with SEM+EDS (Energy-
dispersive spectrometry). A large dataset might also be collected in a semi-
supervised fashion through the development of smart microscopes with
integrated microstructure recognition features.

Furthermore, it is critical to benchmark microstructure-specific tasks against
other popular CNN architectures for semantic segmentation. Our approach of
directly transferring the particle prediction CNN is tenuous, especially due
to the disparity in magnification between the general UHCS and specific
particle segmentation datasets. Rather than training two separate CNNs, it may
be more appropriate train a single CNN in a multi-task setting, so that
microstructures are mapped to a common numerical representation before the
respective microconstituent and particle classification tasks.

Finally, microstructure data science is extremely data-limited in comparison
to most general computer vision tasks. Though outside the scope of the present
report, a detailed follow-on study to fully characterize the training data
requirements of deep learning based microstructure segmentation models would
be a valuable tool to enable experimental planning before significant
investment for industrial application. In parallel, collaboration with
computer scientists working on low-data deep learning, semi-supervised, and
unsupervised techniques could also open the door to applicability in many more
microstructure systems, especially where pixel-level annotations are expensive
or difficult to consistently obtain.

##  4 Conclusions

We demonstrate microstructural segmentation and quantitative analysis at a
high level of abstraction by applying an off-the-shelf deep neural network
architecture for pixel-wise prediction tasks. We also present two new open
microstructure segmentation benchmark datasets featuring the microstructures
in ultra-high carbon steel at different length scales. This data-driven
approach to microstructure segmentation expands the reach of traditional
quantitative microstructure characterization to more complex industrially-
relevant microstructure features that have, until now been, difficult to treat
in an automated fashion. Combined with emerging automated microscopy
capabilities, data-driven microstructure segmentation systems will enable
future applications in high-throughput microstructure studies, including
investigations of structure/processing relationships, microstructure design
and optimization, and microstructure-based material qualification.

## Acknowledgments

We gratefully acknowledge funding for this work through National Science
Foundation grants DMR-1507830 and CMMI-1826218, and through the John and
Claire Bertucci Foundation. We also recognize the support of the National
Institute of Standards and Technology and the National Research Council
Research Associate Program. The UHCS micrographs were graciously provided by
Matthew Hecht, Yoosuf Picard, and Bryan Webler (CMU) [1]. Semantic
microstructure annotations were performed by B.D. The spheroidite annotations
were graciously provided by Matthew Hecht and Txai Sibley. The open source
software projects Scikit-Learn [39], scikit-image [40], MITK [33], and keras
[41] were essential to this work.

## References

  * [1] Brian L. DeCost et al.  “UHCSDB (UltraHigh Carbon Steel micrograph DataBase): tools for exploring large heterogeneous microstructure datasets”  In _Integrating Materials and Manufacturing Innovation_ 6 Springer, 2017, pp. 197–205  DOI: 10.1007/s40192-017-0097-0
  * [2] Matthew D. Hecht et al.  “Ultrahigh Carbon Steel Micrographs”, https://hdl.handle.net/11256/940, 2017
  * [3] Brian L. DeCost et al.  “Ultrahigh Carbon Steel Micrographs: Microconstituent annotations”, https://hdl.handle.net/11256/964, 2018
  * [4] EO Hall  “The deformation and ageing of mild steel: III discussion of results”  In _Proceedings of the Physical Society. Section B_ 64.9 IOP Publishing, 1951, pp. 747
  * [5] C Zener  “quoted by CS Smith” quoted by CS Smith  In _Trans. AIME_ 175, 1948, pp. 15
  * [6] Yann LeCun, Yoshua Bengio and Geoffrey Hinton  “Deep learning”  In _Nature_ 521.7553 Springer Nature, 2015, pp. 436–444  DOI: 10.1038/nature14539
  * [7] Ian Goodfellow, Yoshua Bengio and Aaron Courville  “Deep Learning” http://www.deeplearningbook.org MIT Press, 2016
  * [8] Brian L. DeCost and Elizabeth A. Holm  “A computer vision approach for automated analysis and classification of microstructural image data”  In _Computational Materials Science_ 110 Elsevier, 2015, pp. 126–133  DOI: 10.1016/j.commatsci.2015.08.011
  * [9] Aritra Chowdhury, Elizabeth Kautz, Bülent Yener and Daniel Lewis  “Image driven machine learning methods for microstructure recognition”  In _Computational Materials Science_ 123 Elsevier, 2016, pp. 176–187
  * [10] Nicholas Lubbers, Turab Lookman and Kipton Barros  “Inferring Low-Dimensional Microstructure Representations Using Convolutional Neural Networks”  In _CoRR_ , 2016  arXiv: http://arxiv.org/abs/1611.02764v1
  * [11] Brian L. DeCost, Toby Francis and Elizabeth A. Holm  “Exploring the microstructure manifold: image texture representations applied to ultrahigh carbon steel microstructures”  In _Acta Materialia_ 133 Elsevier, 2017, pp. 30–40  DOI: 10.1016/j.actamat.2017.05.014
  * [12] Seyed Majid Azimi et al.  “Advanced Steel Microstructural Classification by Deep Learning Methods”  In _Scientific Reports_ 8.1 Springer Nature, 2018  DOI: 10.1038/s41598-018-20037-5
  * [13] Aayush Bansal et al.  “Pixelnet: Representation of the Pixels, By the Pixels, and for the Pixels”  In _CoRR_ , 2017  arXiv: http://arxiv.org/abs/1702.06506v1
  * [14] Matthew D Hecht, Bryan A Webler and Yoosuf N Picard  “Digital image analysis to quantify carbide networks in ultrahigh carbon steels”  In _Materials Characterization_ 117 Elsevier, 2016, pp. 134–143
  * [15] Matthew D Hecht, Yoosuf N Picard and Bryan A Webler  “Coarsening of inter and intragranular proeutectoid cementite in an initially pearlitic 2C-4Cr ultrahigh carbon steel”  In _Metallurgical and Materials Transactions A_ Springer, 2017, pp. 1–16  DOI: 10.1007/s11661-017-4012-2
  * [16] Panqu Wang et al.  “Understanding Convolution for Semantic Segmentation”  In _CoRR_ , 2017  arXiv: http://arxiv.org/abs/1702.08502v1
  * [17] V. Badrinarayanan, A. Kendall and R. Cipolla  “SegNet: A Deep Convolutional Encoder-Decoder Architecture for Scene Segmentation”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ PP.99, 2017, pp. 1–1  DOI: 10.1109/TPAMI.2016.2644615
  * [18] A. Bansal, B. Russell and A. Gupta  “Marr Revisited: 2D-3D Alignment via Surface Normal Prediction”  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 5965–5974  DOI: 10.1109/CVPR.2016.642
  * [19] Bharath Hariharan, Pablo Arbelaez, Ross Girshick and Jitendra Malik  “Hypercolumns for Object Segmentation and Fine-Grained Localization”  In _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015
  * [20] Jonathan Long, Evan Shelhamer and Trevor Darrell  “Fully convolutional networks for semantic segmentation”  In _2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ IEEE, 2015  DOI: 10.1109/cvpr.2015.7298965
  * [21] Alex Kendall, Vijay Badrinarayanan and Roberto Cipolla  “Bayesian Segnet: Model Uncertainty in Deep Convolutional Encoder-Decoder Architectures for Scene Understanding”  In _CoRR_ , 2015  arXiv: http://arxiv.org/abs/1511.02680v2
  * [22] Olaf Ronneberger, Philipp Fischer and Thomas Brox  “U-Net: Convolutional Networks for Biomedical Image Segmentation”  In _Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015_ Springer International Publishing, 2015, pp. 234–241  DOI: 10.1007/978-3-319-24574-4_28
  * [23] Simon Jégou et al.  “The One Hundred Layers Tiramisu: Fully Convolutional DenseNets for Semantic Segmentation”  In _CoRR_ , 2016  arXiv: http://arxiv.org/abs/1611.09326v2
  * [24] Karen Simonyan and Andrew Zisserman  “Very Deep Convolutional Networks for Large-Scale Image Recognition”  In _CoRR_ abs/1409.1556, 2014  URL: http://arxiv.org/abs/1409.1556
  * [25] Vinod Nair and Geoffrey E. Hinton  “Rectified Linear Units Improve Restricted Boltzmann Machines”  In _Proceedings of the 27th International Conference on International Conference on Machine Learning_ , ICML’10  Haifa, Israel: Omnipress, 2010, pp. 807–814  URL: http://dl.acm.org/citation.cfm?id=3104322.3104425
  * [26] Sergey Ioffe and Christian Szegedy  “Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift”  In _Proceedings of the 32nd International Conference on Machine Learning_ 37, Proceedings of Machine Learning Research  Lille, France: PMLR, 2015, pp. 448–456  URL: http://proceedings.mlr.press/v37/ioffe15.html
  * [27] Olga Russakovsky et al.  “ImageNet Large Scale Visual Recognition Challenge”  In _International Journal of Computer Vision_ 115.3 Springer Nature, 2015, pp. 211–252  DOI: 10.1007/s11263-015-0816-y
  * [28] Kaiming He, Xiangyu Zhang, Shaoqing Ren and Jian Sun  “Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification”  In _2015 IEEE International Conference on Computer Vision (ICCV)_ IEEE, 2015  DOI: 10.1109/iccv.2015.123
  * [29] Nitish Srivastava et al.  “Dropout: a simple way to prevent neural networks from overfitting”  In _Journal of Machine Learning Research_ 15.1, 2014, pp. 1929–1958  URL: http://dl.acm.org/citation.cfm?id=2670313
  * [30] Ilya Loshchilov and Frank Hutter  “Fixing Weight Decay Regularization in Adam”  In _CoRR_ , 2017  arXiv: http://arxiv.org/abs/1711.05101v2
  * [31] Diederik P. Kingma and Jimmy Ba  “Adam: A Method for Stochastic Optimization”  In _CoRR_ , 2014  arXiv: http://arxiv.org/abs/1412.6980v9
  * [32] Tsung-Yi Lin et al.  “Focal Loss for Dense Object Detection”  In _CoRR_ , 2017  arXiv: http://arxiv.org/abs/1708.02002v2
  * [33] Marco Nolden et al.  “The Medical Imaging Interaction Toolkit: Challenges and Advances”  In _International Journal of Computer Assisted Radiology and Surgery_ 8.4, 2013, pp. 607–620  DOI: 10.1007/s11548-013-0840-8
  * [34] JOHN CANNY  “A Computational Approach to Edge Detection”  In _Readings in Computer Vision_ Elsevier, 1987, pp. 184–203  DOI: 10.1016/b978-0-08-051581-6.50024-6
  * [35] Trevor Hastie, Jerome Friedman and Robert Tibshirani  “The Elements of Statistical Learning”  In _Springer Series in Statistics_ Springer New York, 2001  DOI: 10.1007/978-0-387-21606-5
  * [36] Frank J. Massey  “The Kolmogorov-Smirnov Test for Goodness of Fit”  In _Journal of the American Statistical Association_ 46.253 Informa UK Limited, 1951, pp. 68–78  DOI: 10.1080/01621459.1951.10500769
  * [37] Nobuyuki Otsu  “A Threshold Selection Method from Gray-Level Histograms”  In _IEEE Transactions on Systems, Man, and Cybernetics_ 9.1 Institute of ElectricalElectronics Engineers (IEEE), 1979, pp. 62–66  DOI: 10.1109/tsmc.1979.4310076
  * [38] L. Vincent and P. Soille  “Watersheds in digital spaces: an efficient algorithm based on immersion simulations”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 13.6 Institute of ElectricalElectronics Engineers (IEEE), 1991, pp. 583–598  DOI: 10.1109/34.87344
  * [39] Fabian Pedregosa et al.  “Scikit-learn: Machine learning in Python”  In _The Journal of Machine Learning Research_ 12 JMLR, 2011, pp. 2825–2830
  * [40] Stéfan Walt et al.  “Scikit-Image: Image Processing in Python”  In _PeerJ_ 2.nil, 2014, pp. e453  DOI: 10.7717/peerj.453
  * [41] François Chollet  “Keras”  GitHub, https://github.com/fchollet/keras, 2015