IEEEexample:BSTcontrol

# Phase Unwrapping of Color Doppler Echocardiography using Deep Learning

Hang Jung Ling  Olivier Bernard  Nicolas Ducros  and Damien Garcia Manuscript
received 14 April 2023; revised 5 June 2023; accepted 20 June 2023.This work
was supported in part by the MEGA Doctoral School (ED 162); in part by the
French National Research Agency (ANR) through the “4D-iVFM” Project under
Grant ANR-21-CE19-0034-01; in part by LABEX PRIMES, Université de Lyon, under
Grant ANR-11-LABX-0063; in part by LABEX CELYA within the program
“Investissements d’Avenir” operated by the French ANR under Grant
ANR-10-LABX-0060; and in part by GENCI-IDRIS (HPC resources) under Grant 2022-
[AD010313603]. (Corresponding author: Hang Jung Ling.)Hang Jung Ling, Olivier
Bernard, Nicolas Ducros and Damien Garcia are with CREATIS, CNRS UMR 5220,
INSERM U1294, INSA Lyon, University of Lyon 1, 69100 Villeurbanne, France
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>lyon.fr).Nicolas Ducros is also with Institut Universitaire de France (IUF),
France.

###### Abstract

Color Doppler echocardiography is a widely used non-invasive imaging modality
that provides real-time information about the intracardiac blood flow. In an
apical long-axis view of the left ventricle, color Doppler is subject to phase
wrapping, or aliasing, especially during cardiac filling and ejection. When
setting up quantitative methods based on color Doppler, it is necessary to
correct this wrapping artifact. We developed an unfolded primal-dual network
to unwrap (dealias) color Doppler echocardiographic images and compared its
effectiveness against two state-of-the-art segmentation approaches based on
nnU-Net and transformer models. We trained and evaluated the performance of
each method on an in-house dataset and found that the nnU-Net-based method
provided the best dealiased results, followed by the primal-dual approach and
the transformer-based technique. Noteworthy, the primal-dual network, which
had significantly fewer trainable parameters, performed competitively with
respect to the other two methods, demonstrating the high potential of deep
unfolding methods. Our results suggest that deep learning-based methods can
effectively remove aliasing artifacts in color Doppler echocardiographic
images, outperforming DeAN, a state-of-the-art semi-automatic technique.
Overall, our results show that deep learning-based methods have the potential
to effectively preprocess color Doppler images for downstream quantitative
analysis.

Color Doppler, Flow imaging, Echocardiography, Phase unwrapping, Dealiasing,
Deep learning, Deep unfolding, U-Net, Transformer, Primal-dual.

##  1 Introduction

Color Doppler ultrasound is a widely accepted clinical imaging modality for
non-invasive, real-time analysis of cardiovascular blood flow. While two-
dimensional color Doppler is commonly used for qualitative mapping of flow
characteristics, its applications for quantitative analysis are limited.
Common cardiovascular applications of color Doppler include detection of
valvular diseases [1] and septal defects [2], or guiding the positioning of
the pulsed-wave sample volume for spectral Doppler [3]. Among methods to make
color Doppler quantitative, Vector Flow Mapping has been introduced for
intracardiac flow dynamics. This method allows the computation of 2D or 3D
intraventricular velocity vector maps from color Doppler fields, using a
physically constrained optimization approach [4, 5]. Intracardiac vector flow
mapping from color Doppler requires two prerequisite steps: 1) delineation of
the endocardial inner wall, and 2) correction of wrapped (aliased) Doppler
regions.

With respect to the first step, Painchaud et al.[6] recently introduced a
2D+time deep learning architecture to enforce temporal consistency and
smoothness from one frame to the next. The second step of correction is
necessary due to the occurrence of aliasing, which is an artifact resulting
from insufficient slow-time sampling. This issue arises when the pulse
repetition frequency (PRF) is unable to capture high axial velocities
effectively. This causes the Doppler velocity to be wrapped to the opposite
side of the Doppler spectrum when its absolute value exceeds the Nyquist
velocity. Experienced clinicians can easily identify zones of aliasing in most
color Doppler images, where the color-coded velocities shift from red to blue
or vice versa. Aliasing can be removed in Doppler echocardiography by
designing multi-PRF sequences, as described by Posada et al.[7]. However, this
approach requires control of the ultrasound machine and is primarily suitable
for high-frame-rate echocardiography. When clinical scanners are used,
aliasing must be corrected by post-processing the color Doppler fields. While
a number of unwrapping algorithms have been proposed for dealiasing data maps
in atmospheric science, geodesy, and optical interferometry [8, 9, 10], this
problem has received less attention in color Doppler imaging.

Inspired by traditional radar approaches, Muth et al.[11] developed a
segmentation-based method for color Doppler dealiasing using statistical
region merging, called DeAN. This unsupervised method uses a scalar
hyperparameter to control the segmentation process. An optimal parameter was
determined from a supervised analysis of 50 color Doppler data. However, it
turns out that the DeAN method fails in difficult cases as shown in Fig. 11 in
[11] and that supervised corrections are still necessary in some situations.
With the goal of developing imaging tools that quantify blood flow from color
Doppler, we propose a deep learning (DL) approach to correct the aliased areas
of echocardiographic color Doppler maps. DL has been proposed for color
Doppler dealiasing in vascular flow imaging by Nahas et al.[12]. Their
approach aimed to solve the double aliasing problem using two U-Nets. The
first U-Net detected the presence of single aliasing while the second U-Net
was trained to identify and segment double-aliased pixels. They evaluated the
performance of their model by training it with different types of ultrasound
information. They found that the model trained with a combination of Doppler
frequency, power, and bandwidth performed the best for dealiasing in the
femoral bifurcation.

In our work, we focused on Doppler echocardiography. In contrast to vascular
flow imaging, cardiac color flow imaging can be subject to substantial clutter
signals originating from the myocardium and tending to spread the aliased
patterns. With the goal of proposing a robust DL method that correctly handles
aliasing in most situations, we developed and compared several architectures.
Our main contributions are:

  1. 1.

We designed a primal-dual network based on the idea of deep unfolding, and
compared it with state-of-the-art DL segmentation methods and DeAN.

  2. 2.

We used a private color Doppler echocardiographic dataset acquired in apical
three-chamber view (45 patients, 1,338 aliased and 2,379 non-aliased frames)
to train the neural networks and analyze their performance.

  3. 3.

We investigated the value of adding Doppler power as input information to
improve dealiasing.

  4. 4.

We introduced a data augmentation strategy that generates synthetic aliasing,
which solved the class imbalance problem and improved dealiasing performance
on difficult color Doppler images.

Highlights
---
$\bullet$ |  Our deep-unfolding-based primal-dual network (PDNet) incorporated the forward operator as prior information and had only 0.03M parameters.
$\bullet$ |  Our deep learning (DL) models outperformed a state-of-the-art non-DL approach in phase unwrapping of color Doppler echocardiography, with nnU-Net being the best candidate, followed by PDNet with 233 times fewer parameters.
$\bullet$ |  Automatic and accurate color Doppler echocardiographic phase unwrapping ensures correct visualization and enables quantification of intracardiac blood flow.

##  2 Methodology

Aliasing artifacts occur when axial blood speeds (velocity magnitudes) exceed
the Nyquist velocity $V_{N}$. The acquired Doppler velocity $V_{D}$ can be
written as a function of the unwrapped or alias-free Doppler velocity $V_{u}$
as follows:

$V_{D}=V_{u}-2\times n_{N}\,V_{N},$ (1)

where $n_{N}$ is an integer called the Nyquist number, which represents the
number of times the signal wraps around the Nyquist limit. The Nyquist number
reads (see [7] for the demonstration)

$n_{N}=\text{floor}\left(\frac{V_{u}+V_{N}}{2V_{N}}\right).$ (2)

Except for highly turbulent flows that may occur in transvalvular or
transseptal jets, there is no multiple aliasing in the adult left ventricle
scanned in the apical long-axis view, i.e., the integer $n_{N}$ belongs to
$\\{-1,0,1\\}$. Indeed, in adult echocardiography with a 3 MHz phased array,
Nyquist velocities typically range from 0.55 to 0.7 m/s. Thus, single (i.e.,
$n_{N}=-1\,\text{or}\,1$) or no (i.e., $n_{N}=0$) wrapping occurs as long as
the actual blood speed is less than 1.65-2.1 m/s (see (4)). It follows that
double aliasing does not occur in the left ventricle in most situations
without valvular disease or cardiac shunt. Equations (1) and (2) can be
rewritten to express $V_{D}$ as a wrapped version of $V_{u}$:

$V_{D}=K(V_{u})=(V_{u}+V_{N})\,\text{mod}\,(2V_{N})-V_{N},$ (3)

where mod is the modulo operation. In particular, for $n_{N}\in\\{-1,0,1\\}$,
the wrapping function $K$ becomes

$V_{D}=K(V_{u})=\small{\begin{cases}V_{u}-2V_{N}&\text{if
}V_{N}<V_{u}<3\,V_{N}\\\ V_{u}&\text{if }-V_{N}\leq V_{u}\leq V_{N}\\\
V_{u}+2V_{N}&\text{if }-3\,V_{N}<V_{u}<-V_{N}\end{cases}.}$ (4)

This representation implies that the dealiasing problem can be approached in
two different ways: _i)_ a problem that inverts the wrapping function (3) and
recovers $V_{u}$ from the Doppler velocities $V_{D}$ by changing absolute
jumps greater than $V_{N}$ to their $2{\times}V_{N}$ complement; _ii)_ a
multi-class segmentation approach that assigns a Nyquist number $n_{N}$ (2) to
each pixel of the input image then computes the actual unwrapped velocities
using (1). We investigated three deep learning (DL) models for dealiasing
color Doppler. The first method was derived from the deep unfolding/unrolling
framework and solved the inverse problem defined in (3) to estimate the actual
velocities. We faced a nonlinear inverse problem on non-trivial data, whose
solution may contain phase jumps at the blood/myocardium interfaces. Unrolled
methods are well suited for solving inverse problems. Primal-dual
optimization, on the other hand, is useful for nonlinear problems. For these
reasons, we tested a learned primal-dual algorithm inspired by Adler et
al.[13], as described in the following subsection. The other two methods were
state-of-the-art networks that we have adapted for determining Nyquist numbers
in color Doppler images. Fig. 1 illustrates the pipeline used for all three
methods, whose input data were the Doppler velocity multiplied by the Doppler
power before scan-conversion.

Figure 1: Pipeline of the deep learning-based methods for color Doppler
dealiasing.

###  2.1 OriPDNet: A Primal-Dual-based Deep Unfolding Network to Solve Inverse
Problems

To solve our nonlinear inverse problem (3), we used OriPDNet (refer to Fig. 2
in [13] for the network architecture), a deep unfolding network based on a
primal-dual optimization scheme [13]. Given a general inverse problem aiming
to obtain the solution $f$ from the measurement $g$ with the forward operator
$K$:

$g=K(f),$ (5)

the outline of OriPDNet to solve this problem is presented in Algorithm 1.

Algorithm 1 OriPDNet: Original primal-dual network

Initialize $\boldsymbol{f_{0}}$,
$\boldsymbol{h_{0}}=[\textbf{0},\textbf{0},\textbf{0},\textbf{0},\textbf{0}]\in\mathbb{R}^{M\times
N\times 5}$

for $i=1,\cdots,I$ do

$\boldsymbol{h_{i}}\leftarrow\Gamma_{\theta^{d}_{i}}(\boldsymbol{h_{i-1}},K(\boldsymbol{f^{(2)}_{i-1}}),\boldsymbol{g})$

$\boldsymbol{f_{i}}\leftarrow\Lambda_{\theta^{p}_{i}}(\boldsymbol{f_{i-1}},[\partial
K(\boldsymbol{f^{(1)}_{i-1})}]^{*}\,(\boldsymbol{h^{(1)}_{i}}))$

end for

return $\boldsymbol{f^{(1)}_{I}}$

OriPDNet involves several variables and operators, including the forward
operator $K$, the adjoint of its Fréchet derivative $[\partial K]^{*}$, the
input measured data $g$, the primal and dual variables $f_{i}$ and $h_{i}$,
and the learned primal and dual proximal operators $\Lambda_{\theta^{p}_{i}}$
and $\Gamma_{\theta^{d}_{i}}$. Convolutional layers are used to learn these
proximal operators. The hyperparameter $I$, which determines the number of
iterations, requires careful tuning for each specific problem. The primal and
dual variables, $f_{i}$ and $h_{i}$, are initialized then iteratively updated
using the learned primal and dual proximal operators,
$\Lambda_{\theta^{p}_{i}}$ and $\Gamma_{\theta^{d}_{i}}$. The solution to the
inverse problem (5) is obtained by extracting the first element of the primal
variables, $f_{I}^{(1)}$. In [13], the authors recommended setting the
dimension of the primal and dual spaces to five as the best compromise between
memory usage and reconstruction quality, i.e.,
$f_{i}=[f_{i}^{(1)},f_{i}^{(2)},f_{i}^{(3)},f_{i}^{(4)},f_{i}^{(5)}]\in\mathbb{R}^{M\times
N\times 5}$ and
$h_{i}=[h_{i}^{(1)},h_{i}^{(2)},h_{i}^{(3)},h_{i}^{(4)},h_{i}^{(5)}]\in\mathbb{R}^{M\times
N\times 5}$, where $(M\times N)$ is the size of the input data. We conducted
preliminary testing and validated the use of five spaces, in accordance with
their suggestion.

###  2.2 PDNet: A Deep Unfolding Network for Color Doppler Dealiasing

To deal with our specific inverse problem (3) for color Doppler dealiasing, we
adapted OriPDNet. The modified version was named as PDNet and summarized in
Algorithm 2, with the main changes highlighted in blue.

Algorithm 2 PDNet: Proposed primal-dual network

Initialize $\boldsymbol{{\color[rgb]{0,0,1}V_{0}}}$,
$\boldsymbol{h_{0}=[0,0,0,0,0]}\in\mathbb{R}^{M\times N\times 5}$

for $i=1,\cdots,I$ do

$\boldsymbol{h_{i}}\leftarrow{\color[rgb]{0,0,1}\Gamma_{\theta^{d}}}(\boldsymbol{h_{i-1}},{\color[rgb]{0,0,1}K}(\boldsymbol{{\color[rgb]{0,0,1}V^{(2)}_{i-1}}}),\boldsymbol{{\color[rgb]{0,0,1}V_{D}}})$

$\boldsymbol{{\color[rgb]{0,0,1}V_{i}}}\leftarrow{\color[rgb]{0,0,1}\Lambda_{\theta^{p}}}(\boldsymbol{{\color[rgb]{0,0,1}V_{i-1}}},\boldsymbol{h^{(1)}_{i}})$

end for

return
${\color[rgb]{0,0,1}C_{\theta}}(\boldsymbol{{\color[rgb]{0,0,1}V^{(1)}_{I}}})$

with $\boldsymbol{V_{u}}=\boldsymbol{{\color[rgb]{0,0,1}V^{(1)}_{I}}}$ and
$\boldsymbol{n_{N}}={\color[rgb]{0,0,1}C_{\theta}}(\boldsymbol{{\color[rgb]{0,0,1}V^{(1)}_{I}}})$.

Specifically, we defined the forward operator $K$ as a wrapping function given
by (3). Despite the discontinuity of this function at each $V=V_{N}\pm
2kV_{N}$ (with $k\in\mathbb{Z}^{*}$), its derivative was an identity function,
i.e., $\partial K(V)=\text{id}$. Thus, its adjoint $[\partial K(V)]^{*}$ was
also an identity function. Unlike the original approach (see Algorithm 1), we
used the same feature maps for each iteration of the main loop, which
significantly reduced the number of parameters to learn (30,000 instead of
30,000$\times I$, with $I$ being the number of iterations) while maintaining
the same accuracy. We made this change to avoid training instabilities that we
observed while experimenting with OriPDNet. We also added a convolutional
layer $C_{\theta}$ at the end of the network to output the Nyquist number from
the estimated velocities $V_{I}^{(1)}$. The main reason for this was to avoid
non-integer Nyquist numbers due to the regressed velocities. For a fair
comparison between PDNet and OriPDNet, the same convolution layer $C_{\theta}$
was also applied to the output of OriPDNet.

###  2.3 Segmentation Networks for Color Doppler Dealiasing

Table 1: Main configurations of the three methods evaluated in this study. Lowest resolution: Size of the lowest resolution of feature maps in pixels. Down. scheme: Downsampling scheme. Up. scheme: Upsampling scheme. Optimization scheme: Optimizer + initial learning rate (+ learning rate scheduler used). # param.: Number of trainable parameters. Methods | Number of | Lowest | Down. | Up. | Normalization | Batch | Optimization | Loss | # param.
---|---|---|---|---|---|---|---|---|---
feature maps | resolution | scheme | scheme | scheme | size | scheme | function
PDNet | $32\rightarrow 32\rightarrow 5$ | 192$\times$40 | - | - | - | 4 |  | ADAM[14] +
---
0.001 + cosine
annealing
| Cross entropy
---
\+ Dice
0.03M
nnU-Net | $32\downarrow 480\uparrow 32$ | 12$\times$5 |  | Stride
---
pooling
| Transposed
---
conv.
| Instance
---
norm.
4 |  | SGD + 0.01 +
---
polynomial
decay
| Cross entropy
---
\+ Dice
7M
BATFormer | $16\downarrow 256\uparrow 16$ | 16$\times$16 |  | Max
---
pooling
| 2 $\times$ 2
---
repeats
| Batch
---
norm.
4 |  | ADAM +
---
0.001
| Cross entropy
---
\+ Dice +
smooth L1
1.2M

nnU-Net is currently one of the best performing approaches for medical image
segmentation [15]. This model is based on the U-Net architecture and
implements several successful DL tricks, such as automatic hyperparameter
search of the U-Net architecture to increase accuracy, a patch-wise approach
to preserve image resolution, a deep supervision strategy to maintain accuracy
at all scales, and data augmentation during both training and inference to
enforce generalization. In this study, we addressed the dealiasing of color
Doppler as a 3-class segmentation problem with nnU-Net, where each class
corresponded to a Nyquist number $n_{N}\in\\{-1,0,1\\}$. Our network included
four stages in the encoder/decoder parts and had an input size of 192 $\times$
40 pixels, which was the median image size of our dataset. Table 1 provides
more details about the architecture and the training scheme (see Fig. C.1 in
[15] for an illustration of the nnU-Net architecture).

Recently, transformer-based approaches have been shown to outperform the nnU-
Net model in some medical challenges [16]. These models attempt to solve the
segmentation problem in a different way, using attention mechanisms with
receptive fields that cover the entire image. Among the best-performing
models, we chose to train the BATFormer architecture [17] for the color
Doppler dealiasing task, using a segmentation-based technique. This model
employs a multiscale approach based on a U-Net architecture, with transformer
blocks added to the decoder. This strategy results in an efficient and
lightweight architecture (1.2M parameters) that is suitable for learning from
small to medium-sized datasets. The main configurations of BATFormer are
listed in Table 1, and an illustration of its architecture can be found in
Fig. 2 in [17].

###  2.4 Input Data Strategy

Color Doppler echocardiography produces two types of information: _i)_ Doppler
velocity, which can be corrupted by aliasing in regions of high blood speed,
and _ii)_ Doppler power, which provides insight into the regions where
velocity measurements are reliable. Using both Doppler power and velocity as
input to DL models allowed them to learn how to limit the dealiasing process
in regions of interest and identify ambiguous areas. Therefore, we performed
an ablation study to evaluate the potential improvement provided by Doppler
power. This study was conducted using the nnU-Net architecture, known for its
stability in training and optimal configurations. Specifically, three nnU-Net
models were trained with three combinations of input data: 1) nnU-Net #1
trained with Doppler velocity only, 2) nnU-Net #2 trained with the
concatenation of Doppler velocity with Doppler power, and 3) nnU-Net #3
trained with the multiplication of Doppler velocity by Doppler power. After
determining the best candidate for the input data, we trained the three
models, PDNet, nnU-Net, and BATFormer, using this input combination to compute
the Nyquist numbers $n_{N}$ (2) from which the unwrapped Doppler velocities
$V_{u}$ (1) were derived. Our goal was to increment or decrement the Doppler
velocities by $2n_{N}V_{N}$, not to modify them by smoothing, for example.

###  2.5 Artificial Aliasing Augmentation Strategy

Color Doppler images may exhibit aliasing only in localized regions or frames,
resulting in datasets that are often imbalanced, with most pixels belonging to
the background class (i.e., without aliasing). To address this issue, we used
standard data augmentation techniques such as rotation, flipping etc., during
training. We also proposed an additional data augmentation technique, which we
called artificial aliasing augmentation, to improve the generalizability of
our algorithms. This technique involved identifying regions with high Doppler
velocity and power on alias-free Doppler images, and applying a wrapping
function defined by (3) with a lower Nyquist velocity to create artificial
aliasing artifacts, followed by a normalization. The ground-truth references
of these artificially aliased frames were created on the fly by comparing the
Doppler velocities before and after this augmentation. By creating realistic
artificially aliased frames, as shown in Fig. 2, this strategy enabled us to
balance the classes in training batches. To evaluate the potential benefits of
artificial aliasing augmentation, we conducted an additional ablation study
where we tested the three DL models with and without this technique during
training.

Figure 2: Generation of synthetic images with artificial aliasing artifacts
(right column) from non-aliased images (left column).

###  2.6 DeAN: State-of-the-Art non DL-Based Dealiasing Method for CDI

DeAN, mentioned in Section 1, is currently one of the most powerful non DL-
based methods for color Doppler dealiasing. It is a semi-supervised method
with a hyperparameter, $Q$. To unwrap aliased pixels, DeAN first segments
color Doppler images using a region-merging scheme based on the Hoeffding’s
probability inequality. Then, DeAN compares each segmented region with its
nearest neighbors and performs dealiasing if necessary. This step is based on
the assumption that the largest segment is not aliased and is repeated until
all the segments have been analyzed. The main drawback of this method is the
need to manually search for the optimal $Q$ hyperparameter for each frame to
obtain the best dealiasing results. We compared the dealiasing performance of
DeAN with both the default $Q=10$ and with the manually optimized $Q$
hyperparameter, against the three DL methods.

###  2.7 Evaluation metrics

All three DL models were designed to output the Nyquist numbers $n_{N}$; the
dealiased velocity maps, $V_{u}$, were recovered using (1). To evaluate the
accuracy of the dealiased velocity maps and the Nyquist numbers outputted by
each method, we computed four evaluation metrics.

We compared the dealiased Doppler velocity maps $V_{u}$ with the ground-truth
alias-free Doppler velocity maps $V_{ref}$ by computing the cosine similarity
index:

$\text{Cosim}(V_{u},V_{ref})=\frac{V_{u}\cdot V_{ref}}{\lVert
V_{u}\rVert\cdot\lVert V_{ref}\rVert}$ (6)

Cosine similarity is a commonly used similarity measure for comparing text
data or images. We used this similarity index in our previous work on color
Doppler dealiasing [11]. In addition, we computed three classification metrics
to verify whether each pixel was classified correctly on color Doppler images.
The first classification metric was the balanced accuracy score, which is more
suitable for unbalanced datasets than the classical accuracy score. It was
calculated using the following formula:

$\text{Accuracy}=\frac{1}{2}\times\left(\frac{\text{TP}}{\text{TP}+\text{FN}}+\frac{\text{TN}}{\text{TN}+\text{FP}}\right)$
(7)

where TP, FN, TN and FP refer to true positives, false negatives, true
negatives, and false positives, respectively. Besides, the classical recall
($\frac{\text{TP}}{\text{TP}+\text{FN}}$) and precision
($\frac{\text{TP}}{\text{TP}+\text{FP}}$) metrics were also computed to
evaluate the overall performance of the methods.

To ensure the reliability and relevance of the results, we conducted a 9-fold
cross-validation to compute the scores presented in each table in Section 3.2.
For each fold, we split the dataset into training, validation, and test sets
using a ratio of 36/4/5 patients. This resulted in an average of 2,974, 330,
and 413 color Doppler frames for the training, validation, and test sets,
respectively.

##  3 Experiment Setup and Results

###  3.1 Dataset and Training Strategies

Figure 3: Generation of the ground-truth masks from the Doppler velocities.
The red, black, and green segmentation masks correspond to Nyquist numbers
$n_{N}=1$, $0$, and $-1$, respectively.

####  3.1.1 Color Doppler dataset

To evaluate the performance of our methods, we used a color Doppler
echocardiographic dataset of 45 patients that were acquired using a Vivid 7
ultrasound system (GE Healthcare, USA) and a GE 5S cardiac sector probe
(bandwidth = 2-5 MHz). Doppler velocity and power data prior to scan
conversion were extracted into HDF formats using EchoPAC software (GE
Healthcare). The EchoPAC software returned unitless power data in the range of
1-100. The power data, $P$, were compressed by taking the logarithm and then
scaled to [0,1]: $log(P)/2\leftarrow P$. The cardiologist used the default
settings (center frequency, pulse length, pulse repetition frequency, packet
size, clutter filter, etc.). These proprietary parameters are masked and could
not be extracted. In most cases, the Nyquist velocity was about 0.6 m/s.
Assuming a center frequency of 3 MHz, the pulse repetition frequency was
approximately 4500 Hz. These anonymized data came from a previous published
study [18]. The patients were examined in the echocardiography laboratory
under standard medical conditions. As a result, some patients had significant
heart disease, while others had no visible pathology. The random selection of
patients and their anonymization prevent us from knowing their demographic and
pathological status. The sequences were acquired in the apical three-chamber
view and included both Doppler velocity and power information. Each sequence
covered at least one complete cardiac cycle, resulting in a total of 1,338
aliased frames and 2,379 non-aliased frames. Since color Doppler has a
relatively low frame rate of 10 to 15 frames per second in clinical
echocardiography, we considered each frame to be independent. To avoid
interpolation artifacts, the data were collected and processed in polar
coordinates (i.e., before scan-conversion), but for better visualization, the
results were presented in Cartesian coordinates.

The training, validation, and test data sets, i.e., the pairs of original and
alias-free Doppler velocity maps, were generated by an experienced analyzer.
For this task, the non-scan-converted clinical Doppler maps were oversegmented
and labeled using a statistical growing region method (see Fig. 1.B in [11]).
The analyst manually identified the aliased regions, specifically focusing on
those related to intraventricular blood flow, which were then corrected by
applying $\pm 2V_{N}$. The noisy regions associated with low Doppler power
were left unchanged. Examples of color Doppler frames along with their
reference segmentations are illustrated in Fig. 3.

####  3.1.2 Training strategies

To train the deep learning (DL) methods described in Section 2, we performed
supervised learning using the ground-truth segmentations from our in-house
dataset. Besides applying the data augmentation strategies mentioned in
Section 2.5, we further addressed the class imbalance of our dataset by
ensuring that each batch contained at least one aliased image, whether real or
synthetic. The BATFormer model was designed using the official implementation
proposed in its GitHub repository111https://github.com/xianlin7/BATFormer.
This model took color Doppler images resized to 256 $\times$ 256 pixels as
input and was trained for 400 epochs. On the other hand, nnU-Net and PDNet
were implemented using the ASCENT framework222https://github.com/creatis-
myriad/ASCENT. For these two approaches, we used a patch-wise approach to
preserve the resolution of the input data. The models were trained for 1000
epochs to prevent any potential under/overfitting. More details on the
BATFormer, nnU-Net, and PDNet architectures are provided in Table 1.

###  3.2 Experimental Results

####  3.2.1 Doppler power information was useful in dealiasing difficult case

Table 2 reports the results of the ablation study aimed at identifying the
optimal combination of input data. The results indicate that the three nnU-Net
models performed similarly across all metrics, implying that incorporating
Doppler power information in the input data did not substantially improve the
models’ performance. However, upon evaluation on a challenging test set (right
part of Table 2), the model that was trained with the multiplication of
Doppler velocity and Doppler power (nnU-Net #3) demonstrated better
performance for all metrics except precision. The last two columns of Fig. 4
show two samples taken from the difficult fold, where the aliased and non-
aliased regions had similar hues. This made the correction of the aliased
velocities difficult. Thus, although not critical, using the Doppler velocity-
Doppler power product as input data is recommended as it can enhance the
models’ generalization ability, especially for challenging data. For
subsequent experiments and results, we trained all DL methods with this input
combination.

Table 2: Dealiasing by nnU-Net trained with different combinations of input data using 9-fold cross-validation. #1, #2, and #3 correspond to the use of Doppler velocity only, the concatenation of Doppler velocity with Doppler power, and the multiplication of the Doppler velocity by Doppler power as input data, respectively. The Difficult fold column shows the evaluation results of different nnU-Net models on a challenging test set containing color Doppler frames with aliased and non-aliased regions of similar hue. Methods | Full dataset (Number of frames = 3,717) | Difficult fold (Number of frames = 413)
---|---|---
Cosim | Accuracy | Recall | Precision | Cosim | Accuracy | Recall | Precision
$\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$
nnU-Net #1 velocity only | 0.99 | 0.95 | 0.90 | 0.92 | 0.98 | 0.93 | 0.86 | 0.95
$\pm$0.01 |  $\pm$0.03 |  $\pm$0.07 |  $\pm$0.06 |  $\pm$0.01 |  $\pm$0.06 |  $\pm$0.13 |  $\pm$0.01
nnU-Net #2 velocity-power concatenation | 0.99 | 0.95 | 0.91 | 0.91 | 0.98 | 0.94 | 0.88 | 0.92
$\pm$0.01 |  $\pm$0.03 |  $\pm$0.07 |  $\pm$0.05 |  $\pm$0.01 |  $\pm$0.07 |  $\pm$0.13 |  $\pm$0.05
nnU-Net #3 velocity-power multiplication | 0.99 | 0.95 | 0.89 | 0.90 | 0.99 | 0.94 | 0.89 | 0.92
$\pm$0.01 |  $\pm$0.06 |  $\pm$0.13 |  $\pm$0.11 |  $\pm$0.01 |  $\pm$0.07 |  $\pm$0.14 |  $\pm$0.02

####  3.2.2 PDNet outperformed its original counterpart

We conducted a study to determine the optimal number of iterations for
updating the primal and dual variables in both OriPDNet and PDNet, given the
sensitivity of this type of method to this parameter. The results are shown in
Table 3. From this table, we can see that OriPDNet reached a plateau after 10
iterations, and beyond 20 iterations, it became unstable during training and
failed to produce results. On the other hand, PDNet exhibited greater training
stability and reached a performance plateau after 20 iterations. Moreover,
PDNet achieved better optimal results than OriPDNet for all metrics except for
the cosine similarity index, where both methods performed equally well. These
results suggest that using the same feature maps in the primal-dual approach
is more suitable for the dealiasing task and support the use of our deep
unfolding network. Additionally, it is worth mentioning that PDNet had only
30,000 parameters, making it the lightest of the three DL models tested, as
detailed in Table 1. The significant reduction in parameters of PDNet as
compared with other models was due to the inclusion of the forward operator as
prior information and the use of the same feature maps per iteration.

Table 3: 9-fold cross-validation dealiasing results of PDNet trained with different number of iterations (# iter.) for updating the primal and dual variables. The results on the left were obtained with OriPDNet, i.e., different feature maps per iteration, while the results on the right were obtained with the proposed PDNet using the same feature maps for each iteration. # iter. |  OriPDNet: Different feature maps per iteration |  PDNet: Same feature maps for each iteration
---|---|---
Cosim | Accuracy | Recall | Precision | Cosim | Accuracy | Recall | Precision
$\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$
1 | 0.95 | 0.86 | 0.71 | 0.51 | 0.95 | 0.87 | 0.74 | 0.52
$\pm$0.03 |  $\pm$0.06 |  $\pm$0.13 |  $\pm$0.09 |  $\pm$0.03 |  $\pm$0.06 |  $\pm$0.12 |  $\pm$0.10
10 | 0.98 | 0.92 | 0.84 | 0.78 | 0.98 | 0.92 | 0.84 | 0.80
$\pm$0.02 |  $\pm$0.06 |  $\pm$0.13 |  $\pm$0.11 |  $\pm$0.02 |  $\pm$0.06 |  $\pm$0.12 |  $\pm$0.10
20 | 0.98 | 0.92 | 0.84 | 0.80 | 0.98 | 0.94 | 0.87 | 0.83
$\pm$0.02 |  $\pm$0.06 |  $\pm$0.12 |  $\pm$0.11 |  $\pm$0.01 |  $\pm$0.06 |  $\pm$0.13 |  $\pm$0.11
30 | - | - | - | - | 0.97 | 0.91 | 0.82 | 0.77
$\pm$0.02 |  $\pm$0.07 |  $\pm$0.13 |  $\pm$0.14

####  3.2.3 Artificial aliasing augmentation improved the performance of
segmentation-based networks

The results presented in Table 4 show that the use of artificial aliasing
augmentation during training had varying effects on the performance of the
different DL models. For nnU-Net, there was a slight improvement in all
metrics except precision. In contrast, BATFormer showed significant
improvement in accuracy, recall, and precision metrics, with values increasing
from 0.88, 0.76, and 0.85 to 0.91, 0.81, and 0.91, respectively. However, for
PDNet, the use of artificial aliasing augmentation resulted in degraded
performance, with accuracy and recall decreasing from 0.94 and 0.87 to 0.88
and 0.77, respectively. These results highlight the challenge that the primal-
dual-based regression methods face when generalizing to different types of
aliasing. Based on the findings shown in Tables 2 to 4, we determined that the
multiplication of Doppler velocity and power should be used as input for all
methods, while artificial aliasing augmentation should be applied during
training only for the segmentation-based techniques, i.e., nnU-Net and
BATFormer.

Table 4: 9-fold cross-validation dealiasing results of the three implemented deep learning solutions trained with and without the proposed artificial aliasing augmentation strategy. Methods | Without artificial aliasing augmentation | With artificial aliasing augmentation
---|---|---
Cosim | Accuracy | Recall | Precision | Cosim | Accuracy | Recall | Precision
$\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$
PDNet | 0.98 | 0.94 | 0.87 | 0.83 | 0.98 | 0.88 | 0.77 | 0.84
$\pm$0.02 |  $\pm$0.06 |  $\pm$0.12 |  $\pm$0.10 |  $\pm$0.01 |  $\pm$0.08 |  $\pm$0.16 |  $\pm$0.09
nnU-Net | 0.99 | 0.95 | 0.89 | 0.90 | 0.99 | 0.96 | 0.91 | 0.89
$\pm$0.01 |  $\pm$0.03 |  $\pm$0.09 |  $\pm$0.06 |  $\pm$0.01 |  $\pm$0.03 |  $\pm$0.06 |  $\pm$0.06
BATFormer | 0.98 | 0.88 | 0.76 | 0.85 | 0.98 | 0.91 | 0.81 | 0.91
$\pm$0.02 |  $\pm$0.08 |  $\pm$0.17 |  $\pm$0.11 |  $\pm$0.02 |  $\pm$0.07 |  $\pm$0.14 |  $\pm$0.07

####  3.2.4 nnU-Net gave the best dealiasing results

Table 5 presents the final results of our study, where we compare the three DL
methods with their optimal configurations against the DeAN algorithm. We
observe that all three DL methods outperformed the DeAN algorithm, even the
version with the manually chosen optimal $Q$ hyperparameter. This outcome
confirms the potential of DL methods for color Doppler dealiasing. Among the
DL methods, nnU-Net achieved the highest scores overall, with a cosine
similarity close to 1, an accuracy of 0.96, a recall of 0.91, and a precision
of 0.89. Therefore, we conclude that nnU-Net is the best DL approach currently
available for dealiasing tasks in echocardiography. Additionally, it is
interesting to note that BATFormer showed a clear improvement when we
increased the amount of synthetic data, indicating that this type of approach
requires a larger dataset to improve its performance for the dealiasing task.
Finally, it is worth noting that PDNet achieved promising results with 233
times fewer parameters compared to nnU-Net, highlighting the potential of
incorporating analytical context into the DL framework to regularize the
solution space.

We also provide a visual inspection of the performance of the various methods
on aliased images with different degrees of difficulty in Fig. 4. We can see
that the DL methods performed similarly well on the easy and moderate cases
(first two columns), but nnU-Net produced the closest results to the reference
on the more challenging case (third column). This finding is consistent with
the quantitative results presented in Table 5. The last column in Fig. 4 shows
an example where no method was able to handle aliasing correctly. This example
is similar to the one in the first column, but with a more pronounced level of
aliasing. In this particular case, the DeAN method gave the best results. This
suggests that it would be advisable to supplement the training set with
challenging configurations.

Table 5: 9-fold cross-validation final dealiasing results of DL methods with their best configurations against non-DL DeAN method. Methods | Cosim | Accuracy | Recall | Precision
---|---|---|---|---
$\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$ | $\pm\sigma$
DeAN (Q=10) | 0.95 | 0.81 | 0.62 | 0.53
$\pm$0.03 |  $\pm$0.08 |  $\pm$0.16 |  $\pm$0.20
DeAN (Optimized Q) | 0.98 | 0.91 | 0.83 | 0.80
$\pm$0.01 |  $\pm$0.04 |  $\pm$0.08 |  $\pm$0.13
PDNet | 0.98 | 0.94 | 0.87 | 0.83
$\pm$0.02 |  $\pm$0.06 |  $\pm$0.12 |  $\pm$0.10
nnU-Net | 0.99 | 0.96 | 0.91 | 0.89
$\pm$0.01 |  $\pm$0.03 |  $\pm$0.06 |  $\pm$0.06
BATFormer | 0.98 | 0.91 | 0.81 | 0.91
$\pm$0.02 |  $\pm$0.07 |  $\pm$0.14 |  $\pm$0.07
Figure 4: Color Doppler images acquired during: from left to right column,
late diastolic filling, systolic ejection, early diastolic filling, and
diastole (a failed case). They were dealiased by DL-based methods and by DeAN
with optimized $Q$ hyperparameters. First row: aliased raw color Doppler.
Second row: alias-free ground truth (GT).

##  4 Discussions

Color Doppler imaging takes high-pass-filtered in-phase/quadrature (I/Q) data
of the same region of interest acquired along the slow-time axis and
differentiates them pairwise using a lag-1 autocorrelator. The resulting maps
show blood displacements between two consecutive slow-time samples. By its
very nature, color Doppler imaging is an interferometric technique that
enables the measurement of displacements with a precision that can reach
fractions of the center wavelength. Similarly, synthetic aperture radar
interferometry (InSAR), a remote sensing technique used to map the Earth’s
surface deformations, generates interferograms that display ground-surface
displacements. Like color Doppler, most interferometric imaging techniques in
fields such as medical imaging, remote sensing, and optical metrology (e.g.,
phase-contrast MRI, InSAR, holographic interferometry) are subject to
aliasing, i.e., jumps that occur whenever the phase shift equals $\pm\pi$. Our
study aimed to address the issue of phase jumps. Among the traditional methods
for phase unwrapping, one can mention: i) graph cuts [19, 20], which involve
representing the wrapped phase data as a graph and determining the minimum cut
that separates the known and unknown phase values; ii) least-squares
approaches, which minimize the differences between partial derivatives of the
wrapped phase and those of the unwrapped solution [21, 8]. Specifically for
color Doppler echocardiography, Muth et al.developed DeAN, a dealiasing
algorithm based on statistical region merging [11], which was used in this
study for comparative purposes.

Recently, deep learning (DL) techniques have been used to improve traditional
methods in phase unwrapping [22, 23, 12, 24]. Our goal was to obtain alias-
free color Doppler echocardiography by applying DL to dealias clinical Doppler
velocity fields. DL-based approaches have been introduced for 2D phase
unwrapping in InSAR [24]. Unlike echocardiographic images, InSAR
interferograms are subject to multiple wraps, making the networks proposed in
our study not suited since we focused only on single aliasing. On the other
hand, InSAR images are not subject to significant clutter, whereas clutter in
Doppler echocardiography generates substantial noise near moving tissues,
making 2D phase unwrapping challenging. As a result, non-DL approaches such as
graph cuts or least-squares methods, which work well for InSAR interferograms,
are not effective for echocardiographic Doppler fields. Although the DeAN
technique largely solved the problem, it still fails in some situations, as
shown in our study. Therefore, we turned to DL and conducted an in-depth
analysis and comparison of three architectures, including PDNet, which
utilizes an unfolding framework. In addition, we illustrated the potential
benefits of incorporating Doppler power information, since low power generally
indicates poor blood Doppler signal. To better balance the aliased and non-
aliased input data during training, we resorted to data augmentation by
generating synthetic aliasing.

###  4.1 Comparison of the DL Methods

Our study found that the three DL methods we tested (PDNet, nnU-Net, and
BATFormer) outperformed the non-DL DeAN method for color Doppler dealiasing.
Notably, we observed that nnU-Net had the best performance, suggesting that
the 2D U-Net architecture used in nnU-Net may be particularly well-suited for
this task due to its ability to effectively capture spatial features. For
example, in a challenging case where the aliased and non-aliased regions had
similar hues (Fig. 4, third column), nnU-Net was able to unwrap correctly
while other DL methods failed or were less successful. Similar structures
corrected by an expert were part of the training dataset, which implies that
nnU-Net probably learned the flow patterns and leveraged this knowledge to
achieve successful outcomes. PDNet also performed well, requiring $>$ 200
times fewer parameters than nnU-Net, highlighting the potential for simpler DL
models to achieve competitive results in color Doppler processing. Further
exploration of this type of unfolding approach, including more complex
modeling of the forward operator, is needed.

Although the third input strategy (velocity-power multiplication) contained
less information than the second (velocity-power concatenation), it performed
slightly better in the difficult fold (last row of Table 2). The
multiplication strategy largely suppressed velocity discontinuities in noisy
regions, making the training task easier. In contrast, the concatenation of
Doppler power and velocity allowed the model to learn the best strategy for
combining these two inputs, which could be beneficial in larger datasets.

While we did not observe significant performance gains with BATFormer, the
addition of synthetic data improved the outcomes, indicating that BATFormer
also has potential for color Doppler dealiasing, especially when more data is
available. Our results demonstrate that DL methods can significantly improve
upon traditional methods for color Doppler processing. Additionally, they
underscore the importance of further investigating the performance of
different DL architectures for this task and finding ways to effectively
exploit the strengths of each architecture.

###  4.2 Limitations and Perspectives

Color Doppler aliasing in the left ventricle mainly occurs in the mitral jet
during early and late filling, as well as during ejection into the ventricular
outflow tract. As depicted in the figures, the aliasing in our study was
single. However, in certain valvular diseases, such as mitral stenosis or
aortic regurgitation, multiple aliasing can occur in the intraventricular
cavity due to the high fluctuating velocities of the turbulent jet. The nature
of multiple aliasing in this context differs from that observed in InSAR,
requiring specific studies to assess the feasibility of removing aliasing in
areas with significant local flow perturbations. Although this remains to be
verified, it is likely that a similar strategy could also work with disturbed
flows, provided that we have access to Doppler data with their alias-free
references. Such ground truths could be obtained by supervised correction, as
in this study, and by simulations [25].

Since we used a clinical ultrasound system with a color Doppler rate of 10 to
15 frames per second, our study did not exploit temporal information. In the
context of high-frame-rate echocardiography [26], neural networks with
enforced temporal consistency [6] or 3D U-Net could potentially improve
dealiasing performance by leveraging temporal information. This approach would
be especially relevant as high-frame-rate color Doppler is subject to more
noise related to clutter signals.

###  4.3 Applications in Quantitative Color Doppler

Once corrected, color Doppler images can be used to visualize and quantify
intracardiac blood flow. As mentioned in the introduction, intraventricular
vector flow mapping (iVFM) is an approach to obtain comprehensive flow
information, from which hemodynamic parameters can be estimated. Using a color
M-mode, it is also possible to estimate the pressure difference between the
apex and the mitral base, which reflects the cardiac filling [27, 28].
However, prior dealiasing is required for this method [28]. To this end, the
approaches outlined in this study could be used with color M-mode images. In a
more ambitious perspective, it would be conceivable to develop neural networks
that can directly infer velocity vector fields or relative pressure fields
from color Doppler images, once properly trained. In this case, the dealiasing
process would be intrinsically integrated into the network. The main
difficulty would lie in obtaining paired input data that provide the reference
values. Simulations combining flows and acoustics could provide a relevant
avenue for this purpose [29, 25].

##  5 Conclusion

While traditional methods are effective for interferometric imaging, they are
limited for color Doppler echocardiography due to high noise generated by
clutter signals in moving tissues surrounding the blood flow. We have
demonstrated that deep learning (DL) techniques can achieve alias-free color
Doppler echocardiography. Our proposed DL methods outperformed the non-DL DeAN
method, with nnU-Net achieving the best performance, followed by PDNet. In
addition, the incorporation of power information and artificial aliasing
augmentation improved the results. The application of DL techniques to color
Doppler echocardiography is a promising approach that could enhance the
clinical utility of this widely used imaging modality.

## References

  * [1] W. A. Zoghbi, D. Adams, R. O. Bonow _et al._ , “Recommendations for noninvasive evaluation of native valvular regurgitation: a report from the american society of echocardiography developed in collaboration with the society for cardiovascular magnetic resonance,” _Journal of the American Society of Echocardiography_ , vol. 30, no. 4, pp. 303–371, 2017.
  * [2] K. Hanumansetty, A. M. Aspari, and D. J. Palamattam, “Ventricular septal defects in two dimensional imaging seen as a single defect on colour flow doppler imaging,” _Annals of Cardiac Anaesthesia_ , vol. 25, no. 2, p. 240, 2022.
  * [3] M. V. Revzin, A. Imanzadeh, C. Menias _et al._ , “Optimizing image quality when evaluating blood flow at doppler us: a tutorial,” _Radiographics_ , vol. 39, no. 5, pp. 1501–1523, 2019.
  * [4] F. Vixège, A. Berod, Y. Sun _et al._ , “Physics-constrained intraventricular vector flow mapping by color Doppler,” _Physics in Medicine & Biology_, vol. 66, no. 24, p. 245019, 2021.
  * [5] F. Vixège, A. Berod, P.-Y. Courand _et al._ , “Full-volume three-component intraventricular vector flow mapping by triplane color Doppler,” _Physics in Medicine and Biology_ , vol. 67, no. 9, 2022.
  * [6] N. Painchaud, N. Duchateau, O. Bernard _et al._ , “Echocardiography segmentation with enforced temporal consistency,” _IEEE Transactions on Medical Imaging_ , vol. 41, no. 10, pp. 2867–2878, 2022.
  * [7] D. Posada, J. Porée, A. Pellissier _et al._ , “Staggered Multiple-PRF Ultrafast Color Doppler,” _IEEE Transactions on Medical Imaging_ , vol. 35, 2016.
  * [8] D. C. Ghiglia and L. A. Romero, “Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods,” _Journal of the Optical Society of America A_ , vol. 11, no. 1, p. 107, 1994.
  * [9] A. P. Shanker and H. Zebker, “Edgelist phase unwrapping algorithm for time series InSAR analysis,” _Journal of the Optical Society of America A_ , vol. 27, no. 3, pp. 605–612, 2010.
  * [10] J. Martinez-Carranza, K. Falaggis, and T. Kozacki, “Fast and accurate phase-unwrapping algorithm based on the transport of intensity equation,” _Applied Optics_ , vol. 56, no. 25, pp. 7079–7088, 2017.
  * [11] S. Muth, S. Dort, I. A. Sebag _et al._ , “Unsupervised dealiasing and denoising of color-doppler data,” _Medical Image Analysis_ , vol. 15, no. 4, pp. 577–588, 2011.
  * [12] H. Nahas, J. S. Au, T. Ishii _et al._ , “A deep learning approach to resolve aliasing artifacts in ultrasound color flow imaging,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 67, no. 12, pp. 2615–2628, 2020.
  * [13] J. Adler and O. Öktem, “Learned primal-dual reconstruction,” _IEEE Transactions on Medical Imaging_ , vol. 37, no. 6, pp. 1322–1332, 2018.
  * [14] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” in _Proc. ICLR_ , 2015.
  * [15] F. Isensee, P. F. Jaeger, S. A. A. Kohl _et al._ , “nnU-Net: a self-configuring method for deep learning-based biomedical image segmentation,” _Nature Methods_ , vol. 18, no. 2, pp. 203–211, 2021\.
  * [16] M. Antonelli, A. Reinke, S. Bakas _et al._ , “The Medical Segmentation Decathlon,” _Nature Communications_ , vol. 13, no. 1, p. 4128, 2022.
  * [17] X. Lin, L. Yu, K.-T. Cheng _et al._ , “Batformer: Towards boundary-aware lightweight transformer for efficient medical image segmentation,” _IEEE Journal of Biomedical and Health Informatics_ , pp. 1–13, 2023.
  * [18] F. Mehregan, F. Tournoux, S. Muth _et al._ , “Doppler vortography: a color Doppler approach for quantification of the intraventricular blood flow vortices,” _Ultrasound in medicine & biology_, vol. 40, no. 1, pp. 210–221, 2014.
  * [19] J. Dong, F. Chen, D. Zhou _et al._ , “Phase unwrapping with graph cuts optimization and dual decomposition acceleration for 3d high-resolution mri data,” _Magnetic resonance in medicine_ , vol. 77, no. 3, pp. 1353–1358, 2017.
  * [20] J. M. Bioucas-Dias and G. Valadao, “Phase unwrapping via graph cuts,” _IEEE Transactions on Image processing_ , vol. 16, no. 3, pp. 698–709, 2007\.
  * [21] B. Hunt, “Matrix formulation of the reconstruction of phase values from phase differences,” _JOSA_ , vol. 69, no. 3, pp. 393–399, 1979.
  * [22] G. Spoorthi, S. Gorthi, and R. K. S. S. Gorthi, “Phasenet: A deep convolutional neural network for two-dimensional phase unwrapping,” _IEEE Signal Processing Letters_ , vol. 26, no. 1, pp. 54–58, 2018.
  * [23] K. Wang, Y. Li, Q. Kemao _et al._ , “One-step robust deep learning phase unwrapping,” _Optics express_ , vol. 27, no. 10, pp. 15 100–15 115, 2019\.
  * [24] H. Zhou, C. Cheng, H. Peng _et al._ , “The phu-net: A robust phase unwrapping method for mri based on deep learning,” _Magnetic Resonance in Medicine_ , vol. 86, no. 6, pp. 3321–3333, 2021.
  * [25] Y. Sun, F. Vixege, K. Faraz _et al._ , “A pipeline for the generation of synthetic cardiac color doppler,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 69, no. 3, pp. 932–941, 2022.
  * [26] J. Faurie, M. Baudet, K. C. Assi _et al._ , “Intracardiac vortex dynamics by high-frame-rate doppler vortography—in vivo comparison with vector flow mapping and 4-d flow mri,” _IEEE transactions on ultrasonics, ferroelectrics, and frequency control_ , vol. 64, no. 2, pp. 424–432, 2016.
  * [27] K. Sørensen, S. Fadnes, L. Mertens _et al._ , “Assessment of early diastolic intraventricular pressure difference in children by blood speckle-tracking echocardiography,” _Journal of the American Society of Echocardiography_ , 2023.
  * [28] A. Hodzic, D. Garcia, E. Saloux _et al._ , “Echocardiographic evidence of left ventricular untwisting-filling interplay,” _Cardiovascular Ultrasound_ , vol. 18, no. 1, pp. 1–12, 2020.
  * [29] S. Shahriari and D. Garcia, “Meshfree simulations of ultrasound vector flow imaging using smoothed particle hydrodynamics,” _Physics in Medicine & Biology_, vol. 63, no. 20, p. 205011, 2018.


Present address: ]Paul Scherrer Institute, Villigen 5232, Switzerland

# Electrical transport signatures of metallic surface state
formation in the strongly-correlated insulator FeSb2

A. G. Eaton  N. J. M. Popiel  Cavendish Laboratory, University of Cambridge,
JJ Thomson Avenue, Cambridge, CB3 0HE, United Kingdom  K.-J. Xu  Department of
Applied Physics and Physics, Stanford University,
Stanford, California 94305, USA  A. J. Hickey  H. Liu  Cavendish Laboratory,
University of Cambridge,
JJ Thomson Avenue, Cambridge, CB3 0HE, United Kingdom  M. Ciomaga Hatnean  [
G. Balakrishnan  Department of Physics, University of Warwick,
Coventry, CV4 7AL, United Kingdom  G. F. Lange  R.-J. Slager  TCM Group,
Cavendish Laboratory, University of Cambridge,
JJ Thomson Avenue, Cambridge, CB3 0HE, United Kingdom  Z.-X. Shen  Department
of Applied Physics and Physics, Stanford University,
Stanford, California 94305, USA  Stanford Institute for Materials and Energy
Sciences,
SLAC National Accelerator Laboratory, Menlo Park, California 94025, USA  S. E.
Sebastian<EMAIL_ADDRESS>Cavendish Laboratory, University of
Cambridge,
JJ Thomson Avenue, Cambridge, CB3 0HE, United Kingdom

###### Abstract

We present local and nonlocal electrical transport measurements of the
correlated insulator FeSb2. By employing wiring configurations that delineate
between bulk- and surface-dominated conduction, we reveal the formation of a
metallic surface state in FeSb2 for temperatures $\lessapprox 5$ K. This
result is corroborated by an angular rotation study of this material’s
magnetotransport, which also shows signatures of the transition from bulk- to
surface-dominated conduction over the same temperature interval as the
local/nonlocal transport divergence. Notable similarities with the topological
Kondo insulator candidate SmB6 are discussed.

Exploring the manifestation of strong electron-electron correlation effects in
the bulk of electrically insulating materials has led to several significant
advancements in our understanding of quantum matter [1, 2, 3]. Recently, the
exploration of topological insulators – combining the juxtaposition of a
conducting metallic surface state atop a non-correlated insulating bulk – has
also yielded notable results [4, 5]. The interplay of these two phenomena –
metallic surface states in combination with strong bulk interactions – may
thereby prove to be fertile ground for exploring novel physical behavior.

Iron diantimonide (FeSb2) is a strongly correlated $d$-electron insulator [6]
that crystallizes in the orthorhombic P${}_{\text{nnm}}$ crystal structure,
and exhibits several similarities in common with $f$-electron Kondo insulators
[7]. Evidence of electron-electron interactions in this material has been
observed in a variety of physical properties, including electrical and thermal
transport [8], magnetic susceptibility [9], optical conductivity [10],
specific heat [11], and x-ray absorption spectroscopy [12]. At low
temperatures $T\approx 10$ K a colossal Seebeck effect has been reported, of
$S\approx-45$ mV/K [13, 14]. While the microscopic origin of this effect is
not well understood, it has been conjectured to arise from strongly
renormalized carrier masses, the transport of which may be assisted through
the phonon-drag effect [15, 16, 17].

Further effects of the strongly interacting ground state of this material have
been uncovered by introducing small doping quantities of Te [18], Co [19], and
Sn [20] which reveal a proximate transition to a heavy-fermion metallic state.
A similar doping-induced evolution, from strongly correlated insulator to
disordered Fermi liquid, has previously been observed in the related compound
iron monosilicide (FeSi) [21, 22, 23]. FeSi also exhibits several
characteristics of Kondo insulator-like strong correlations [7, 24, 25].
Recent electrical transport measurements on FeSi have indicated the presence
of conductive metallic surface states [26] reminiscent of previous findings on
samarium hexaboride (SmB6) [27, 28], a 3D topological Kondo insulator
candidate [29, 30, 31]. While recent angle-resolved photoemission spectroscopy
(ARPES) measurements of FeSb2 have indicated the presence of metallic surface
states on this material as well [32], complementary electrical transport
studies are required in order to probe the interplay between bulk and surface
effects.

In this Letter, we present local and nonlocal electrical transport
measurements on the correlated insulator FeSb2. Samples were grown by the
chemical vapor transport (CVT) method, from starting materials of 99.995%at Fe
and 99.999%at Sb. Measurements were performed in a Quantum Design Physical
Properties Measurement System (PPMS) over a temperature range of 1.7–300 K and
up to applied magnetic field strengths of 14 T. Samples were contacted with 25
$\upmu$m diameter gold wire using DuPont 4929N silver paint thinned with
2-butoxyethyl acetate to obtain a suitable consistency with low contact
resistance. Some samples also had these connections spot-welded, with no
discernible differences between these two preparatory techniques. All contact
resistances were measured at room temperature using two-point followed by
four-point measurements, indicating combined lead and contact resistances
$\sim 1$ $\Omega$ in all cases.

Figure 1:  Electrical transport anomalies. Resistivity, $\rho$, as a function
of temperature for (a) SmB6, and (b) FeSb2, with the Hall coefficient,
-$R_{H}$, plotted on the same logarithmic temperature scale for each in (c)
and (d). Both materials show the anomalous formation of a resistance plateau
for low temperatures $T\lessapprox 5$ K, which in the case of SmB6 has
previously been associated with the development of the low-temperature surface
state [30, 31]. Upper insets show resistivity against inverse temperature,
with arrows indicating the location of the peak in -$R_{\text{H}}$, which
occurs for both materials in close proximity to the resistivity plateau. Red
lines show regions fitted to extract activated charge gap energies
$\varDelta$. FeSb2 exhibits two gaps, with the lower temperature one,
$\varDelta_{2}$, being of very similar magnitude to the sole gap of SmB6.
Lower insets show the respective crystal structures [33, 34, 35].

Figure 1 compares the temperature dependence of electrical transport for SmB6
and FeSb2. In general, the resistivity of a Kondo insulator should increase
logarithmically as T $\xrightarrow{~{}}0$ K [36]; however, both SmB6 [28, 37]
and FeSb2 exhibit low-temperature resistance plateaux, a hallmark of the
development of a surface state, for $T$ $\lessapprox$ 4 K. All sample batches
examined in this study exhibited a low temperature resistance plateau (see
Supplementary Materials). The formation of this plateau is found to be
accompanied by a steep peak in the Hall coefficient, -$R_{\text{H}}$.

Fitting the FeSb2 resistivity data in Fig. 1b to a doubly gapped Arrhenius
activation model, of the form
$\rho\propto\exp(\nicefrac{{\varDelta}}{{k_{\text{B}}T}})$, yields a high
temperature indirect bandgap, $\varDelta_{1}$, of $\approx$ 12 meV for $T>$ 30
K, comparable with prior reports [7]. However, at lower temperatures another
gap – $\varDelta_{2}\approx$ 3.7 meV – opens up, of comparable magnitude to
the 4.0 meV gap in SmB6 attributed to Kondo-induced band hybridisation [31].
Very similar observations, of low-temperature resistivity saturation
concomitant with a dramatic reduction in $R_{\text{H}}$, have previously been
attributed with topologically protected surface state formation in the non-
correlated topological insulators ${\text{Bi}}_{2}{\text{Te}}_{2}\text{Se}$
[38] and Bi1.08Sn0.02Sb0.9Te2S [39]. The notable similarities between these
materials and FeSb2 are strongly suggestive of similar microscopic mechanisms
at play in this system.

Figure 2:  Magnetotransport evolution from bulk- to surface-dominated
conduction on cooling. All traces performed with a constant dc current of 25
$\mu$A for $T\leq 5$ K and 100 $\mu$A for $T>5$ K. (a) Magnetoresistance of
SmB6 at successive temperatures as a function of rotation angle, $\theta$. In
all panels, the right hand side schematic gives the applied current direction,
j, and rotational plane, while purple curves in each plot give symmetrised
resistance for applied magnetic fields of $\pm$ 14 T, to exclude any
transverse contributions. At 15 K, away from the resistivity plateau, SmB6
exhibits $C_{4}$ rotational symmetry, whereas at 2 K, on the plateau, this has
reduced to $C_{2}$ symmetry. (b), Similar temperature progression on an FeSb2
crystal as for SmB6 in panel (a), clearly exhibiting a very similar peanut-
shaped trace at 2 K.

To further investigate the potential surface state origin of the low-
temperature resistivity plateau of FeSb2, a rotational magnetotransport study
was performed. Figure 2 shows the angular-dependent magnetoresistance (ADMR)
of two FeSb2 samples and one SmB6 sample as a function of rotation angle
$\theta$, at various temperatures in applied magnetic fields of $\pm$ 14 T.
Samples with large surface to volume ratios were chosen, with
$\theta=0\degree$ corresponding to field applied normal to the two largest
(parallel) faces, the top one of which was contacted with gold wires for
standard four-point electrical transport measurements as indicated in the
schematics in Fig. 2. For $T\geq$ 15 K, well away from the plateau, both
materials’ longitudinal magnetoresistances exhibit almost isotropic rotational
symmetry, as expected for simple bulk conduction. However, on cooling to
$T\leq$ 3 K, within the plateau region, the rotational magnetoconductance
profiles of both SmB6 and FeSb2 undergo a pronounced evolution. In Fig. 2A at
5 K the studied SmB6 sample exhibits nodes (antinodes) in the magnitude of the
magnetoresistance at 0$\degree$ and 180$\degree$ (90$\degree$ and
270$\degree$); however, at 2 K this has swapped to give nodes (antinodes) at
90$\degree$ and 270$\degree$ (0$\degree$ and 180$\degree$). Comparatively, the
FeSb2 sample in Fig. 2b evolves similarly – with a stark difference between
the 5 K and 4 K temperature points, whereby the nodes at 90$\degree$ and
270$\degree$ at 5 K have shifted to 0$\degree$ and 180$\degree$ at 4 K.

Interestingly, this evolution goes through a drastic change around the
temperature of the Hall effect peak (Fig. 2), with a convolution of 3D and 2D
transport signatures observed in the ‘butterfly’ shape of the 5 K trace in
Fig. 2b – similar behaviour has previously been reported in studies of Weyl
semi-metals [40]. At the lowest temperatures, a distinct peanut-shaped trace
is observed in the rotational magnetoresistance of FeSb2, in marked similarity
with prior SmB6 reports [41, 42].

This symmetry evolution fits well within the surface state model, whereby at
the lowest temperatures electrical conduction is dominated by these surface
states, as their conductivity is much higher than that of the insulating bulk.
Hence, for the samples measured here with two large, dominant flat faces (as
depicted in Fig. 2), we can conclude that a model comprising the presence of
surface-dominated conduction for $T\lessapprox 6$ K, and bulk-dominated
conduction at elevated temperatures, very well captures the symmetry evolution
observed in the ADMR data as a function of temperature for FeSb2 (and for
SmB6).

Figure 3:  Bulk- and surface-dominated conduction. Finite element simulations
of the current density distributions for electrical transport through (a),
bulk-dominated channels, and (b), for surface-dominated conduction. Figure 4:
Local and nonlocal 4-wire electrical transport measurements performed on an
FeSb2 single crystal. Current and voltage leads were attached as per the
configurations indicated in Fig. 3. At high temperatures the nonlocal
resistance is greater than the local resistance, as expected for simple bulk
conduction. However, at low temperatures the converse is true, with the local
resistance exceeding that of the nonlocal resistance. This inversion is
characteristic of the formation of a conductive surface state.

To further probe the character of low-temperature electrical conduction in
FeSb2, a nonlocal resistance study was undertaken. Previously, an inversion
between local and nonlocal electrical transport measurements upon cooling to
the resistivity plateau was taken as a key piece of evidence in support of
metallic surface state formation in SmB6 [31, 27]. In the general case of
electrical conduction through a section of bulk 3D conductor, electrical
transport is expected to take place predominantly through the bulk of the
material, as indicated by the current density distribution depicted in Fig.
3a. However, for a sample cooled to $T$ $\lessapprox$ 4 K, the surface state
model predicts that band hybridisation will cause these bulk conduction
channels to become strongly gapped out as the surface state forms. This
results in electrical transport being constrained almost entirely to the
metallic surface of a sample (Fig. 3b). By performing four-wire resistance
measurements of local and non-local configurations, as shown by Figure 3, the
applicability of this bulk/surface model to FeSb2 may be investigated, and
compared to the excellent correspondence previously found for such
measurements in SmB6 [31, 27].

Fig. 4 shows the measured electrical resistance of an FeSb2 crystal for local
(blue curve) and nonlocal (red curve) wiring configurations. At elevated
temperatures ($T\gg$ 4 K) the resistance of the non-local transport channel is
greater than the local channel, as expected by simple geometric resistivity
arguments for bulk conduction. However, for $T$ $\lessapprox$ 5 K the reverse
is true, with $R_{\text{local}}>R_{\text{non-local}}$. This crossover occurs
just before the resistivity plateau, over the same temperature range as the
Hall effect peak in Fig. 1. Such a remarkable observation would be difficult
to reconcile within a bulk impurity-band picture, or even with an elaborate
low-temperature highly-anisotropic conductivity model. However, these
observations may be easily explained by the surface state picture as applied
to very similar observations in SmB6 [31, 27], and thus provide compelling
evidence that the same phenomenon is occurring in FeSb2. We note that these
measurements are in strong agreement with a recent corbino disc study [43].

Furthermore, the fact that the local/non-local resistance inversion should
occur over the same small temperature range (approximately 4 to 6 K) as the
peak in the Hall coefficient, the opening of the resistance plateau, and the
magnetoresistive symmetry evolution, each in concurrence with previous reports
on SmB6, gives strong confidence that these observations all point to the same
correlation-driven surface state mechanism, previously investigated in the
case of SmB6, also being present in FeSb2.

To investigate the possibility of non-trivial topology manifesting the
observed surface states in FeSb2, we performed DFT calculations to evaluate
the crystalline symmetry (see Supplemental Materials). However, from these
calculations we find that $z_{4}=0$, indicating that the bandstructure of
FeSb2 is topologically trivial. We note that our DFT calculations are limited
to assuming a nonmagnetic groundstate as previously reported [44] – however,
the presence of localized Fe 3d states may necessitate a more detailed
calculation, fully accounting for the strong electronic correlations, to fully
model the true bandstructure. Regardless of the possibility of a topologically
trivial or nontrivial origin, the metallic surface state in FeSb2 represents a
fascinating platform on which to explore the interplay between strongly
correlated physics and reduced dimensionality.

It should be remarked that a clear-cut consensus on SmB6 being a topological
Kondo insulator is yet to reach unanimous agreement throughout the community;
although, recently this appears to be nearing [31]. That these measurements on
FeSb2 presented here yield a plethora of striking similarities with SmB6 opens
a number of pertinent questions, especially regarding the role of $f$-electron
physics – previously thought to be responsible for driving the surface state
mechanism in SmB6, but obviously not present in the case of FeSb2.

Studies on SmB6 into what physical features are fundamental properties
required for the manifestation of the (putatively topological) surface state,
and what features are not relevant, have been severely hampered by drastic
differences between sample preparation methods [31]. Samples grown in aluminum
flux are reasonably easy to produce, however the presence of aluminum
inclusions, inhomogeneous growth rates of different crystalline facets, and
contamination from the melt container, have led to notable inconsistencies in
observed properties between different research groups [45]. The presence of
small quantities of rare-earth impurities, such as gadolinium, can have a
pronounced effect on the physical properties [46] at the low temperatures
where the surface state forms, and the purification of elemental lanthanides
from extraneous impurities of other rare-earth contaminants is notoriously
difficult.

Higher purity SmB6 samples may be obtained via the floating zone technique
[47], as evidenced by these crystals’ superior inverse residual resistance
ratios, magnetisation, quantum oscillation and thermal transport properties
[48]. Samples grown by this method have also exhibited some of the most
compelling ARPES evidence for a topologically non-trivial nature of the
surface state [49]. However, this method is considerably more difficult to
perform than flux growth, with the high melting point of SmB6 requiring the
use of mirror furnaces to reach temperatures approaching $\sim$ 3,000 K [47].
Furthermore, the attainment of such high temperatures increases the likelihood
of introducing samarium vacancies, which have also been reported to be a
significant variable factor in the correlated electronic properties of this
material [50]. These variations in physical properties between differing
growth batches and sample preparation methodologies of SmB6 has hampered the
identification of what constituent properties are prerequisite to the
development of the correlation-induced surface state.

FeSb2, by contrast, has a much lower melting point of just over 1,000 K. This
makes the acquisition of high-purity samples much more accessible through use
of the CVT method, such as for the crystals investigated in this study. The
CVT method, like the floating zone method, is immune to extraneous
substitutive impurities from a crucible origin, or from flux inclusions;
importantly, the considerably lower temperatures involved significantly
minimize the probability of introducing vacancies, dislocations, or other
sources of crystalline disorder that higher temperature procedures suffer
from. Hence, obtaining excellent purity FeSb2 crystals is much more realizable
for a broader number of experimental research groups than is the attainment of
similarly high quality SmB6.

Therefore, given the striking similitude in low-temperature surface state
properties between FeSb2 and SmB6, an exciting possibility lies in the
identification of FeSb2 as a model $d$-electron analogue of SmB6, with
extremely similar low-temperature behaviour, but with chemical properties that
make the attainment of pristine high-purity specimens significantly easier to
realize.

In conclusion, we have experimentally observed evidence of low-temperature
surface-dominated conduction in the strongly correlated $d$-electron insulator
FeSb2 through measurements of local and nonlocal electrical transport. This
interpretation is corroborated by an angular rotation study of this material’s
magnetotransport, which shows an associated evolution from bulk- to surface-
dominated symmetry on cooling to temperatures $\lessapprox 5$ K. Complementary
measurements, such as scanning tunneling microscopy and spin-polarised ARPES,
would be instructive in resolving any possible topological character of these
surface states, and may reveal further manifestations of exotic correlation-
driven effects previously observed in similar systems.

## References

  * Mott [1949] N. F. Mott, The Basis of the Electron Theory of Metals, with Special Reference to the Transition Metals, Proc. Phys. Soc., London, Sect. A 62, 416 (1949).
  * Anderson [1958] P. W. Anderson, Absence of Diffusion in Certain Random Lattices, Phys. Rev. 109, 1492 (1958).
  * Kondo [1964] J. Kondo, Resistance Minimum in Dilute Magnetic Alloys, Progress of Theoretical Physics 32, 37 (1964).
  * Moore [2010] J. E. Moore, The birth of topological insulators., Nature 464, 194 (2010).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Sun  _et al._ [2011] P. Sun, M. Søndergaard, B. Iversen, and F. Steglich, Strong electron correlations in FeSb2, Ann. Phys. 523, 612 (2011).
  * Tomczak [2018] J. M. Tomczak, Thermoelectricity in correlated narrow-gap semiconductors, J. Phys. Condens. Matter 30, 183001 (2018).
  * Sun  _et al._ [2009a] P. Sun, N. Oeschler, S. Johnsen, B. B. Iversen, and F. Steglich, ${\text{FeSb}}_{2}$: Prototype of huge electron-diffusion thermoelectricity, Phys. Rev. B 79, 153308 (2009a).
  * Petrovic  _et al._ [2005] C. Petrovic, Y. Lee, T. Vogt, N. D. Lazarov, S. L. Bud’ko, and P. C. Canfield, Kondo insulator description of spin state transition in $\mathrm{Fe}{\mathrm{Sb}}_{2}$, Phys. Rev. B 72, 045103 (2005).
  * Herzog  _et al._ [2010] A. Herzog, M. Marutzky, J. Sichelschmidt, F. Steglich, S. Kimura, S. Johnsen, and B. B. Iversen, Strong electron correlations in FeSb2: An optical investigation and comparison with RuSb2, Phys. Rev. B 82, 245205 (2010).
  * Sun  _et al._ [2009b] P. Sun, N. Oeschler, S. Johnsen, B. B. Iversen, and F. Steglich, Huge Thermoelectric Power Factor: FeSb2 versus FeAs2 and RuSb2, Appl. Phys. Express 2, 091102 (2009b).
  * [12] P. Sun, N. Hollmann, W. H. Zhi, W. Schnelle, M. Baenitz, S. Johnsen, M. Søndergaard, B. Iversen, C. C. Cezar, N. B. Brookes, N. H. Tjeng, and F. Steglich,  _Magnetism and Spin States of FeSb 2_ (Max Planck Institute for Chemical Physics of Solids Scientific Report, Dresden, 2011) pp. 56–58, https://www.cpfs.mpg.de/en/publications/scientific_report .
  * Bentien  _et al._ [2007] A. Bentien, S. Johnsen, G. K. Madsen, B. B. Iversen, and F. Steglich, Colossal Seebeck coefficient in strongly correlated semiconductor FeSb 2, EPL 80, 17008 (2007).
  * Sun  _et al._ [2010] P. Sun, N. Oeschler, S. Johnsen, B. B. Iversen, and F. Steglich, Narrow band gap and enhanced thermoelectricity in FeSb2, Dalton Trans. 39, 1012 (2010).
  * Hong  _et al._ [2013] S. Hong, P. Ghaemi, J. E. Moore, and P. W. Phillips, Tuning thermoelectric power factor by crystal-field and spin-orbit couplings in Kondo-lattice materials, Phys. Rev. B 88, 075118 (2013).
  * Homes  _et al._ [2018] C. C. Homes, Q. Du, C. Petrovic, W. H. Brito, S. Choi, and G. Kotliar, Unusual electronic and vibrational properties in the colossal thermopower material FeSb2, Scientific Reports 8, 1 (2018).
  * Takahashi  _et al._ [2016] H. Takahashi, R. Okazaki, S. Ishiwata, H. Taniguchi, A. Okutani, M. Hagiwara, and I. Terasaki, Colossal Seebeck effect enhanced by quasi-ballistic phonons dragging massive electrons in FeSb2, Nat. Commun. 7, 12732 (2016).
  * Hu  _et al._ [2012] R. Hu, K. Wang, H. Ryu, H. Lei, E. S. Choi, M. Uhlarz, J. Wosnitza, and C. Petrovic, Electronic Griffiths Phase in the Te-Doped Semiconductor ${\mathrm{FeSb}}_{2}$, Phys. Rev. Lett. 109, 256401 (2012).
  * Kassem  _et al._ [2020] M. A. Kassem, Y. Tabata, T. Waki, and H. Nakamura, Crystal growth and metallic ferromagnetism induced by electron doping in FeSb2 (2020), arXiv:2011.10343 .
  * Bentien  _et al._ [2006] A. Bentien, G. K. H. Madsen, S. Johnsen, and B. B. Iversen, Experimental and theoretical investigations of strongly correlated ${\mathrm{FeSb}}_{2-x}{\mathrm{Sn}}_{x}$, Phys. Rev. B 74, 205105 (2006).
  * DiTusa  _et al._ [1997] J. F. DiTusa, K. Friemelt, E. Bucher, G. Aeppli, and A. P. Ramirez, Metal-Insulator Transitions in the Kondo Insulator FeSi and Classic Semiconductors Are Similar, Phys. Rev. Lett. 78, 2831 (1997).
  * DiTusa  _et al._ [1998] J. F. DiTusa, K. Friemelt, E. Bucher, G. Aeppli, and A. P. Ramirez, Heavy fermion metal–Kondo insulator transition in ${\mathrm{FeSi}}_{1-x}{\mathrm{Al}}_{x}$, Phys. Rev. B 58, 10288 (1998).
  * Manyala  _et al._ [2008] N. Manyala, J. F. DiTusa, G. Aeppli, and A. P. Ramirez, Doping a semiconductor to create an unconventional metal, Nature 454, 976 (2008).
  * Schlesinger  _et al._ [1993] Z. Schlesinger, Z. Fisk, H.-T. Zhang, M. B. Maple, J. DiTusa, and G. Aeppli, Unconventional charge gap formation in FeSi, Phys. Rev. Lett. 71, 1748 (1993).
  * Mandrus  _et al._ [1995] D. Mandrus, J. L. Sarrao, A. Migliori, J. D. Thompson, and Z. Fisk, Thermodynamics of FeSi, Phys. Rev. B 51, 4763 (1995).
  * Fang  _et al._ [2018] Y. Fang, S. Ran, W. Xie, S. Wang, Y. S. Meng, and M. B. Maple, Evidence for a conducting surface ground state in high-quality single crystalline FeSi, Proc. Natl. Acad. Sci. U.S.A. 115, 8558 (2018).
  * Wolgast  _et al._ [2013] S. Wolgast, Ç. Kurdak, K. Sun, J. W. Allen, D.-J. Kim, and Z. Fisk, Low-temperature surface conduction in the Kondo insulator SmB6, Phys. Rev. B 88, 180405(R) (2013).
  * Kim  _et al._ [2013] D. J. Kim, S. Thomas, T. Grant, J. Botimer, Z. Fisk, and J. Xia, Surface Hall Effect and Nonlocal Transport in SmB6: Evidence for Surface Conduction, Scientific Reports 3, 1 (2013).
  * Dzero  _et al._ [2010] M. Dzero, J. Xia, V. Galitski, and P. Coleman, Topological Kondo Insulators, Phys. Rev. Lett. 104, 106408 (2010).
  * Dzero  _et al._ [2016] M. Dzero, J. Xia, V. Galitski, and P. Coleman, Topological Kondo Insulators, Ann. Rev. Cond. Matt. Phys. 7, 249 (2016).
  * Li  _et al._ [2020] L. Li, K. Sun, C. Kurdak, and J. W. Allen, Emergent mystery in the Kondo insulator samarium hexaboride, Nat. Rev. Phys. 2, 463 (2020).
  * Xu  _et al._ [2020] K.-J. Xu, S.-D. Chen, Y. He, J. He, S. Tang, C. Jia, E. Yue Ma, S.-K. Mo, D. Lu, M. Hashimoto, T. P. Devereaux, and Z.-X. Shen, Metallic surface states in a correlated d-electron topological Kondo insulator candidate FeSb2, Proc. Natl. Acad. Sci. U.S.A. 117, 15409 (2020).
  * Momma and Izumi [2011] K. Momma and F. Izumi, VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data, Journal of Applied Crystallography 44, 1272 (2011).
  * spr [2016a] SmB6 Crystal Structure Datasheet (2016a), Springer-Verlag Berlin Heidelberg & Material Phases Data System (MPDS), Switzerland & National Institute for Materials Science (NIMS), Japan.
  * spr [2016b] FeSb2 Crystal Structure Datasheet (2016b), Springer-Verlag Berlin Heidelberg & Material Phases Data System (MPDS), Switzerland & National Institute for Materials Science (NIMS), Japan.
  * Coleman [2015] P. Coleman,  _Introduction to Many-Body Physics_ (Cambridge Univ. Press, 2015) pp. 582–655.
  * Kim  _et al._ [2014] D. J. Kim, J. Xia, and Z. Fisk, Topological surface state in the Kondo insulator samarium hexaboride, Nature Mater. 13, 466 (2014).
  * Ren  _et al._ [2010] Z. Ren, A. A. Taskin, S. Sasaki, K. Segawa, and Y. Ando, Large bulk resistivity and surface quantum oscillations in the topological insulator Bi2Te2Se, Phys. Rev. B 82, 241306 (2010).
  * Kushwaha  _et al._ [2016] S. K. Kushwaha, I. Pletikosic, T. Liang, A. Gyenis, S. H. Lapidus, Y. Tian, H. Zhao, K. S. Burch, J. Lin, W. Wang, H. Ji, A. V. Fedorov, A. Yazdani, N. P. Ong, T. Valla, and R. J. Cava, Sn-doped Bi1.1Sb0.9Te2S bulk crystal topological insulator with excellent properties, Nat. Commun. 7, 11456 (2016).
  * Ali  _et al._ [2016] M. N. Ali, L. M. Schoop, C. Garg, J. M. Lippmann, E. Lara, B. Lotsch, and S. S. P. Parkin, Butterfly magnetoresistance, quasi-2D Dirac Fermi surface and topological phase transition in ZrSiS, Sci. Adv. 2, 12 (2016).
  * Chen  _et al._ [2015] F. Chen, C. Shang, Z. Jin, D. Zhao, Y. P. Wu, Z. J. Xiang, Z. C. Xia, A. F. Wang, X. G. Luo, T. Wu, and X. H. Chen, Magnetoresistance evidence of a surface state and a field-dependent insulating state in the Kondo insulator ${\mathrm{SmB}}_{6}$, Phys. Rev. B 91, 205133 (2015).
  * Yue  _et al._ [2015] Z. Yue, X. Wang, D. Wang, J. Wang, D. Culcer, and S. Dou, Crossover of Magnetoresistance from Fourfold to Twofold Symmetry in SmB6 Single Crystal, a Topological Kondo Insulator, J. Phys. Soc. Japan 84, 044717 (2015).
  * Eo  _et al._ [2023] Y. S. Eo, K. Avers, J. A. Horn, H. Yoon, S. R. Saha, A. Suarez, M. S. Fuhrer, and J. Paglione, Extraordinary bulk-insulating behavior in the strongly correlated materials FeSi and FeSb2, Applied Physics Letters 122, 233102 (2023).
  * Petrovic  _et al._ [2003] C. Petrovic, J. W. Kim, S. L. Bud’ko, A. I. Goldman, P. C. Canfield, W. Choe, and G. J. Miller, Anisotropy and large magnetoresistance in the narrow-gap semiconductor FeSb2, Physical Review B - Condensed Matter and Materials Physics 67, 155205 (2003).
  * Phelan  _et al._ [2016] W. A. Phelan, S. M. Koohpayeh, P. Cottingham, J. A. Tutmaher, J. C. Leiner, M. D. Lumsden, C. M. Lavelle, X. P. Wang, C. Hoffmann, M. A. Siegler, N. Haldolaarachchige, D. P. Young, and T. M. McQueen, On the Chemistry and Physical Properties of Flux and Floating Zone Grown SmB6 Single Crystals, Sci. Rep. 6, 20860 (2016).
  * Fuhrman  _et al._ [2018] W. T. Fuhrman, J. R. Chamorro, P. A. Alekseev, J. M. Mignot, T. Keller, J. A. Rodriguez-Rivera, Y. Qiu, P. Nikolić, T. M. McQueen, and C. L. Broholm, Screened moments and extrinsic in-gap states in samarium hexaboride, Nat. Commun. 9, 1539 (2018).
  * Hatnean  _et al._ [2013] M. C. Hatnean, M. R. Lees, D. M. Paul, and G. Balakrishnan, Large, high quality single-crystals of the new Topological Kondo Insulator, SmB6, Sci. Rep. 3, 3071 (2013).
  * Hartstein  _et al._ [2020] M. Hartstein, H. Liu, Y. T. Hsu, B. S. Tan, M. Ciomaga Hatnean, G. Balakrishnan, and S. E. Sebastian, Intrinsic Bulk Quantum Oscillations in a Bulk Unconventional Insulator SmB6, iScience 23, 101632 (2020).
  * Ohtsubo  _et al._ [2019] Y. Ohtsubo, Y. Yamashita, K. Hagiwara, S. ichiro Ideta, K. Tanaka, R. Yukawa, K. Horiba, H. Kumigashira, K. Miyamoto, T. Okuda, W. Hirano, F. Iga, and S. ichi Kimura, Non-trivial surface states of samarium hexaboride at the (111) surface, Nat. Commun. 10, 2298 (2019).
  * Eo, Yun Suk and Rakoski, Alexa and Sinha, Shriya and Mihaliov, Dmitri and Fuhrman, Wesley T. and Saha, Shanta R. and Rosa, Priscila F. S. and Fisk, Zachary and Hatnean, Monica Ciomaga and Balakrishnan, Geetha and Chamorro, Juan R. and Phelan, W. Adam and Koohpayeh, Seyed M. and McQueen, Tyrel M. and Kang, Boyoun and Song, Myung-suk and Cho, Beongki and Fuhrer, Michael S. and Paglione, Johnpierre and Kurdak, Cagliyan [2021] Eo, Yun Suk and Rakoski, Alexa and Sinha, Shriya and Mihaliov, Dmitri and Fuhrman, Wesley T. and Saha, Shanta R. and Rosa, Priscila F. S. and Fisk, Zachary and Hatnean, Monica Ciomaga and Balakrishnan, Geetha and Chamorro, Juan R. and Phelan, W. Adam and Koohpayeh, Seyed M. and McQueen, Tyrel M. and Kang, Boyoun and Song, Myung-suk and Cho, Beongki and Fuhrer, Michael S. and Paglione, Johnpierre and Kurdak, Cagliyan, Bulk transport paths through defects in floating zone and Al flux grown ${\mathrm{SmB}}_{6}$, Phys. Rev. Materials 5, 055001 (2021).